Equations of motion
1. Kinetics for a rigid body moving in water
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 bodyfixed system of axes will be chosen.
The linear and angular momenta are given by:
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 bodyfixed system of reference, Equation (1) can be written:
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:
the vector crossmultiplication occurring in Equation (2) is equivalent to multiplication from the left by the matrix
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
Equation (2) can be written:
and Equation (3) becomes:
The 6 x 6 matrix =
will be referred to as the bodymass 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 addedmass 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.:
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:
where \(\mathrm {\boldsymbol{A_{ij}},i=1,2}\) and \(\mathrm {j=1,2}\) are addedmass 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
will be referred to as the addedmass matrix.
The added inertia force and moment are given by
which is completely analogous to Equation (8). For a derivation of this result, (Newman, 1977) or (MilneThomson, 1955) can be consulted. Equation (13) is seen to agree with the wellknown 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 virtualmass matrix
and the virtual linear and angular momentum
the equation of motion becomes
where \(\mathrm {\boldsymbol{P}}\) and \(\mathrm {\boldsymbol{L}}\) are obtained from
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
Substituting Equation (18) in acceleration gives:
The addedmass 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 addedmass 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 addedmass force does not hold in real life. As regards the assumption of idealfluid properties, it is commonly assumed that the effect of viscosity can be treated separately and simply added to the idealfluid 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 dragforce relationships obtained from experiments are used. It is then likely that \(\mathrm {\boldsymbol{M_{_{1}}}}\) will include the idealfluid 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 wavemaking, 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. seabottom, 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 righthand 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 zaxis. This alters the direction of the body xaxis and yaxis.

Second, the body is rotated by an angle \(\mathrm {\theta }\) about its yaxis. This changes the orientation of the xaxis and the zaxis, but the body yaxis still lies in the global xyplane.

Finally, the body is rotated by an angle \(\mathrm {\phi }\) about its xaxis, altering the orientation of the yaxis and zaxis.
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:
i.e. if a vector is represented by \(\mathrm {\boldsymbol{x^{(B)}}}\) in the body system, the representation in the global system is
Conversely,
since \(\mathrm {\boldsymbol{\Lambda^{1}}=\boldsymbol{\Lambda^T}}\)
3. Position vector. Model formulation
The position vector of the body is defined as
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:
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 statespace formulation, since a wider class of numerical methods for solution is available.
Consider the vectors:
The subvectors \(\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 globalsystem representation, i.e.
Let the vector \(\mathrm {\boldsymbol{\omega }}\) in Equation (19) be expressed as
in the body system. The following relationship can be shown to exist, cf. (Norrbin, 1971):
Equation (29) can be expressed in vector form:
where:
Differentiating Equation (30) with respect to time gives
where
Assuming, \(\mathrm {\boldsymbol{x}}\) ,
\(\mathrm {\boldsymbol{\xi }}\)
and \(\mathrm {\boldsymbol{\dot x}}\) are known,
Equation (24) is realized as follows:

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

\(\mathrm {\boldsymbol{\omega }}\) is calculated from \(\mathrm {\dot \gamma }\) using the inverse of Equation (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 }\) .

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

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

\(\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:
4.1. Modified Euler method
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. 3rdorder RungeKuttalike method
4.3. NewmarkBeta predictorcorrector method
Predictor (cf. Equation (35))
Corrector:
Steps (a)  (c) are repeated a predefined number of times, and as a future extension until
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}\): FoxGoodwins method

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

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