The general dynamic equation of motion in the time domain is:

(1–1)

where:

,, and are the acceleration, velocity, and displacement vectors.

, , are the mass, damping, and stiffness matrices.

is the vector of applied (external) excitation.

is the vector of internal nonlinear forces.

Periodic excitation and response may be assumed and enforced by the following equations [1].

(1–2)

(1–3)

(1–4)

where:

represents the real part of a complex quantity.

is the time domain harmonic index.

is the number of harmonics (specified as NH on HROPT,HBM,NH).

is the circular frequency of the period.

is time.

(¯) represents the complex Fourier coefficient of a quantity.

Substituting Equation 1–2 – Equation 1–4 into Equation 1–1 and equating coefficients of gives the balance equations describing the dynamics of a system in the frequency domain for time harmonic :

(1–5)

(1–6)

Collating the time harmonics in Equation 1–6 formulates the set of algebraic equations:

(1–7)

The expanded version of the multiharmonic Equation 1–7 expressed with complex quantities is:

(1–8)

The corresponding expanded version of Equation 1–7 expressed with real quantities is:

(1–9)

where:

(1–10)

(1–11)

(1–12)

and represents the imaginary part of a complex quantity.

Newton-Raphson Method and Numerical Continuation

Equation 1–9 can be solved to obtain the forced response at a given frequency . A Newton-Raphson procedure is used to iterate from an initial guess (specified via HBMOPT,UINIT) [5, 6]. The nonlinear equation solver uses a trust-region algorithm based on a gradient-descent method [5,8].

The residual of the balance equation (Equation 1–9) is given by:

(1–13)

The solver tries to satisfy the condition iteratively making guesses of using a Jacobian described as [1, 2]:

(1–14)

A numerical continuation algorithm may be optionally used with the Newton-Raphson type solver. In this case, the frequency is treated as an additional solution variable, and the solver tries to satisfy the following condition:

where:

is the augmented multiharmonic solution vector.

is the circular frequency (rad/s).

Numerical continuation enables automatic determination of the frequency step during solution [7, 8].

You can improve the convergence rate by scaling the augmented multiharmonic solution vector using the HBMOPT,SCAL command. The following scaling strategies are available:

You can also scale the circular frequency by a scalar value (FSCAL) using the HBMOPT,SCAL command:

Scaling of the multiharmonic solution vector only, of the frequency only, and of both (augmented multiharmonic vector) are supported.

At any given iteration of the nonlinear equation solver, the solution is considered to be converged if the norm of the residual is lower than a tolerance criterion :

(1–15)

The solver also stops at the k^th iteration, irrespective of convergence if the norm of the difference between consecutive solution estimates is below a tolerance criterion :

(1–16)

where:

is the solution estimate at the k^th iteration, and

for continuation.

The solver also stops at the k^th iteration, irrespective of convergence if the norm of the gradient estimate is below a tolerance criterion :

(1–17)

where:

is the gradient estimate at the k^th iteration, .

These conditions on the norm of the difference between consecutive solution estimates (step size) , residual norm , and gradient norm are stopping criteria. Specify their values as UTOL, RTOL, and GTOL on HBMOPT,TRTOL.

The numerical continuation algorithm stops if one of the following conditions occurs:

and in Equation 1–13 and Equation 1–14 are obtained using an alternating frequency-time (AFT) procedure. At any given iteration step, the solution estimates and parameters are used to calculate the equivalent periodic time domain displacements using backward fast Fourier transforms (backward FFTs). Note that this is equivalent to enforcing Equation 1–2. The initial velocity at the first time step of the period may also be calculated using similar equations.

A transient analysis is performed on a nonlinear submodel with the calculated quantities applied as boundary conditions (BCs) and initial conditions (ICs). The ordinary differential equation representing the dynamics of the nonlinear submodel is given by:

(1–18)

where: , , are the mass, damping, and stiffness matrices of the nonlinear submodel obtained by assembling the nonlinear elements only.

The internal forces at nonlinear elements are assembled to calculate the internal force vector at each time step. This calculation is repeated for time points spanning multiple periods until periodic convergence is achieved or some user-specified limit is reached. Then, forward fast Fourier transforms (FFTs) are used to convert back the periodic time domain internal forces to individual harmonics in the frequency domain. These harmonic terms are the components of the force vector on the right-hand side of Equation 1–8. Similar calculations are performed to obtain the Jacobian .