Equations of motion

1. Kinetics for a rigid body moving in water

\[\begin{array}{ll}\boldsymbol{\dot P_B}&=\boldsymbol{F}\\\\\boldsymbol{\dot L_B}&=\boldsymbol{M}\end{array}\]

where:

\(\boldsymbol{P_B}\)

linear momentum

\(\boldsymbol{L_B}\)

angular momentum about some reference point

\(\boldsymbol{F}\)

external force

\(\boldsymbol{M}\)

external moment about the reference point

As reference point for \(\mathrm {\boldsymbol{L_B}}\) and \(\mathrm {\boldsymbol{M}}\) , the origin of the body-fixed system of axes will be chosen.

The linear and angular momenta are given by:

\[\begin{array}{ll}\boldsymbol{P_B}&=\boldsymbol{m}(\boldsymbol{v}+\boldsymbol{w}\times \boldsymbol{r_c})\\\\\boldsymbol{L_B}&=\boldsymbol{I\omega }+\boldsymbol{mr_c}\times \boldsymbol{v}\end{array}\]

where:

\(\boldsymbol{m}\)

body mass

\(\boldsymbol{v}\)

velocity of body origin

\(\boldsymbol{\omega}\)

angular velocity of body

\(\boldsymbol{r_c}\)

position of body centre of mass relative to body origin

\(\boldsymbol{I}\)

body’s inertia tensor with respect to body origin

Using ’ to denote time derivatives with respect to the body-fixed system of reference, Equation (1) can be written:

\[\begin{array}{ll}\boldsymbol{P'_B}+\boldsymbol{\omega }\times \boldsymbol{P_B}&=\boldsymbol{F}\\\\\boldsymbol{L'_B}+\boldsymbol{\omega }\times \boldsymbol{L_B}+\boldsymbol{v}\times \boldsymbol{P_B}&=\boldsymbol{M}\end{array}\]

Here, the second term in both equations is an effect of the rotating frame of reference, while the third term in the second equation is due to the translation of the body origin. When the centre of mass coincides with the origin, \(\mathrm {\boldsymbol{r_c}=[0,\,0,\,0]^T}\) , this term will vanish since the velocity and the linear momentum then will be collinear.

Let the vectors in Equation (2) be represented by their triples of coordinates with respect to the body coordinate system and let the inertia matrix be represented by a \(\mathrm {3\times 3}\) matrix with respect to the same axes.

Writing:

\[\boldsymbol{r_c}=(x_c,y_c,z_c)^T\]

the vector cross-multiplication occurring in Equation (2) is equivalent to multiplication from the left by the matrix

\[\boldsymbol{R}=\begin{bmatrix}0&-z_c&y_c\\z_c&0&-x_c\\-y_c&x_c&0\end{bmatrix}\]

Similarly, the operation is equivalent to left multiplication by \(\mathrm {-\boldsymbol{R}=\boldsymbol{R}^T}\) .

Using \(\mathrm {\boldsymbol{\Gamma}}\) to denote the \(\mathrm {3\times 3}\) identity matrix and defining the matrices

\[\begin{array}{lll} \boldsymbol{B_{11}} &= m\boldsymbol{\Gamma},\\ \boldsymbol{B_{12}} &= m\boldsymbol{R}^T \\\\ \boldsymbol{B_{21}} &= \boldsymbol{R},\\ \boldsymbol{B_{22}} &= \boldsymbol{I} \end{array}\]

Equation (2) can be written:

\[\begin{array}{ll}\boldsymbol{P_B}&=\boldsymbol{B_{11}v}+\boldsymbol{B_{12}\omega }\\\\\boldsymbol{L_B}&=\boldsymbol{B_{21}v}+\boldsymbol{B_{22}\omega }\end{array}\]

and Equation (3) becomes:

