Harmonic balance solutions and solver

SDT-nlsim Contents Functions

1.7 Harmonic balance solutions and solver

1.7.1 Time/frequency representation of solutions/loads

The solver is meant to support a generic representation of time dependence. The generalized degrees of freedom considered here are components of a vector noted Z and correspond to amplitudes multiplying functions of both space and time. For a scalar function of space (single spatial DOF), multiple harmonic DOF f_k and time varying shape functions h_k(t) fully describe the time dependence

(1.121)

In the specific case of harmonic balance [20], the time dependence of a given DOF is written as

(1.122)

where spatial DOFs are indexed by q and temporal DOFs are indexed by k but need also a representation of the harmonic n. Thus the ordering of generalized space/time DOFs is

(1.123)

with a spatial dependence q (degree of freedom in the time domain equation(1.1)) and a temporal dependence k. Note that in the notation above, the individual components are assumed real, but the problem can also be written in complex form (1.135). As always in SDT, it is expected that the ordering of Z could be changed. Thus a list of DOFs is needed to specify the meaning of each component of the vector. This list is given in the form

hdof={ % List of active DOFs
  1.01   'c0'  % Node1.DOF1(x), constant(B0)
 10.02   's1'  % Node10.DOF2(y), first harmonic (A1)
 10.02   'c1'  % Node10.DOF2(y), first harmonic (B1)
};

where the first column specifies the spatial DOF and the label in the second column the temporal component. hbm_solve harm commands provide utilities to manipulate spatio-temporal DOF definitions.

The definition of a solution is thus made using a data structure with fields

.def matrix with columns being Z_qk values in given configuration (frequency dependence will be obtained with multiple columns).
.hdof cell array giving the meaning of spatio-temporal amplitudes in .def rows (see equation (1.122))
TR since, in general harmonic balance is performed using a reduced model, restitution of the full DOFs is based on (??) with the reduction basis stored in def.TR.

Building of the time response of a given DOF q(t) should, using the convention of summed indices, actually be written as

(1.124)

with H_lj = H_l(t_j) the l^th harmonic function at time t_j. Times are assumed to be sampled into the regular time interval [t₁ : t_N ], so that H can be written in this case as a matrix

(1.125)

In further equations and in the code, the notation Hkt(H_kt) is used even though in the real form (with opt.Opt.complex=0) there are two l lines for each harmonic k (except for harmonic 0). Note also that when using harmonic fractions (for example in engines it is usual to use N/2 harmonics to represent a period over two revolutions), you should declare opt.Opt.nu=2 and the .harm field will contain k=1/2,1, ....

An advantage of the retained definition is that it is not necessary to define all harmonics. Furthermore problems with sub-harmonics of the excitation frequency can also be considered and only require development of appropriate label handling methods to define sub-harmonic functions h_k(t).

The inverse transform is the way to obtain the harmonic amplitudes from a time vector and is associated with the Htk linear operator

(1.126)

It is thus verified that H_ktH_tk=[I]_{Nk× Nk}.

Some of the literature considers linear DOFs (including modal DOFs and some physical DOFs), non-linear DOFs (internal or physical). SDT-HBM does not ask the user for the nature of DOFs since non-linear DOFs can actually be deduced from the consideration of DOFs present in non-linear observation matrices and the notion was not found to be needed.

The spatial DOF nature (physical, modal, internal) is also not necessary. It is however useful to have automated procedures to assign an identifier for every spatial DOF. It is thus expected that generalized (modal) DOFs be assigned a node number at the time of superelement creation. To avoid viewing mistakes DOF associated with modes or internal states are affected to DOF 99 (of the form NodeId.99).

Finally, some non-linearities use internal states. As the HBM solver needs explicit access to these states, the hbm_solve InitHBM command performs a pre-processing step that affects a node number for each internal state so that they are present in q.

1.7.2 Harmonic balance equation (real DOF)

The base iterations are associated with a non-linear least squares minimization problem of the form

(1.127)

The contents of matrices [A(ω)], [b], and [b_nl] depend on the choice of time functions in (1.124).

