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, it will experience added inertia due to the acceleration of the surrounding fluid particles. The added inertia force and moment 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 to the added inertia.

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}v}+\boldsymbol{A_{12}\omega }\\\\\boldsymbol{L_A}&=\boldsymbol{A_{21}v}+\boldsymbol{A_{22}\omega }\end{array}\]

where \(\mathrm {\boldsymbol{A_{ij}},i=1,2}\) and \(\mathrm {j=1,2}\) are added-mass tensors with respect to the body origin. 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 inertia force and moment are given by

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

which is completely analogous to Equation (8). For a derivation of this result, (Newman, 1977) or (Milne-Thomson, 1955) can be consulted. 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{v}\times \boldsymbol{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}\]

and the virtual linear and angular momentum

\[\begin{array}{ll}\boldsymbol{P}&=\boldsymbol{P_A}+\boldsymbol{P_B}\\\\\boldsymbol{L}&=\boldsymbol{L_A}+\boldsymbol{L_B}\end{array}\]

the equation of motion becomes

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

where \(\mathrm {\boldsymbol{P}}\) and \(\mathrm {\boldsymbol{L}}\) are obtained from

\[\begin{bmatrix}\boldsymbol{P}\\\boldsymbol{L}\\\end{bmatrix}=\boldsymbol{V}\begin{bmatrix}\boldsymbol{v}\\\boldsymbol{\omega }\end{bmatrix}\]

The accelerations occurring in Equation (16) 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 (18) in acceleration gives:

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

The added-mass force depends on the body’s velocity relative to the fluid. Thus when a homogeneous steady current is present, this relative velocity must be used in the expressions for the added-mass force Equation (11) and force is not influenced by the addition of a constant vector to the body velocity (in spite of the apparent dependency of velocity through the presence of \(\mathrm {\boldsymbol{v}}\) in \(\mathrm {\boldsymbol{P_{B}}}\) and \(\mathrm {\boldsymbol{L_{B}}}\) ). As a consequence, Equation (19) and equations preceding it can be used as they stand, provided \(\mathrm {\boldsymbol{v}}\) is interpreted as the relative velocity.

The assumptions which led to the expressions for added-mass force does not hold in real life. As regards the assumption of ideal-fluid properties, it is commonly assumed that the effect of viscosity can be treated separately and simply added to the ideal-fluid forces. Hence the viscous forces will be contained in the vectors \(\mathrm {\boldsymbol{F_{1}}}\) and \(\mathrm {\boldsymbol{M_{1}}}\) in Equation (19). Care must be taken, however, if drag-force relationships obtained from experiments are used. It is then likely that \(\mathrm {\boldsymbol{M_{_{1}}}}\) will include the ideal-fluid effects which are already taken care of by the term \(\mathrm {\boldsymbol{v}\times \boldsymbol{P_A}}\) in Equation (13) which is included in \(\mathrm {\boldsymbol{v}\times \boldsymbol{P}}\) in Equation (19). In this case, therefore, the appropriate components of \(\mathrm {\boldsymbol{v}\times \boldsymbol{P_A}}\)should be set to zero. Thus, if quadratic current drag coefficients are given for a body, the yaw component of \(\mathrm {\boldsymbol{v}\times \boldsymbol{P_A}}\) is set to zero.

For surface vessels, the assumption of unboundedness of the fluid is violated, as there will be a free surface. Due to wave-making, memory effects will exist, i.e. the hydrodynamic force will not only be a function of the instantaneous velocity and acceleration. In the frequency domain this is known as frequency dependence of added mass. In this case, the high frequency limit of the added mass should be used in Equation (19) and the effect of frequency dependence be included in \(\mathrm {\boldsymbol{F_{1}}}\) and \(\mathrm {\boldsymbol{M_{1}}}\) as convolution integrals.

When there are rigid boundaries in the vicinity of the body, e.g. sea-bottom, walls or other bodies, the results in this chapter does not apply since the fluid forces will depend on the body’s position and orientation and not only on its linear and angular velocities. A general treatment of this problem is complicated. Some special cases in the context of ship steering are treated by (Norrbin, 1971).

When body coordinates are used to represent the vectors on the right-hand side of Equation (19), 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 (24) 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 (19) 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 (19) 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 (29) 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 (30) 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 (24) 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 (19).

  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 (18) 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 (18).

  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 (32).

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 (35). 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 (35))

\[\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 (38) 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