Time integration

Solfec implements three variants of time integration of rigid body kinematics:

  • RIG_POS – a semi–explicit scheme with a positive energy drift, cf. NEW1 in [1]
  • RIG_NEG – a semi–explicit scheme with a negative energy drift, cf. NEW2 in [1]
  • RIG_IMP – a semi–implicit without energy drift (but rather, oscillation), cf. NEW3 in [1]

And two variants of time integration for pseudo–rigid and finite–element kinematics:

  • DEF_EXP – an explicit leap–frog scheme, cf. Section 5.1 in [2]
  • DEF_LIM – a linearly implicit leap–frog scheme, [3] and TR1

A time integration scheme is selected by modifying the scheme parameter of the BODY object.

Rigid integration

Linear motion is integrated like deformable motion. Rigid rotations are integrated as follows

(31)\[\mathbf{\Lambda}^{t+\frac{h}{2}}=\mathbf{\Lambda}^{t}\exp\left[\frac{h}{2}\mathbf{\Omega}^{t}\right]\]
(32)\[\mathbf{T}^{t+\frac{h}{2}}=\left(\mathbf{\Lambda}^{t+\frac{h}{2}}\right)^{T}\mathbf{t}^{t+\frac{h}{2}}\]
(33)\[\mathbf{\Omega}^{t+\frac{h}{2}}=\mathbf{J}^{-1}\left[\exp\left[-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+\frac{h}{2}\mathbf{T}^{t+\frac{h}{2}}\right]\]
(34)\[\mathbf{\Omega}_{1}^{t+h}=\mathbf{\Omega}^{t}+\mathbf{J}^{-1}h\left[\mathbf{T}^{t+\frac{h}{2}}-\mathbf{\Omega}^{t+\frac{h}{2}}\times\mathbf{J}\mathbf{\Omega}^{t+\frac{h}{2}}\right]\]

If explicit

(35)\[\mathbf{\Lambda}^{t+h}=\mathbf{\Lambda}^{t+\frac{h}{2}}\exp\left[\frac{h}{2}\mathbf{\Omega}_{1}^{t+h}\right]\]
(36)\[\mathbf{\Omega}_{2}^{t+h}=\mathbf{J}^{-1}\exp\left[-\frac{h}{2}\mathbf{\Omega}_{1}^{t+h}\right]\left[\exp\left [-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+h\mathbf{T}^{t+\frac{h}{2}}\right]\]

otherwise

(37)\[\textbf{solve} \left(\exp\left[\frac{h}{2}\mathbf{\Omega}_{3}^{t+h}\right]\mathbf{J}\mathbf{\Omega}_{3}^{t+h}= \exp\left[-\frac{h}{2}\mathbf{\Omega}^{t}\right]\mathbf{J}\mathbf{\Omega}^{t}+h\mathbf{T}^{t+\frac{h}{2}}\right)\]
(38)\[\mathbf{\Lambda}^{t+h}=\mathbf{\Lambda}^{t+\frac{h}{2}}\exp\left[\frac{h}{2}\mathbf{\Omega}_{3}^{t+h}\right]\]

The scheme ending at (35) is DEF_POS, ending at (36) is DEF_NEG, and using instead (37) and (38) is DEF_IMP. Above, \(\exp\left[\cdot\right]\) is the exponential map defined by the Rodrigues formula

\[\exp\left[\mathbf{\Psi}\right]=\mathbf{I}+\frac{\sin\left\Vert \mathbf{\Psi}\right\Vert }{\left\Vert \mathbf{\Psi}\right\Vert } \hat{\mathbf{\Psi}}+\frac{1-\cos\left\Vert \mathbf{\Psi}\right\Vert }{\left\Vert \mathbf{\Psi}\right\Vert ^{2}}\hat{\mathbf{\Psi}}^{2}\]

where \(\mathbf{I}\) is the \(3\times3\) identity operator, \(\hat{\mathbf{\Psi}}\) creates the skew symmetric matrix out of a 3-vector \(\mathbf{\Psi}\), and \(\left\Vert \cdot\right\Vert\) stands for the Euclidean norm. The time step is denoted as \(h\).

Deformable integration

Deformable time integrator reads

\[\mathbf{q}^{t+\frac{h}{2}}=\mathbf{q}^{t}+\frac{h}{2}\mathbf{u}^{t}\]
\[\mathbf{u}^{t+h}=\mathbf{u}^{t}+\mathbf{A}^{-1}h\mathbf{f}\left(\mathbf{q}^{t+\frac{h}{2}},\mathbf{u}^{t}\right)\]
\[\mathbf{q}^{t+h}=\mathbf{q}^{t+\frac{h}{2}}+\frac{h}{2}\mathbf{u}^{t+h}\]

where in the explicit case

(39)\[\mathbf{A}=\mathbf{M}\text{ for DEF_EXP}\]

and in the linearly implicit case

(40)\[\mathbf{A}=\mathbf{M}+\left(\frac{\eta h}{2}+\frac{h^{2}}{4}\right)\mathbf{K}\left(\mathbf{q}^{t+h/2}\right)\text{ for DEF_LIM}\]

The time step is denoted as \(h\). See TR1 for technical details.

Implementation

Time integration is implement in bod.c (rigid, pseudo–rigid) and fem.c (finite–element) files. Inverse generalized inertia matrix \(\mathbf{A}^{-1}\) is declared in bod.h as follows:

struct general_body
{
  /* ... */

  MX *inverse;      /* generalized inverse inertia oprator */

  /* ... */
}

Rigid integration formulae (31)-(34) are in bod.c:BODY_Dynamic_Step_Begin.
Rigid integration formulae (35)-(38) are in bod.c:BODY_Dynamic_Step_End.
Pseudo–rigid integration is included in the same routines: first half–step and second half–step.
Finite–element, total Lagrangian formulation based, integration formulae (5) and (6) are in fem.c:TL_dynamic_step_begin.
Finite–element, total Lagrangian formulation based, integration formula (7) is in fem.c:TL_dynamic_step_end.

[1](1, 2, 3) IJNME, 81(9):1073–1092, 2010.
[2]Koziara, PhD thesis, 2008.
[3]ANM, 25(2–3): 297–302, 1997.