With the states described in (1.122) using (1.125), the rows of the residue matrix correspond to the work equilibrium equation (1.1) integrated over a period. Thus the harmonic 0 (constant term) is given by

(1.128)

The harmonic k leads to a sine (respectively cosine) contribution corresponding to an integral over the N time points of a period

(1.129)

Assuming states ordered with sine contributions first followed by cosine, the generalized Jacobian matrix ∂ R/∂ Z is composed of a linear part

(1.130)

and a series of non-linear contributions

(1.131)

The linear part of A_H in (1.130) is skew-symmetric. The negative sign is applied to the lines corresponding to the sine contributions. Using the notation (1.122), the negative sign is then located on the upper triangular part.

In the implementations derived from [21], the non-linear contributions are computed by finite differences (or possibly analytically although this is not yet implemented). And the inner loop iterations of the non-linear solver seeking the solution of (1.127) are of the form Δ Z=[J]⁻¹R.

In direct frequency solvers by SDTools, one seeks to obtain the ideal single step convergence, where a sequant inter-harmonic stiffness [K_єH] is found such that

(1.132)

The two have the notable major difference that in general [K_єH] ≠ ∂ σ/∂ є. In other words secant and tangent stiffness are different matrices.

In the case of a scalar strain є and a linear complex stiffness of the form k(1+iη). The inter-harmonic coupling is block diagonal and of the form

(1.133)

since −iσ_s+σ_c=k(1+iη)(iє_s+є_c)=ik(є_s + η є_c)+k(є_c−ηє_s).

Another form considers an extended Z with frequency stored as an additional component (this is achieved using opt.Opt.fvar=1)

(1.134)

In this case a last column ∂ R/∂ ω is added to the Jacobian. For the last row, the ideal would be to compute ∂ ω/∂ Z but the evaluation of this quantity is quite difficult, so that arc length techniques are used for continuation.

1.7.3 Complex formulation and equivalent stiffness

To clarify the notion of equivalent stiffness found at the harmonic balance solution, the real harmonic balance states (1.122) can be written in complex form expressing the time response as

(1.135)

where spatial DOFs are indexed by q and temporal DOFs are indexed by k. Thus the ordering of generalized space/time complex DOFs

(1.136)

The harmonic function is then complex,

(1.137)

The usual harmonic balance using Fourier coefficients, one has the equivalence

(1.138)

as e^{i k ω t}=cos(kω t) + i sin(kω t), the equivalence between (1.135) and (1.122) is indeed

(1.139)

For each harmonic k the mean equilibrium over a period (1.129) leads to an equation of the form

(1.140)

where

(1.141)

Since b is a constant matrix, the harmonic contribution of a non-linear stress can be computed at the gauss point level using complex notation

(1.142)

For a linear spring with a loss factor and harmonic strains u_k

(1.143)

Assuming a scalar stress relation, the non-linear harmonic problem (1.140) is thus equivalent to a linear harmonic problem

(1.144)

with the equivalent complex stiffness given for each Gauss point by

(1.145)

This expression of the problem is the basis for fixed point solvers, where a starting K_eq,k is introduced.

1.7.4 Inter-harmonic coupling discussion

Harmonic balance Jacobian estimation based on local non-linear physics is rarely discussed in the literature. Possible for a given harmonic to estimate an equivalent local linear constitutive law. This notion comes close to quasi-Newton methods in transient simulations, that estimate Jacobians with fixed operators based on the non-linear constitutive laws.

The harmonic balance framework adds some complexity as the transient response is decomposed in the time domain. It seems however possible to assess Jacobian relevant topologies using Taylor series expansion when the non-linear forces can be expressed in a functional way (coherent with the existence of a constitutive law).

(1.146)

(1.147)

(1.148)

First order terms induced by the linear and cubic constraint term respectively induced by є=є + δ_cq1 cos(ω t) and є=є + δ_sq1 sin(ω t)

(1.149)

