Implicit trapezoidal rule-based, Newmark and Generalized-alpha solver

This solver represents a class of solvers, which are – in the undamped case – based on the implicit trapezoidal rule (in the view of Runge-Kutta methods). The interpolation of the quantities for one step includes the start and the end value of the time step, thus being called trapezoidal integration rule. In some special cases in Newmark’s method , the interpolation might only depend on the start value or the end value.

For now, all implemented solvers can be viewed as a generalization of Newmark’s method, but there are called differently in the solver interfaces

• Implicit trapezoidal rule (Newmark with \(\beta = \frac 1 4\) and \(\gamma = \frac 1 2\))

• Newmark’s method

• Generalized-\(\alpha\) method (\(=\) generalized Newmark method with additional parameters), see Chung and Hulbert for the original method and Arnold and Brüls for the application to multibody system dynamics.

Newmark and Generalized-alpha method

Newmark’s method has two parameters \(\beta\) and \(\gamma\). The main ideas are given in the following. First, displacements and velocities are linearly interpolated using the accelerations \(\ddot {\mathbf{q}}\) of the beginning of the time step (subindex ‘0’) and the end of the time step (subindex ‘T’). The SON displacements and velocities and for FON coordinates are given by (definition of \(\aalg\) will become clear later):

(37)\[\begin{split}{\mathbf{q}}_T & = & {\mathbf{q}}_0 + h \dot {\mathbf{q}}_0 + h^2 (\frac 1 2 -\beta) \aalg_0 + h^2 \beta \aalg_T \nonumber\ \dot {\mathbf{q}}_T & = & \dot {\mathbf{q}}_0 + h (1-\gamma) \aalg_0 + h\gamma \aalg_T \nonumber\\ {\mathbf{y}}_T & = & {\mathbf{y}}_0 + h (1-\gamma_\FO) \vel^0_\FO + h\gamma_\FO \vel^T_\FO\end{split}\]

Hereafter, the system equations are solved at the end of the time step (\(T\)) for the unknown accelerations as well as for FON and AEN coordinates.

Remarks:

• The system of equations may be solved for accelerations \(\ddot {\mathbf{q}}\), but also for displacements \({\mathbf{q}}\) or even velocities as unknowns while the remaining quantities are reconstructed from Eq. (37). In case of displacements as unknowns, a scaling of the Jacobian is necessary, see later.

• For consistency reasons, one may set \(\gamma_\FO = \gamma\), but currently we use \(\gamma_T = \frac 1 2\), leading to no numerical damping for ODE1 variables \({\mathbf{y}}\).

• In the extension to the so-called generalized-\(\alpha\) method , algorithmic accelerations \(\aalg\) are used in Eq. (37).

• Algorithmic accelerations are no longer equivalent to the time derivatives of displacements, \(\aalg \neq \ddot {\mathbf{q}}\); thus, both sets of variables are used independently. In case of Newmark or the implicit trapezoidal rule just use \(\aalg = \ddot {\mathbf{q}}\).

• Implicit solvers are also available with Lie groups, if according rigid body nodes (NodeRigidBodyRotVecLG) are used, for theory see Holzinger et al. .

For generalized-\(\alpha\), the algorithmic accelerations \(\aalg\) are computed from the recurrence relation

\[(1-\alpha_m){\mathbf{a}}_T + \alpha_m {\mathbf{a}}_0 = (1-\alpha_f) \ddot {\mathbf{u}}_T + \alpha_f \ddot {\mathbf{u}}_0\]

which can be resolved for the unknown \({\mathbf{a}}_T\),

\[{\mathbf{a}}_T = \frac{(1-\alpha_f) \ddot {\mathbf{u}}_T + \alpha_f \ddot {\mathbf{u}}_0 - \alpha_m {\mathbf{a}}_0}{(1-\alpha_m)}\]

For the first step, one can simply use \(\aalg_0 = \ddot {\mathbf{q}}_0\).

Parameter selection for Generalized-alpha

Compared to alternative implicit integration methods (including the Newmark method), the generalized-\(\alpha\) integrator’s parameters break down to one single parameter \(\rho_\infty\), which allows to chose numerical damping in a practical way.