\[\begin{array}{ll}\boldsymbol{B_{11}v'}+\boldsymbol{B_{12}\omega '}+\boldsymbol{\omega }\times \boldsymbol{P_B}&=\boldsymbol{F}\\\\\boldsymbol{B_{21}v'}+\boldsymbol{B_{22}\omega }+\boldsymbol{\omega } \times \boldsymbol{L_B}+\boldsymbol{v}\times \boldsymbol{P_B}&=\boldsymbol{M}\end{array}\]

The 6 x 6 matrix =

\[\boldsymbol{B}=\begin{bmatrix}\boldsymbol{B_{11}}&\boldsymbol{B_{12}}\\\boldsymbol{{B_{21}}}&\boldsymbol{B_{22}}\\\end{bmatrix}\]

will be referred to as the body-mass matrix.

A detailed derivation of the result Equation (8) presented in a different form is given by (Abkowitz, 1964). The form of Equation (8) was chosen to obtain similarity with the expressions for hydrodynamic added-mass force which are presented in the following.

When the body is moving in a fluid, hydrodynamic forces are included via so-called added mass. These forces and moments are included in the external \(\mathrm {\boldsymbol{F}}\) and \(\mathrm {\boldsymbol{M}}\) in Equation (8), i.e.:

\[\begin{array}{ll}\boldsymbol{F}&=\boldsymbol{F_1}+\boldsymbol{F_A}\\\\\boldsymbol{M}&=\boldsymbol{M_1}+\boldsymbol{M_A}\end{array}\]

where subscript \(\mathrm {\boldsymbol{A}}\) refers added mass forces.

Assuming the fluid to be ideal, irrotational and unbounded in all directions, \(\mathrm {\boldsymbol{F_{A}}}\) and \(\mathrm {\boldsymbol{M_{A}}}\) can be associated with added linear and angular momenta which can be expressed as:

\[\begin{array}{ll} \boldsymbol{P_A}& = \boldsymbol{A_{11}} (\boldsymbol{v} - \boldsymbol{U}) + \boldsymbol{A_{12}} \boldsymbol{\omega} \\ \\ \boldsymbol{L_A}& = \boldsymbol{A_{21}} (\boldsymbol{v} - \boldsymbol{U}) + \boldsymbol{A_{22}} \boldsymbol{\omega} \end{array}\]

where \(\boldsymbol{A_{ij}}\) are added-mass tensors with respect to the body origin and \(\boldsymbol{U}\) is a homogeneous steady current velocity. The tensors will be represented by \(\mathrm {3\times 3}\) matrices in the body system of coordinates, and the \(\mathrm {6\times 6}\) matrix

\[\boldsymbol{A}=\begin{bmatrix} \boldsymbol{A_{11}} & \boldsymbol{A_{12}} \\ \boldsymbol{A_{21}} & \boldsymbol{A_{22}} \end{bmatrix}\]

will be referred to as the added-mass matrix.

The added mass forces and moments are given by

\[\begin{array}{ll} \boldsymbol{A_{11}v'} + \boldsymbol{A_{12}\omega '} + \boldsymbol{\tilde{\omega} } \times \boldsymbol{\tilde{P}_A} &= \boldsymbol{-F_A} \\ \\ \boldsymbol{A_{21}v'} + \boldsymbol{A_{22}\omega '} + \boldsymbol{\tilde{\omega} } \times \boldsymbol{\tilde{L}_A} + \boldsymbol{\tilde{v}} \times \boldsymbol{\tilde{P}_A} &= \boldsymbol{-M_A} \end{array}\]

which is completely analogous to Equation (8). When the fluid is unbounded in all directions, we can make the substitutions \(\boldsymbol{\tilde{v}} = \boldsymbol{v} - \boldsymbol{U}\), \(\boldsymbol{\tilde{\omega}} = \boldsymbol{\omega}\), \(\boldsymbol{\tilde{P}_A} = \boldsymbol{P_A}\) and \(\boldsymbol{\tilde{L}_A} = \boldsymbol{L_A}\). However, when a free surface is present and added mass results from linear theory, different pragmatic substitutions depending on the type of analysis is used by SIMO:

\[\begin{array}{ll} \boldsymbol{\tilde{v}} = - \boldsymbol{U} \\ \\ \boldsymbol{\tilde{\omega}} = \mathbf{0} \end{array}\]

The option name Block low pass translation velocities in Hydro Filter Method has a misleading name since both translational and rotational velocities are set to zero when the option is applied. The name will change in a future release.
\[\begin{array}{ll} \boldsymbol{\tilde{v}} = \boldsymbol{v} - \boldsymbol{U} \\ \\ \boldsymbol{\tilde{\omega}} = \begin{bmatrix} 0 & 0 & \omega_{3} \end{bmatrix} \end{array}\]
\[\begin{array}{ll} \boldsymbol{\tilde{v}} = \boldsymbol{v}_{LF} - \boldsymbol{U} \\ \\ \boldsymbol{\tilde{\omega}} = \begin{bmatrix} 0 & 0 & \omega_{3,LF} \end{bmatrix} \end{array}\]

  • If Hydro Filter Method is set to No filtering or Filter activated and Body Type 6 DOF - separated analysis is used, such that \(\boldsymbol{v}\) and \(\boldsymbol{\omega}\) can be interpreted as low-frequency velocities, we have that

\[\begin{array}{ll} \boldsymbol{\tilde{v}} = \boldsymbol{v} - \boldsymbol{U} \\ \\ \boldsymbol{\tilde{\omega}} = \boldsymbol{\omega} \end{array}\]

In all cases we write,

\[\begin{array}{ll} \boldsymbol{\tilde{P}_A}& = \boldsymbol{A_{11}} \boldsymbol{\tilde{v}} + \boldsymbol{A_{12}} \boldsymbol{\tilde{\omega}} \\ \\ \boldsymbol{\tilde{L}_A}& = \boldsymbol{A_{21}} \boldsymbol{\tilde{v}} + \boldsymbol{A_{22}} \boldsymbol{\tilde{\omega}} \end{array}\]

Further, when Body Type 6 DOF - Time domain is used, added mass is taken at the high-frequency limit, and frequency dependency due to free-surface effects is included in \(\mathrm {\boldsymbol{F_{1}}}\) and \(\mathrm {\boldsymbol{M_{1}}}\) via convolution integrals. When Body Type 6 DOF - separated analysis is used, added mass is taken at the low-frequency limit.

Equation (13) is seen to agree with the well-known paradox of D’Alembert, i.e. there is no force on a body of arbitrary shape which moves with constant velocity without rotating in an ideal fluid. However, as is also known, a moment will exist. This moment is given by the term \(\mathrm {\boldsymbol{\tilde{v}}\times \boldsymbol{\tilde{P}_A}}\) in Equation (13).

Adding Equation (8) and Equation (13) and defining the virtual-mass matrix

\[\boldsymbol{V}=\begin{bmatrix}\boldsymbol{V_{11}}&\boldsymbol{V_{12}}\\\boldsymbol{V_{21}}&\boldsymbol{V_{22}}\end{bmatrix}=\begin{bmatrix}\boldsymbol{A_{11}}+\boldsymbol{B_{11}}&\boldsymbol{A_{12}}+\boldsymbol{B_{12}}\\\boldsymbol{A_{21}}+\boldsymbol{B_{21}}&\boldsymbol{A_{22}}+\boldsymbol{B_{22}}\\\end{bmatrix}=\boldsymbol{A}+\boldsymbol{B}\]

the equation of motion becomes

\[\boldsymbol{V} \begin{bmatrix} \boldsymbol{v'} \\ \boldsymbol{\omega '} \end{bmatrix} + \begin{bmatrix} \boldsymbol{\omega } \times \boldsymbol{P_B} + \boldsymbol{\tilde{\omega} } \times \boldsymbol{\tilde{P}_A} \\ \boldsymbol{\omega } \times \boldsymbol{L_B} + \boldsymbol{v} \times \boldsymbol{P_B} + \boldsymbol{\tilde{\omega}} \times \boldsymbol{\tilde{L}_A} + \boldsymbol{\tilde{v}} \times \boldsymbol{\tilde{P}_A} \end{bmatrix} = \begin{bmatrix} \boldsymbol{F_1} \\ \boldsymbol{M_1} \end{bmatrix}\]

The accelerations occurring in Equation (20) give the rate of change of velocity relative to the rotating body system. It is desirable to refer the accelerations to a fixed frame of reference. The relationship is

\[\begin{array}{ll}\boldsymbol{v'}&=\boldsymbol{\dot v}-\boldsymbol{\omega }\times \boldsymbol{v}\\\\\boldsymbol{\omega '}&=\boldsymbol{\dot \omega }-\boldsymbol{\omega }\times \boldsymbol{\omega }=\boldsymbol{\dot \omega }\end{array}\]

Substituting Equation (21) in Equation (20) gives:

\[\begin{bmatrix}\boldsymbol{\dot v}\\\boldsymbol{\dot \omega }\end{bmatrix} = \boldsymbol{V^{-1}} \begin{bmatrix} -\boldsymbol{\omega } \times \boldsymbol{P_B} - \boldsymbol{\tilde{\omega} } \times \boldsymbol{\tilde{P}_A} + \boldsymbol{F_1} \\ -\boldsymbol{\omega } \times \boldsymbol{L_B} - \boldsymbol{v} \times \boldsymbol{P_B} - \boldsymbol{\tilde{\omega}} \times \boldsymbol{\tilde{L}_A} - \boldsymbol{\tilde{v}} \times \boldsymbol{\tilde{P}_A} + \boldsymbol{M_1} \end{bmatrix} + \begin{bmatrix} \boldsymbol{\omega } \times \boldsymbol{v} \\ \boldsymbol{0} \end{bmatrix}\]

When body coordinates are used to represent the vectors on the right-hand side of Equation (22), and the matrix \(\mathrm {\boldsymbol{V}}\) is defined accordingly, the coordinates of \(\mathrm {\boldsymbol{\dot v}}\) and \(\mathrm {\boldsymbol{\dot \omega }}\) also refer to the body system - or rather, the fixed system which at that particular moment coincides with the rotating body system. To express the acceleration vectors in coordinates of the fixed, global system, a coordinate transformation is necessary.

2. Coordinate transformation

The orientation of the body system of axes is thought of as being derived from an original identification with the global system by three rotations in succession.

  • First, the body is rotated by an angle \(\mathrm {\psi}\) about the common z-axis. This alters the direction of the body x-axis and y-axis.

  • Second, the body is rotated by an angle \(\mathrm {\theta }\) about its y-axis. This changes the orientation of the x-axis and the z-axis, but the body y-axis still lies in the global xy-plane.

  • Finally, the body is rotated by an angle \(\mathrm {\phi }\) about its x-axis, altering the orientation of the y-a-xis and z-axis.

The orientation of the body system is now uniquely determined by the angles, \(\mathrm {\phi }\) , \(\mathrm {\theta }\) , \(\mathrm {\psi}\) and the specified order of rotation.

The matrix for transformation of a vector’s coordinates in the body system to the representation in the global system becomes:

\[\boldsymbol{\Lambda}=\begin{bmatrix}\cos\psi\cos\theta &-\sin\psi\cos\phi +\cos\psi\sin\theta \sin\phi &\sin\psi\sin\phi +\cos\psi\sin\theta \cos\phi \\\sin\psi\cos\theta &\cos\psi\cos\phi +\sin\psi\sin\theta \sin\phi &-\cos\psi\sin\phi +\sin\psi\sin\theta \cos\phi \\-\sin\theta &\cos\theta \sin\phi &\cos\theta \cos\phi \end{bmatrix}\]

i.e. if a vector is represented by \(\mathrm {\boldsymbol{x^{(B)}}}\) in the body system, the representation in the global system is

\[\boldsymbol{x^{(G)}}=\boldsymbol{\Lambda x^{(B)}}\]

Conversely,

\[\boldsymbol{x^{(B)}}=\boldsymbol{\Lambda^Tx^{(G)}}\]

since \(\mathrm {\boldsymbol{\Lambda^{-1}}=\boldsymbol{\Lambda^T}}\)

3. Position vector. Model formulation

The position vector of the body is defined as

\[\boldsymbol{x}=(x,y,z,\phi ,\theta ,\psi)^T\]

where \(\mathrm {x,y,z}\) are coordinates of the body origin with respect to the global coordinate system (Strictly, x is not a vector since \(\mathrm {\phi ,\theta }\) and \(\mathrm {\psi}\) are not coordinates)

It is desirable to express the dynamic model in the form:

\[\boldsymbol{\ddot x}=\boldsymbol{f(x,\dot x,\xi )}\]

where \(\mathrm {\boldsymbol{\xi }}\) is a vector of inputs not depending on \(\mathrm {\boldsymbol{x}}\) or \(\mathrm {\boldsymbol{\dot x}}\) , e.g. speed and direction of wind and current, undisturbed surface elevation. When the body is coupled to other bodies, \(\mathrm {x}\) will include the positions or velocities of these bodies as well.

The formulation Equation (27) is preferred to that of the state-space formulation, since a wider class of numerical methods for solution is available.

Consider the vectors:

\[\begin{array}{ll} \boldsymbol{\dot x} &= (\dot x,\dot y,\dot \phi ,\dot \theta ,\dot \psi)^T \\\\ \boldsymbol{\ddot x} &= (\ddot x,\ddot y,\ddot \phi ,\ddot \theta ,\ddot \psi)^T \end{array}\]

The sub-vectors \(\mathrm {(\dot x,\dot y,\dot z)^T}\) and \(\mathrm {(\ddot x,\ddot y,\ddot z)^T}\) are the vectors \(\mathrm {\boldsymbol{v}}\) and \(\mathrm {\boldsymbol{\dot v}}\) of Equation (22) in global-system representation, i.e.

\[\begin{bmatrix}\dot x\\\dot y\\\dot z\\\end{bmatrix}=\boldsymbol{\Lambda v},\quad \begin{bmatrix}\ddot x\\\ddot y\\\ddot z\\\end{bmatrix}=\boldsymbol{\Lambda\dot v}\]

Let the vector \(\mathrm {\boldsymbol{\omega }}\) in Equation (22) be expressed as

\[\boldsymbol{\omega }=(p,q,r)^T\]

in the body system. The following relationship can be shown to exist, cf. (Norrbin, 1971):

\[\begin{array}{ll}\dot \phi &=p+q\sin\phi \tan\theta +r\cos\phi \tan\theta \\\\\dot \theta &=q\cos\phi -r\sin\phi \\\\\dot \psi&=q\sin\phi \sec\theta +r\cos\phi \sec\theta \end{array}\]
\[\boldsymbol{\gamma }=(\phi ,\theta ,\psi)^T\]

Equation (32) can be expressed in vector form:

\[\boldsymbol{\dot \gamma }=\boldsymbol{f}(\boldsymbol{\gamma },\boldsymbol{\omega })=\boldsymbol{M}(\boldsymbol{\gamma })\boldsymbol{\omega }\]

where:

\[\boldsymbol{M}=\begin{bmatrix}1&\sin\phi \tan\theta &\cos\phi \tan\theta \\0&\cos\phi &-\sin\phi \\0&\sin\phi \sec\theta &\cos\phi \sec\theta \end{bmatrix}\]

Differentiating Equation (33) with respect to time gives

\[\boldsymbol{\ddot \gamma }=\displaystyle \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{\gamma }}\boldsymbol{\dot \gamma }+\boldsymbol{M\dot \omega }=\displaystyle \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{\gamma }}\boldsymbol{M\omega }+\boldsymbol{M\dot \omega }=\begin{bmatrix}\ddot \phi \\\ddot \theta \\\ddot \psi\end{bmatrix}\]