First order terms induced by the linear and cubic constraint term respectively induced by є=є + δ_cq3 cos(3ω t) and є=є + δ_sq3 sin(3 ω t)

(1.150)

First order terms induced by the linear and quadratic constraint term respectively induced by є=є + δ_cq1 cos(ω t) and є=є + δ_sq1 sin(ω t)

(1.151)

First order terms induced by the linear and quadratic constraint term respectively induced by є=є + δ_cq3 cos(2ω t) and є=є + δ_sq2 sin(2 ω t)

(1.152)

1.7.5 Harmonic load definition

The SDT definition of loads DofLoad or enforced displacement DofSet use the formalism of input shape matrices. Thus the time dependence of the load is given by

(1.153)

where [b]_{N× NS} describes the spatial content and the harmonic content is fully contained in {u(t)}. u is often scalar but can be a vector if multiple loads are combined. For each u component j, a two dimensional curve is defined giving

(1.154)

The harmonic load frequency dependence is then defined for each load vector b_j by scalar coefficients associated to each harmonic

(1.155)

These terms are declared by the field .curve defined in the Load structure.

The field can be omitted or left empty. It then assumed that the force is a constant value associated to the c1 harmonic.
The field can be a single curve, producing a scalar amplitude value per harmonic. The same amplitude is then applied to all loads.
The field can be a cell array of curves, each curve being associated to a column of Load.def. The result of each curve is then coherent with the harmonics declared in the curve.

A given curve entry will provide the frequency dependency of a given Load vector for a specified set of harmonic shapes.

It can be defined using a Tabular form, or a Functional form by a structure coherent with sdtweb curve formats, with fields

.X A 1x2 cell-array.
- .X{1}The first cell array provides the base frequency vector that will provide the linear interpolation coefficients. This can be left empty for functional definitions.
- .X{2} The second cell is a column cell-array providing the harmonic shape labels to whom the load is applied.
.Y the field providing the amplitude
- Tabular form: a matrix with as many lines a in field .X{1} and as many columns a the number of harmonics provided in field .X{2}.
- Functional form using MATLAB anonymous function handles. A structure with fields
  - .anonymous Provides the inline function. The anonymous function header is set by default if omitted, @(Zf,w). The inline can access the curve structure Zf whose fields will contain the fields declared in curve.Param, and the current frequency w.
  - .csv A parameter declaration string under cingui ParamEdit format. This declares the parameters to be used, with a default value, their type and a possible brief explaination.
  - .Param The current parameters. .Param can be
    - a string defining the parameters declared in the .csv by par1=val1 par2=val2 ...
    - a structure with fields corresponding exactly to the declared parameters struct('par1', val1, 'par2', val2).
    Any omitted parameter will be set to its default declared in the csv. Lack of default values would then results in an error at the function execution.
  - .tex a string providing a tex format of the formula used in .anonymous. Lazy declaration can be done by providing a string using w instead of the structure format.

% tabular formulation
C1=struct('X',{{1, ... % Frequencies can be a column vector if varying
    {'c0';'c1'}}}, ...  % harmonic labels, see hdof
    'Xlab',{{'Frequency','harm'}}, ... % Frequency Hz or rad/s
    'Y',[.1 50]);  % 
% Functional formulation
C2=struct('X',{{[],{'s1'}}},...
'Y',{{struct('anonymous','Zf.k1*w^2',...
'Param','k1=1e-2',....
'csv','k1(1#%g#"")',...
'tex','k_1 \omega^2')}});
% Lazy anonymous
C3=struct('X',{{[],{'c0';'c2'}}},'Y',{{'10*sqrt(w)','w^2'}});

Internally this definition is transformed to use a command matrix [b]_{Nhdof × Nhload}, with N_hdof the number of harmonic DOF defined in field hdof and N_hload the number of harmonic loading. One physical load is replicated by the number of harmonics specified in the curve field .X input to allow distinct amplitudes per harmonic.

At a given pulsation, a vector u(ω)_{Nhload × 1} is generated by parsing the curve inputs,

(1.156)

so that the total external harmonic load vector is expressed as

(1.157)