Based on a simple single DOF mass-spring-damper model , having the eigen frequency \(\omega = 2\pi f\) with frequency \(f\) and period \(T=1/f\), the spectral radius \(\rho\) for the integrator defines the amount of damping for a given step size \(h\) related to \(T\), thus using the dimensionless step size \(\bar h=h/T\).

In Fig. 35 the spectral radius is shown versus \(\bar h\) for various spectral radii at infinity \(\rho_\infty\). Here, \(\rho_\infty\) specifies the numerical damping of very time step for large step sizes (or very high frequencies). An amount of \(\rho_\infty=0.9\) means that high frequency parts of the system ((\(\bar h \gg 1\)); high compared to the step rate) are damped to \(90\%\) in every step, reducing an initial value \(1\) to \(2.66e-5\) after 100 steps, which is already much larger than usual physical damping in many cases.

Furthermore, low frequency parts of the system (\(\bar h \ll 1\)) receive almost no numerical damping, see again Fig. 35. Exemplarily, consider \(\rho\) a low frequency situation with different \(\rho_\infty\):

• \(\rho(\bar h=0.01, \rho_\infty=0.9) = 1 - 1.13\cdot 10^{-9}\)

• \(\rho(\bar h=0.01, \rho_\infty=0.6) = 1 - 1.22\cdot 10^{-7}\)

which shows that numerical damping is very low for moderately small step sizes (100 steps for one oscillation).

Obviously, \(\rho_\infty\) does not have a large influence for very high or low frequencies in the system as long as it is \(\neq 1\) and we could even use \(\rho_\infty=0\). Regarding differential algebraic equations (DAEs), \(\rho_\infty<1\) allows to integrate index 3 DAEs. Typically a value of \(\rho_\infty=0.7\) leads to a stable integration, but values depend on the structure of the multibody system.

Once having chosen \(\rho_\infty\), all other parameters follow automatically , regarding the \(\alpha\)s

\[\alpha_m = \frac{2 \rho_\infty - 1}{\rho_\infty + 1}, \quad \alpha_f = \frac{\rho_\infty}{\rho_\infty + 1}\]

and Newmarks’s parameters,

\[\gamma = \frac{1}{2} - \alpha_m + \alpha_f, \quad \beta = \frac{1}{4}(1- \alpha_m + \alpha_f)^2\]

Fig. 35 Spectral radius for generalized-\(\alpha\) method depending on dimensionless step size \(\bar h=h/T\), in which \(T\) is the period of an equivalent single DOF mass-spring-damper system.

Newton iteration

Thus, the residuals at the end of the time step (\(T\)) read (put all terms to LHS):

(38)\[\begin{split}{\mathbf{r}}^\GA_\SO &=& {\mathbf{M}} \ddot {\mathbf{q}}_T + \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}^\mathrm{T}} \tlambda_T - {\mathbf{f}}_\SO({\mathbf{q}}_T, \dot {\mathbf{q}}_T, t) = 0\\ {\mathbf{r}}^\GA_\FO &=& \dot {\mathbf{y}}_T + \frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}^\mathrm{T}} \tlambda_T - {\mathbf{f}}_\FO({\mathbf{y}}_T, t) = 0\\ {\mathbf{r}}^\GA_\AE &=& {\mathbf{g}}({\mathbf{q}}_T, \dot {\mathbf{q}}_T, {\mathbf{y}}_T, \tlambda_T, t) = 0\end{split}\]

We consider two options for SON: (A) solve for unknown accelerations \(\acc_T\), or (B) for unknown displacements \({\mathbf{q}}_T\).

(A) Solve for unknown accelerations

The unknowns for the Newton method then are

(39)\[\txi^\GA_{k+1} = \vr{\acc_T}{{\mathbf{y}}_T}{\tlambda_T}\]

and at the beginning of the step, we have

(40)\[\txi^\GA_{k} = \vr{\acc_0}{{\mathbf{y}}_0}{\tlambda_0}\]

For the Newton method, we need to compute an update for the unknowns of Eq. (39), using the previous residual \({\mathbf{r}}_{k}\) and the inverse of the Jacobian \({\mathbf{J}}_{k}\) of Newton iteration \(k\),

\[\txi^\GA_{k+1} = \txi^\GA_{k} - {\mathbf{J}}^{-1} \left( \txi^\GA_{k} \right) \cdot {\mathbf{r}}^\GA \left( \txi^\GA_{k} \right)\]