where

\[\displaystyle \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{\gamma }}=\begin{bmatrix}(q\cos\phi -r\sin\phi )\tan\theta &(q\sin\phi +r\cos\phi )\sec^2\theta &0\\-q\sin\phi -r\cos\phi &0&0\\(q\cos\phi -r\sin\phi )\sec\theta &(q\sin\phi +r\cos\phi )\tan\theta \sec\theta &0\\\end{bmatrix}\]

Assuming, \(\mathrm {\boldsymbol{x}}\) , \(\mathrm {\boldsymbol{\xi }}\)
and \(\mathrm {\boldsymbol{\dot x}}\) are known, Equation (27) is realized as follows:

  1. The velocity \(\mathrm {(\dot x,\dot y,\dot z)^T}\) is transformed to body-coordinate representation \(\mathrm {\boldsymbol{v}}\) using the transpose of the matrix \(\mathrm {\boldsymbol{\Lambda}}\) in Equation (22).

  2. \(\mathrm {\boldsymbol{\omega }}\) is calculated from \(\mathrm {\dot \gamma }\) using the inverse of Equation (3).

  3. The external force \(\mathrm {\boldsymbol{F_{1}}}\) and moment \(\mathrm {\boldsymbol{M_{1}}}\) in Equation (21) are calculated using \(\mathrm {\boldsymbol{x}}\), \(\mathrm {\boldsymbol{v}}\), \(\mathrm {\boldsymbol{\omega }}\) and \(\mathrm {\xi }\) .

  4. The accelerations \(\mathrm {\boldsymbol{\dot v}}\) and \(\mathrm {\boldsymbol{\dot \omega }}\) are obtained from Equation (21).

  5. \(\mathrm {\boldsymbol{\dot v}}\) is transformed to a representation in global system using the matrix \(\mathrm {\boldsymbol{\Lambda}}\)

  6. \(\mathrm {\boldsymbol{\ddot \gamma }}\) is calculated using Equation (35).

