Time integration

Solfec-1.0 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

(39)\[\mathbf{q}^{t+\frac{h}{2}}=\mathbf{q}^{t}+\frac{h}{2}\mathbf{u}^{t}\]

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

(41)\[\mathbf{q}^{t+h}=\mathbf{q}^{t+\frac{h}{2}}+\frac{h}{2}\mathbf{u}^{t+h}\]

where in the explicit case

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

and in the linearly implicit case

(43)\[\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 (39) and (40) are in fem.c:TL_dynamic_step_begin.
Finite–element, total Lagrangian formulation based, integration formula (41) 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.