The Jacobian has the following \(3 \times 3\) structure,

\[{\mathbf{J}} = \mr{{\mathbf{J}}_{\SO\SO}}{{\mathbf{J}}_{\SO\FO}}{{\mathbf{J}}_{\SO\AE}} {{\mathbf{J}}_{\FO\SO}}{{\mathbf{J}}_{\FO\FO}}{{\mathbf{J}}_{\FO\AE}} {{\mathbf{J}}_{\AE\SO}}{{\mathbf{J}}_{\AE\FO}}{{\mathbf{J}}_{\AE\AE}} = \mr{{\mathbf{J}}_{\SO\SO}}{\Null}{{\mathbf{J}}_{\SO\AE}} {\Null}{{\mathbf{J}}_{\FO\FO}}{{\mathbf{J}}_{\FO\AE}} {{\mathbf{J}}_{\AE\SO}}{{\mathbf{J}}_{\AE\FO}}{{\mathbf{J}}_{\AE\AE}}\]

in which we consider \({\mathbf{J}}_{\FO\SO}\) and \({\mathbf{J}}_{\SO\FO}\) to vanish in the current implementations, which means that coupling of ODE1 and ODE2 coordinates is only possible due to algebraic equations.

The remaining terms in the Jacobian are currently (or by default settings) evaluated as:

\[\begin{split}{\mathbf{J}}_{\SO\SO}&=&\frac{\partial {\mathbf{r}}^\GA_\SO}{\partial \acc} = \frac{\partial {\mathbf{r}}^\GA_\SO}{\partial {\mathbf{q}}} \frac{\partial {\mathbf{q}}}{\partial \acc} + \frac{\partial {\mathbf{r}}^\GA_\SO}{\partial \dot {\mathbf{q}}} \frac{\partial \dot {\mathbf{q}}}{\partial \acc} = h^2 \beta {\mathbf{K}} + h \gamma {\mathbf{D}} \nonumber \\ {\mathbf{J}}_{\SO\AE}&=&\frac{\partial {\mathbf{r}}^\GA_\SO}{\partial \tlambda} = \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \quad (\mbox{or } \frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}} \mbox{ for constraints at velocity level)} \nonumber \\ {\mathbf{J}}_{\FO\FO}&=&\frac{\partial {\mathbf{r}}^\GA_\FO}{\partial {\mathbf{y}}} \nonumber \\ {\mathbf{J}}_{\AE\SO}&=&\frac{\partial {\mathbf{r}}^\GA_\AE}{\partial \acc} = \frac{\partial {\mathbf{g}}}{\partial \acc} = \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \frac{\partial {\mathbf{q}}}{\partial \acc} + \frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}} \frac{\partial \dot {\mathbf{q}}}{\partial \acc} = h^2 \beta \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} + h \gamma \frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}} \nonumber \\ {\mathbf{J}}_{\AE\FO}&=&\frac{\partial {\mathbf{r}}^\GA_\AE}{\partial {\mathbf{y}}} \nonumber \\ {\mathbf{J}}_{\AE\AE}&=&\frac{\partial {\mathbf{r}}^\GA_\AE}{\partial \tlambda} = \frac{\partial {\mathbf{g}}}{\partial \tlambda}\end{split}\]

Note that some parts of the Jacobian are neglected, such as mass matrix and constraint Jacobian terms in \({\mathbf{J}}_{\SO\SO}\), which are usually of minor influence. Furthermore, Jacobians for state-dependent loads are neglected except for system-wide numerical Jacobians.

Once an update \({\mathbf{q}}^\mathrm{Newton}_{k+1}\) has been computed, the interpolation formulas (37) need to be evaluated before the next residual and Jacobian can be computed.

(B) Solve for unknown displacements

This approach is similar to the previous approach and follows exactly the algorithm given by Arnold and Brüls , however, extended for ODE1 variables, which are integrated by the (undamped) trapezoidal rule. Documentation will be added lateron.

Initial accelerations

For the solvers based on the implicit trapezoidal rule, initial accelerations are necessary in order to significantly increase the accuracy of the first time step. For this reason, the constraints \({\mathbf{g}}({\mathbf{q}}_0, \dot {\mathbf{q}}_0, {\mathbf{y}}_0, \tlambda_0, t) = 0\) in Eq. (33) are differentiated w.r.t. time,