Thus \(\mathrm {\boldsymbol{\ddot x}=(\boldsymbol{\dot v^T},\boldsymbol{\ddot \gamma ^T})^T}\)
is established.

4. Methods for numerical integration

Three methods are available. In the following description \(\mathrm {T}\) is the time step. The following shortened notation is used:

\[\begin{array}{ll}\boldsymbol{x_k}&=\boldsymbol{x}(t_k)\\\\\boldsymbol{\dot x_k}&=\boldsymbol{\dot x}(t_k)\\\\\boldsymbol{\xi _k}&=\boldsymbol{\xi }(t_k)\\\\\boldsymbol{f_k}&=\boldsymbol{f}(\boldsymbol{x_k,\dot x_k,\xi _k})\end{array}\]

4.1. Modified Euler method

\[\begin{array}{lll}\boldsymbol{\dot x_{k+1}}&=\boldsymbol{\dot x_k}+T\boldsymbol{f_k}&\text{(a)}\\\\\boldsymbol{x_{k+1}}&=\boldsymbol{x_k}+T\boldsymbol{\dot x_{k+1}}&\text{(b)}\end{array}\]

The method differs from the ordinary method in that \(\mathrm {\boldsymbol{\dot x_{k+1}}}\) is used instead of \(\mathrm {\boldsymbol{\dot x_k}}\) in Equation (38). This ensures stability when the method is applied to linear models with no damping.

