Time integration: Generalized-Alpha
Kynema time-integration is achieved with a Lie-group generalized-\(\alpha\) algorithm, described in detail by [@Bruls-etal:2012]. This algorithm is second-order time accurate and summarized here for completeness. For a system represented as Eq. (1) and with linearized form shown in Eq. (3), the generalized-\(\alpha\) algorithm can be written as in Algorithm 1, where an \(n\) superscript denotes evaluation at the \(n\)th time step, and \(\rho_\infty \in[0,1]\) is the numerical damping with \(\rho_\infty = 1\) being no damping and \(\rho_\infty=0\) being maximum damping. The term \(\underline{a}^{n} \in \mathbb{R}^k\) is an auxiliary variable [@Arnold-Bruls:2007] for computation and is initialized as \(\underline{a}^{0}= \dot{\underline{v}}^0\). The algorithm includes the left and right preconditioning matrices described in [@Bottasso-etal:2008],
respectively, where \(\underline{\underline{D}}_L,\underline{\underline{D}}_R \in \mathbb{R}^{k\times k}\), \(\underline{\underline{I}}_p\) is the \(p \times p\) identity matrix.
Important to the time-integration scheme is the exponential map, which can be written for \(k_n\) nodes
where
\(\underline{\underline{I}} \in \mathbb{R}^{3 \times 3}\) is the identity matrix, \(\underline{\psi}\) is the Cartesian rotation vector, \(\phi = || \underline{\psi} ||\), and the tilde operator is defined for a vector \(\underline{a}\in \mathbb{R}^3\) with entries \(a_i\) as
Algorithm 1: Lie group generalized-α time integration for one time step [@Bruls-etal:2012]
Require: \(\underline{q}^n\), \(\underline{v}^n\), \(\dot{\underline{v}}^n\), \(\underline{a}^n\), \(t^n\), \(\Delta t\), \(\rho_\infty\), atol, rtol, \(i_\mathrm{max}\), \(k\), \(m\)
\(\alpha_m = \frac{2 \rho_\infty - 1}{\rho_\infty+1}\)
\(\alpha_f = \frac{\rho_\infty}{\rho_\infty+1}\)
\(\gamma = 0.5 + \alpha_f - \alpha_m\)
\(\beta = 0.25 \left( \gamma + 0.5\right)^2\)
\(\beta^\prime = \frac{1-\alpha_m}{\Delta t^2 \beta (1-\alpha_f)}\)
\(\gamma^\prime = \frac{\gamma}{\Delta t \beta}\)
\(t^{n+1} := t^n + \Delta t\)
\(\dot{\underline{v}}^{n+1} := \underline{0}\)
\(\underline{\lambda}^{n+1} := \underline{0}\)
\(\underline{a}^{n+1} := (\alpha_f \dot{\underline{v}}^{n} - \alpha_m \underline{a}^n)/(1-\alpha_m)\)
\(\underline{v}^{n+1} := \underline{v}^n + \Delta t (1-\gamma) \underline{a}^n + \gamma \Delta t \underline{a}^{n+1}\)
\(\Delta \underline{q}^n := \underline{v}^n+(0.5-\beta) \Delta t \underline{a}^n + \beta \Delta t \underline{a}^{n+1}\)
\(\mathrm{err} := 2.0\) Initialize with any \(\mathrm{err} > 1.0\)
\(i := 0\)
While \(i \le i_\mathrm{max}\) and \(\mathrm{err} > 1.0\) do
\(i := i+1\)
\(\underline{q}^{n+1} := \exp_{\underline{\underline{R}}^3\times \mathrm{SO(3)}} ( \Delta t \Delta \underline{q}^n ) \circ \underline{q}^n\)
Calculate \(\underline{r}^{n+1}\); See Eq. (2)
Calculate \(\underline{\underline{S}}_t^{n+1}\); See Eq. (4)
Solve \(\underline{\underline{D}}_L(\beta \Delta t^2) \underline{\underline{S}}_t^{n+1} \underline{\underline{D}}_R(\beta \Delta t^2) \begin{bmatrix} \Delta \underline{x} \Delta \underline{\lambda} \end{bmatrix}= -\underline{\underline{D}}_L(\beta \Delta t^2) \underline{r}^{n+1}\)
\(\begin{bmatrix} \Delta \underline{x} \Delta \underline{\lambda} \end{bmatrix} := \underline{\underline{D}}_R(\beta \Delta t^2) \begin{bmatrix} \Delta \underline{x} \Delta \underline{\lambda} \end{bmatrix}\)
\(\Delta \underline{q}^n := \Delta \underline{q}^n + \Delta \underline{x}/\Delta t^n\)
\(\underline{v}^{n+1} := \underline{v}^{n+1} + \gamma^\prime \Delta \underline{x}\)
\(\dot{\underline{v}}^{n+1} := \dot{\underline{v}}^{n+1} + \beta^\prime \Delta \underline{x}\)
\(\underline{\lambda}^{n+1} := \underline{\lambda}^{n+1} + \Delta \underline{\lambda}\)
\(\mathrm{err} := \sqrt{ \frac{1}{k + m} \left( \sum_{i=1}^{k} \left( \frac{ \Delta x_i }{ \mathrm{atol} + \mathrm{rtol} \left| \Delta t \Delta q_i^n \right| } \right)^2 + \sum_{i=1}^{m} \left( \frac{ \Delta \lambda_i }{ \mathrm{atol} + \mathrm{rtol} \left| \lambda_i^{n + 1} \right| } \right)^2 \right) }\) See [@Arnold-Hante:2017]
end while
\(\underline{a}^{n+1} := \underline{a}^{n+1} + \dot{\underline{v}}^{n+1}\left( 1 - \alpha_f\right) / \left( 1 - \alpha_m\right)\)
Return: \(t^{n+1}\), \(\underline{q}^{n+1}\), \(\underline{v}^{n+1}\), \(\dot{\underline{v}}^{n+1}\), \(\underline{\lambda}^{n+1}\), \(\underline{a}^{n+1}\)
The so-called iteration matrix, \(\underline{\underline{S}}_t \in \mathbb{R}^{(k+m)\times (k+m)}\), for Eqs. (1) and (3) can be written
where the tangent matrix, \(\underline{\underline{T}}_{\mathbb{R}^3\times \mathrm{SO(3)}}(\underline{\psi}) \in \mathbb{R}^{k\times k}\), is
and the variation of the virtual rotation is related to the variation of the Cartesian rotation vector as
with [@Geradin-Cardona:2001]
Arnold, M., and O. Brüls. 2007. “Convergence of the Generalized-\(\alpha\) Scheme for Constrained Mechanical Systems.” Multibody System Dynamics 18: 185–202.
Arnold, M., and S. Hante. 2017. “Implementation Details of a Generalized-\(\alpha\) Differential-Algebraic Equation Lie Group Method.” Journal of Computational and Nonlinear Dynamics 2: 021002.
Bottasso, C. L., D. Dopico, and L. Trainelli. 2008. “On the Optimal Scaling of Index Three DAEs in Multibody Dynamics.” Multibody System Dynamics 19: 3–20.
Brüls, O., A. Cardona, and M. Arnold. 2012. “Lie Group Generalized-\(\alpha\) Time Integration Fo Constrained Flexible Multibody Systems.” Mechanism and Machine Theory, 121–37.
Géradin, M., and A. Cardona. 2001. Flexible Multibody Dynamics: A Finite Element Approach. Chichester: John Wiley & Sons.