(41)\[\dot {\mathbf{g}}({\mathbf{q}}_0, \dot {\mathbf{q}}_0, {\mathbf{y}}_0, \tlambda_0, t) = \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \dot {\mathbf{q}}_0 + \frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}}\ddot {\mathbf{q}}_0 + \frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}} \dot {\mathbf{y}}_0 + \frac{\partial {\mathbf{g}}}{\partial \tlambda} \dot \tlambda + \frac{\partial {\mathbf{g}}}{\partial t} = 0 .\]

Currently, we assume \(\frac{\partial {\mathbf{g}}}{\partial \tlambda} = 0\) for all further derivations on initial accelerations. For velocity level constraints, Eq. (41) is used to extract initial accelerations \(\ddot {\mathbf{q}}_0\),

\[\frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}}\ddot {\mathbf{q}}_0 = -\frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \dot {\mathbf{q}}_0 -\frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}} \dot {\mathbf{y}}_0 - \frac{\partial {\mathbf{g}}}{\partial t} .\]

Finally, the equations for the computation of the initial accelerations read for velocity level constraints, note that \({\mathbf{y}}_{init}\) are the nodal initial values for \({\mathbf{y}}\),

(42)\[\mr{{\mathbf{M}}}{\Null}{\frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}^\mathrm{T}}} {\Null}{{\mathbf{I}}}{\Null} {\frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}}}{\Null}{\Null} \vr{\ddot {\mathbf{q}}_0}{{\mathbf{y}}_0}{\tlambda_0} = \vr{{\mathbf{f}}_\SO({\mathbf{q}}_T, \dot {\mathbf{q}}_T, t)}{{\mathbf{y}}_{init}} {-\frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \dot {\mathbf{q}}_0-\frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}} \dot {\mathbf{y}}_0 - \frac{\partial {\mathbf{g}}}{\partial t}} ,\]

The term \(\frac{\partial {\mathbf{g}}}{\partial t}\) can only occur in case of user functions and therefore currently not implemented, and the ODE1 term \(\frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}} = 0\) is not used yet in constraints.

For position level constraints, we assume \(\frac{\partial {\mathbf{g}}}{\partial \dot {\mathbf{q}}} = 0\) and \(\frac{\partial {\mathbf{g}}}{\partial {\mathbf{y}}} = 0\) in Eq. (41) and perform a second derivation w.r.t. time,

(43)\[\ddot {\mathbf{g}}({\mathbf{q}}_0, \dot {\mathbf{q}}_0, {\mathbf{y}}_0, \tlambda_0, t) = \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}}^2} \dot {\mathbf{q}}_0^2 + 2 \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}} \partial t} \dot {\mathbf{q}}_0 + \frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \ddot {\mathbf{q}}_0 + \frac{\partial^2 {\mathbf{g}}}{\partial t^2} = 0 .\]

For position level constraints, Eq. (43) is used to extract initial accelerations \(\ddot {\mathbf{q}}_0\),

\[\frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} \ddot {\mathbf{q}}_0 = - 2 \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}} \partial t} \dot {\mathbf{q}}_0 - \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}}^2} \dot {\mathbf{q}}_0^2 - \frac{\partial^2 {\mathbf{g}}}{\partial t^2} .\]

Finally, the equations for the computation of the initial accelerations for position level constraints read

(44)\[\mr{{\mathbf{M}}}{\Null}{\frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}^\mathrm{T}}} {\Null}{{\mathbf{I}}}{\Null} {\frac{\partial {\mathbf{g}}}{\partial {\mathbf{q}}} }{\Null}{\Null} \vr{\ddot {\mathbf{q}}_0}{{\mathbf{y}}_0}{\tlambda_0} = \vr{{\mathbf{f}}_\SO({\mathbf{q}}_T, \dot {\mathbf{q}}_T, t)}{{\mathbf{y}}_{init}} {- 2 \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}} \partial t} \dot {\mathbf{q}}_0 - \frac{\partial^2 {\mathbf{g}}}{\partial {\mathbf{q}}^2} \dot {\mathbf{q}}_0^2 - \frac{\partial^2 {\mathbf{g}}}{\partial t^2}} ,\]

The linear system of equations, either Eq. (42) or Eq. (44), is solved prior to an implicit time integration if

simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True,

which is the default value.