4.2. 3rd-order Runge-Kutta-like method

\[\begin{array}{lll} \Delta \boldsymbol{\dot x_1}&=T\boldsymbol{f_k}&\text{(a)}\\\\ \Delta \boldsymbol{x_1}&=T\boldsymbol{\dot x_k}+\frac{1}{4}T^2\boldsymbol{f_k}&\text{(b)}\\\\ \boldsymbol{f'_{k+1}}&=\boldsymbol{f}(\boldsymbol{x_k}+\Delta \boldsymbol{x_1},\boldsymbol{\dot x_k}+\Delta \boldsymbol{\dot x_1},\boldsymbol{\xi _{k+1}})&\text{(c)}\\\\ \Delta \boldsymbol{\dot x_2}&=T\boldsymbol{f'_{k+1}}&\text{(d)}\\\\ \Delta \boldsymbol{x_2}&=T(\boldsymbol{x_k}+\frac{5}{12}\Delta \boldsymbol{x_1})+\frac{1}{3}T^2\boldsymbol{f'_{k+1}}&\text{(e)}\\\\ \boldsymbol{\dot x_{k+1}}&=\boldsymbol{\dot x_k}+\frac{1}{2}(\Delta \boldsymbol{\dot x_1}+\Delta \boldsymbol{\dot x_2})&\text{(f)}\\\\ \boldsymbol{x_{k+1}}&=\boldsymbol{x_k}+\frac{1}{2}(\Delta \boldsymbol{x_1}+\Delta \boldsymbol{x_2})&\text{(g)} \end{array}\]

4.3. Newmark-Beta predictor-corrector method

Predictor (cf. Equation (38))

\[\begin{array}{ll}\boldsymbol{\dot x_{k+1}^{(0)}}=\boldsymbol{\dot x_k}+T\boldsymbol{f_k}\\\\\boldsymbol{x_{k+1}^{(0)}}=\boldsymbol{x_k}+T\boldsymbol{\dot x_{k+1}^{(0)}}\end{array}\]

Corrector:

\[\begin{array}{lll} \boldsymbol{f_{k+1}^{(i)}}&=\boldsymbol{f}(\boldsymbol{x_{k+1}^{(i)}},\boldsymbol{\dot x_{k+1}^{(i)}},\boldsymbol{\xi _{k+1}})&\text{(a)}\\\\ \boldsymbol{\dot x_{k+1}^{(i+1)}}&=\boldsymbol{\dot x_k}+T(1-\gamma )\boldsymbol{f_k}+\gamma \boldsymbol{f_{k+1}^{(i)}}&\text{(b)}\\\\ \boldsymbol{x_{k+1}^{(i+1)}}&=\boldsymbol{x_k}+T\boldsymbol{\dot x_k}+(\frac{1}{2}-\beta )T^2\boldsymbol{f_k}+\beta T^2\boldsymbol{f_{k+1}^{(i)}}&\text{(c)} \end{array}\]

Steps (a) - (c) are repeated a predefined number of times, and as a future extension until

\[\begin{array}{ll}|\boldsymbol{x_{k+1}^{(i+1)}}-\boldsymbol{x_{k+1}^{(1)}}|&<\boldsymbol{\epsilon}\\\\ |\boldsymbol{\dot x_{k+1}^{(i+1)}}-\boldsymbol{\dot x_{k+1}^{(1)}}|&<\displaystyle \frac{\boldsymbol{\epsilon}}{T}\end{array}\]

where \(\mathrm {\boldsymbol{\epsilon}}\) is a vector of specified numbers.

The parameter \(\mathrm {\gamma }\) governs the algorithmic damping in the numerical integration:

  • \(\mathrm {\gamma >\frac{1}{2}}\) gives positive damping

  • \(\mathrm {\gamma =\frac{1}{2}}\) gives no damping

  • \(\mathrm {\gamma <\frac{1}{2}}\) gives negative damping

The parameter \(\mathrm {\beta }\) in Equation (41) should be chosen in the interval \(\mathrm {[0,0.5]}\) . With \(\mathrm {\gamma =1/2}\) , the following known integration methods can be obtained:

  • \(\beta = 0\): Second central difference

  • \(\beta =\frac{1}{12}\): Fox-Goodwins method

  • \(\beta =\frac{1}{6}\): Linear acceleration

  • \(\beta =\frac{1}{4}\): Constant average acceleration (trapezoidal method), unconditionally stable