Motion of any general rigid body is governed by conservation of linear and angular momenta (Newton’s second law of motion). The derivation of the dynamic equations of motion for a rigid body can be separated into two steps:
The application of conservation of linear momemtum yields the translational equations of motion (also known as Newton’s equations for a rigid body). In this step the vehicle is assumed to be a point mass whose mass is assumed to be concentrated at the center of gravity \(G\).
The application of conservation of angular momemtum yields the rotational euations of motion (also known as Euler’s equations for a rigid body). In this step the mass distribution of the rigid body plays a role in the dynamics.
Consider a rigid body that is both translating and rotating as shown in Figure 2.1, whose center of gravity is denoted by \(G\) and the BCS origin is denoted by \(O\). The GCS frame is depicted as \(X_0Y_0Z_0\) with origin at \(O_0\) and the BCS frame is depicted as \(X'Y'Z'\) with origin \(O\). The three unit vectors along GCS are denoted as \((\hat{i}_0, \hat{j}_0, \hat{k}_0)\) and the corresponding unit vectors along BCS are denoted as \((\hat{i}, \hat{j}, \hat{k})\). The position vector from GCS origin to BCS origin expressed in GCS is denoted as shown in \(\eqref{eq-bcs-origin-vec-gcs}\). The position vector from GCS origin to center of gravity \(G\) expressed in GCS is denoted as shown in \(\eqref{eq-cog-vec-gcs}\). Let vector from BCS origin to the center of gravity \(G\) expressed in BCS be denoted as shown in \(\eqref{eq-cog-vec-bcs}\).
Let the velocity and acceleration of the center of gravity \(G\) expressed in GCS and BCS be denoted by \((\vec{v}_{0G}, \vec{a}_{0G})\) and \((\vec{v}_{G}, \vec{a}_{G})\) respectively. Similarly the angular velocity of the body expressed in GCS and BCS be denoted by \((\vec{\omega}_0, \vec{\alpha}_0)\) and \((\vec{\omega}, \vec{\alpha})\) respectively. Let the velocity and acceleration of the BCS origin \(O\) expressed in GCS and BCS be denoted by \((\vec{v}_{0O}, \vec{a}_{0O})\) and \((\vec{v}, \vec{a})\) respectively. Note that from the previous chapter, the vectors \(\vec{v}\) and \(\vec{\omega}\) are denoted as shown in \(\eqref{eq-vel-bcs}\) and \(\eqref{eq-ang-vel-bcs}\) respectively.
It can be seen that the position vector \(\vec{r}_{0G}\) can be expressed as shown in \(\eqref{eq-position-vec-triangle-law}\) where \(\boldsymbol{R}\) is the rotation matrix as defined in \(\eqref{eq-rotm-expansion}\).
Note that \(\dot{\vec{r}}_G = 0\) as in the BCS frame, the location of center of gravity \(G\) will not vary with time due to the rigid body assumption. The vector \(\vec{r}_G\) can be expressed from \(\eqref{eq-position-vec-triangle-law}\) as shown in \(\eqref{eq-rG-recast}\).
The product \(\dot{\boldsymbol{R}} \boldsymbol{R}^T\) can be evaluated using symbolic python (sympy package) as seen in the code give below and the final expression in summarized in a compact form in \(\eqref{eq-RdotRT-compact-1}\). Note that \(s_1 = \sin(\phi)\), \(c_1 = \cos(\phi)\), \(s_2 = \sin(\theta)\), \(c_2=\cos(\theta)\), \(s_3 = \sin(\psi)\) and \(c_3 = \cos(\psi)\).
Code
import sympy as sp# Define symbols for angles and their time derivativesphi, theta, psi = sp.symbols('phi theta psi', cls=sp.Function)t = sp.symbols('t')phi = phi(t)theta = theta(t)psi = psi(t)# Define derivatives of angles with respect to timephi_dot = sp.diff(phi, t)theta_dot = sp.diff(theta, t)psi_dot = sp.diff(psi, t)# Define trigonometric functionss1, c1 = sp.sin(phi), sp.cos(phi)s2, c2 = sp.sin(theta), sp.cos(theta)s3, c3 = sp.sin(psi), sp.cos(psi)# Define the rotation matrix RR = sp.Matrix([ [c2*c3, -c1*s3 + s1*s2*c3, s1*s3 + c1*s2*c3], [c2*s3, c1*c3 + s1*s2*s3, -s1*c3 + c1*s2*s3], [-s2, s1*c2, c1*c2]])# Compute the time derivative of RR_dot = R.diff(t)# Compute R^TR_T = R.transpose()# Compute R_dot * R_Tresult = R_dot * R_T# Simplify the resultsimplified_result = sp.simplify(result)simplified_result
From the previous chapter it is seen that the angular velocity of the body \(\vec{\omega}\) is related to the Euler rates by \(\eqref{eq-ang-vel-bcs-repeat}\).
The angular velocity of the body expressed in GCS will be given by \(\eqref{eq-ang-vel-gcs}\). Note that the following code has been used to evaluate the final expression shown in \(\eqref{eq-ang-vel-gcs}\).
Code
import sympy as sp# Define symbols for angles and their time derivativesphi, theta, psi = sp.symbols('phi theta psi', cls=sp.Function)t = sp.symbols('t')phi = phi(t)theta = theta(t)psi = psi(t)# Define derivatives of angles with respect to timephi_dot = sp.diff(phi, t)theta_dot = sp.diff(theta, t)psi_dot = sp.diff(psi, t)# Define trigonometric functionss1, c1 = sp.sin(phi), sp.cos(phi)s2, c2 = sp.sin(theta), sp.cos(theta)s3, c3 = sp.sin(psi), sp.cos(psi)# Define the rotation matrix RR = sp.Matrix([ [c2*c3, -c1*s3 + s1*s2*c3, s1*s3 + c1*s2*c3], [c2*s3, c1*c3 + s1*s2*s3, -s1*c3 + c1*s2*s3], [-s2, s1*c2, c1*c2]])# Define angular velocity omega in BCSomega = sp.Matrix([[phi_dot - s2*psi_dot], [c1*theta_dot + s1*c2*psi_dot], [-s1*theta_dot + c1*c2*psi_dot]])# Angular velocity in GCS omega_0omega_0 = R * omega# Simplify the resultsimplified_result = sp.simplify(omega_0)simplified_result
It can be seen that the right hand side of \(\eqref{eq-RdotRT-compact-1}\) is a skew symmetric matrix constructed from \(\vec{\omega}_0\) shown in \(\eqref{eq-ang-vel-gcs}\). Using the property of expressing cross product as a matrix multiplication, \(\eqref{eq-vel-triangle-law-1}\) can be written as shown in \(\eqref{eq-vel-triangle-law-2}\).
Substituting \((\vec{v}_{0G} - \vec{v}_{0O})\) from \(\eqref{eq-vel-triangle-law-2}\) into \(\eqref{eq-acc-triangle-law-1}\) yields \(\eqref{eq-acc-triangle-law-2}\).
Newton’s second law for a point mass can be written as \(\eqref{eq-newtons-2nd-law}\) where \(\vec{F}_{0}\) is the force acting on the body expressed in GCS. \(m\) denotes the mass of the rigid body.
\[\begin{align}
m \vec{a}_{0G} = \vec{F}_{0}
\label{eq-newtons-2nd-law}
\end{align}\]
Sustituting the acceleration of \(G\) expressed in GCS from \(\eqref{eq-acc-triangle-law-2}\) leads to \(\eqref{eq-newtons-eom}\).
Notice that all the quantities in \(\eqref{eq-newtons-eom}\) are expressed in GCS frame. Premultiplying \(\boldsymbol{R}^T\) on both sides of \(\eqref{eq-newtons-eom}\) leads to translational equations of motion expressed in BCS frame as shown in \(\eqref{eq-newtons-eom-1}\).
Note that \(\boldsymbol{R}^T(\vec{\alpha}_0 \times (\vec{r}_{0G} - \vec{r}_{0O})) = \vec{\alpha} \times \vec{r}_{G}\) and \(\boldsymbol{R}^T(\vec{\omega}_0 \times (\vec{\omega}_0 \times (\vec{r}_{0G} - \vec{r}_{0O}))) = \vec{\omega} \times (\vec{\omega} \times \vec{r}_{G})\). Transforming vectors to BCS and taking a cross product is equivalent to taking a cross product of vectors in GCS and transforming the resultant vector to BCS. \(\vec{F} = X \hat{i} + Y \hat{j} + Z \hat{k}\) in \(\eqref{eq-newtons-eom-1}\) represents the force acting of the rigid body expressed in BCS frame.
In BCS \(\vec{a} \neq \dot{\vec{v}}\) as BCS is a non-inertial frame of reference. Differentiation of velocity will yield acceleration only in an inertial frame. Instead \(\vec{a}\) can be expressed in terms of \(\dot{\vec{v}}\) and \(\vec{v}\) as shown in \(\eqref{eq-acc-vel-der-bcs}\).
Substituting \(\eqref{eq-acc-vel-der-bcs}\) into \(\eqref{eq-newtons-eom-1}\) results in the translational equations of motion in the vector form and are shown in \(\eqref{eq-newtons-eom-2}\).
Let the angular momentum of the vessel about the center of gravity \(G\) expressed in GCS be denoted by \(\vec{L}_{0G}\). Let the angular momentum of the vessel about the center of gravity \(G\) expressed in BCS be denoted by \(\vec{L}_G\). Let the external moment on the vessel about the center of gravity \(G\) expressed in GCS be denoted by \(\vec{M}_{0G}\). The conservation of angular momentum for the vessel is shown by \(\eqref{eq-ang-momentum-conservation}\).
It is well known from rigid body mechanics that the angular momentum of the body in BCS \(\vec{L}_{G}\) can be expressed as the product of inertia tensor and the angular velocity as shown in \(\eqref{eq-ang-mom-inertia}\).
\(I_x^G\), \(I_y^G\) and \(I_z^G\) are the second mass moments of inertia. \(I_{xy}^G\), \(I_{xz}^G\), \(I_{yz}^G\), \(I_{yx}^G\), \(I_{zx}^G\) and \(I_{zy}^G\) are cross mass moments of inertia. The terms in the inertia tensor \(\boldsymbol{I}_G\) are defined with respect to a coordinate system that is parallel to the BCS but with its origin at the center of gravity \(G\).
The second mass moment of inertia \(I_x^G\) is defined as shown in \(\eqref{eq-IxG}\) where \(\mu\) represents the mass density of the body and the volume integral is performed over the volume \(V\) of the rigid body. Note that \((x - x_G)\), \((y - y_G)\) and \((z - z_G)\) represent the distance between the elemental volume \(dV\) and the center of gravity of the rigid body along directions parallel to the y and z axes of the BCS respectively. \(I_y^G\) and \(I_z^G\) are also defined in a similar manner.
\[\begin{align}
I_x^G = I_{yy}^G + I_{zz}^G = \iiint_V \mu \left[(y - y_G)^2 + (z - z_G)^2\right] dV
\label{eq-IxG}
\end{align}\]
The cross mass moment of inertia \(I_{xy}^G\) is defined as shown in \(\eqref{eq-IxyG}\). The other terms \(I_{xz}^G\), \(I_{yz}^G\), \(I_{yx}^G\), \(I_{zx}^G\) and \(I_{zy}^G\) are defined in a similar manner.
Note that for a rigid body \(\boldsymbol{I}_G\) is time invariant. The advantage of formulating the equations of motion in BCS is that the mass properties of the body are time invariant. On the other hand, when the equations of motion are formulated in GCS, the inertia tensor would be a time varying quantity.
Substituting \(\eqref{eq-ang-mom-inertia}\) into \(\eqref{eq-ang-momentum-conservation-2}\) results in \(\eqref{eq-rotational-eqn-cog}\).
Let \(\vec{M} = K \hat{i} + M \hat{j} + N \hat{k}\) denote the external moment acting on the vessel about BCS origin \(O\) expressed in BCS. \(\vec{M}\) and \(\vec{M}_{0G}\) are related as shown in \(\eqref{eq-ext-moment}\) where \(\vec{F}\) is the external force acting on the body expressed in BCS.
Some of the terms can be simplified. Taking the first two terms together and applying the properties of vector cross product and expressing cross products as matrix multiplications yeilds \(\eqref{eq-inertia-simplify-1}\) and \(\eqref{eq-inertia-simplify-1a}\) respectively.
Note that \(\boldsymbol{I}_O\) in the last step is the inertia tensor defined with respect to BCS and is related to \(\boldsymbol{I}_G\) through the parallel axes theorem shown in \(\eqref{eq-parallel-axes-theorem}\).
\[\begin{align}
\boldsymbol{I}_O = \boldsymbol{I}_G - m \boldsymbol{S}(\vec{r}_G) \boldsymbol{S}(\vec{r}_G)
\label{eq-parallel-axes-theorem}
\end{align}\]
Similarly, the third and fourth terms of \(\eqref{eq-euler-eom-2}\) can also be simplified using the Jacobi identity for triple cross products as shown in \(\eqref{eq-inertia-simplify-2}\), \(\eqref{eq-inertia-simplify-3}\) and \(\eqref{eq-inertia-simplify-4}\).
Substituting final expressions from \(\eqref{eq-inertia-simplify-1a}\) and \(\eqref{eq-inertia-simplify-4}\) into \(\eqref{eq-euler-eom-2}\) yeilds the rotational equations in the vector form as shown in \(\eqref{eq-euler-eom-3}\).
When expanded, these form six equations, with one corresponding to each degree of freedom and are shown in (\(\ref{eq-eom-expanded-1}\) - \(\ref{eq-eom-expanded-6}\)).
These equations can be represented in a compact matrix form as shown in \(\eqref{eq-eom-matrix-form}\), where \(\boldsymbol{M}_{RB}\) is a \(6 \times 6\) rigid body mass matrix, \(\boldsymbol{C}_{RB}\left(\{v\}\right)\) is the \(6 \times 6\) velocity dependent Coriolis matrix and \(\{\tau_{RB}\}\) is the \(6 \times 1\) external force and moment vector.
The rigid body mass matrix \(\boldsymbol{M}_{RB}\) is shown in \(\eqref{eq-mass-matrix}\), where \(\boldsymbol{I}_{3 \times 3}\) is a \(3 \times 3\) identity matrix.
The Coriolis matrix can be parameterized in terms of the mass matrix as shown in \(\eqref{eq-coriolis-mass-param}\), where \(\boldsymbol{0}_{3 \times 3}\) is a \(3 \times 3\) matrix of zeros. Note that \(\boldsymbol{M}_{11} = m\boldsymbol{I}_{3 \times 3}\), \(\boldsymbol{M}_{12} = -\boldsymbol{M}_{21} = -m \boldsymbol{S}(\vec{r}_G)\) and \(\boldsymbol{M}_{22} = \boldsymbol{I}_O\).
Defining a state vector of size \(12 \times 1\) as shown in \(\eqref{eq-state}\), the kinematic and dynamic equations can be expressed together as shown in \(\eqref{eq-kinematics-dynamics}\).
\[\begin{align}
\{x\} = \begin{bmatrix}
u & v & w & p & q & r & x_{0O} & y_{0O} & z_{0O} & \phi & \theta & \psi
\end{bmatrix}^T
\label{eq-state}
\end{align}\]
If the external forces and moments \(\{\tau_{RB}\}\) are known, then the velocities and pose of the body can be calculated by solving the matrix differential equation shown in \(\eqref{eq-kinematics-dynamics}\).
2.4 Exercises
For a rigid body with mass \(m = 100\) kg and center of gravity located at \(\vec{r}_G = [0.1, -0.15, 0.2]^T\) m in BCS, calculate both the upper right and lower left \(3 \times 3\) blocks of the rigid body mass matrix \(\boldsymbol{M}_{RB}\). Verify that they are negative transposes of each other.
Given a rigid body with principal moments of inertia \(I_x^G = 20\) kg⋅m², \(I_y^G = 25\) kg⋅m², \(I_z^G = 30\) kg⋅m² about its center of gravity, and products of inertia \(I_{xy}^G = 2\), \(I_{xz}^G = -3\), \(I_{yz}^G = 1\) kg⋅m², construct the complete inertia tensor \(\boldsymbol{I}_G\). Is the resulting matrix symmetric?
A rigid body has angular velocity \(\vec{\omega} = [0.1, 0.2, 0.3]^T\) rad/s. Calculate the skew-symmetric matrix \(\boldsymbol{S}(\vec{\omega})\) and verify the properties:
Determine the additional inertia terms needed to transform \(\boldsymbol{I}_G\) to \(\boldsymbol{I}_O\),
Show that the resulting matrix is symmetric.
A rigid body has linear velocity \(\vec{v} = [1, 2, -1]^T\) m/s and angular velocity \(\vec{\omega} = [0.1, 0.2, -0.15]^T\) rad/s. Calculate:
The acceleration term \(\vec{\omega} \times \vec{v}\) using cross product,
The same term using the skew-symmetric matrix \(\boldsymbol{S}(\vec{\omega})\),
Verify both methods give the same result.
A rigid body has mass \(m = 200\) kg and center of gravity at \(\vec{r}_G = [0.2, -0.1, 0.15]^T\) m in BCS. If it experiences forces \(\vec{F}_1 = [100, -50, 25]^T\) N at point \(\vec{r}_1 = [0, 0, 0]^T\) m and \(\vec{F}_2 = [-30, 40, -20]^T\) N at point \(\vec{r}_2 = [0.1, 0.1, 0.1]^T\) m (both in BCS), calculate:
The total force,
The moment about O,
The moment about G.
For a rigid body with angular velocity \(\vec{\omega} = [0.5, -0.3, 0.2]^T\) rad/s and inertia tensor \(\boldsymbol{I}_G\) with principal moments \([20, 25, 30]\) kg⋅m² and products of inertia \(I_{xy}^G = I_{yz}^G = I_{xz}^G = 5\) kg⋅m², calculate:
The rotational kinetic energy \(\frac{1}{2}\vec{\omega}^T\boldsymbol{I}_G\vec{\omega}\).
Given a rigid body with mass \(m = 300\) kg, \(\vec{r}_G = [0.25, -0.15, 0.1]^T\) m, and \(\boldsymbol{I}_G = \text{diag}(40, 45, 50)\) kg⋅m², calculate:
The complete \(6 \times 6\) rigid body mass matrix \(\boldsymbol{M}_{RB}\),
Show that it is symmetric,
Calculate its determinant.
A rigid body has linear velocity \(\vec{v} = [2, 1, -1]^T\) m/s, angular velocity \(\vec{\omega} = [0.1, -0.2, 0.15]^T\) rad/s, and center of gravity at \(\vec{r}_G = [1, -0.5, 0.3]^T\) m. Calculate:
The velocity of G in BCS,
The acceleration terms \(\vec{\omega} \times \vec{v}\) and \(\vec{\omega} \times (\vec{\omega} \times \vec{r}_G)\),
The total acceleration of G if \(\dot{\vec{v}} = [0.1, 0.2, -0.1]^T\) m/s² and \(\dot{\vec{\omega}} = [0, 0.1, 0]^T\) rad/s².
For a rigid body with \(\vec{\omega} = [0.2, -0.1, 0.3]^T\) rad/s and \(\vec{r}_G = [0.5, -0.2, 0.1]^T\) m, calculate:
The centripetal acceleration term \(\vec{\omega} \times (\vec{\omega} \times \vec{r}_G)\) using cross products,
The same term using skew-symmetric matrices,
Show that both methods yield identical results.
A rigid body has mass \(m = 400\) kg, center of gravity at \(\vec{r}_G = [0.3, -0.1, 0.2]^T\) m, and inertia tensor \(\boldsymbol{I}_G\) with principal moments \([50, 60, 45]\) kg⋅m² and products of inertia \(I_{xy}^G = I_{yz}^G = I_{xz}^G = 8\) kg⋅m². If it experiences force \(\vec{F} = [200, -100, 50]^T\) N and moment \(\vec{M} = [20, -10, 5]^T\) N⋅m at O, with initial velocities \(\vec{v} = [1, -0.5, 0.2]^T\) m/s and \(\vec{\omega} = [0.1, 0, -0.1]^T\) rad/s, calculate:
The complete mass matrix \(\boldsymbol{M}_{RB}\),
All acceleration-dependent terms,
Solve for \(\dot{\vec{v}}\) and \(\dot{\vec{\omega}}\).
For a rigid body of mass \(m = 100\) kg, center of gravity at \(\vec{r}_G = [0.1, -0.1, 0.2]^T\) m, and inertia tensor \[\begin{align*}
\boldsymbol{I}_G = \begin{bmatrix}
30 & 4 & -2 \\
4 & 35 & 3 \\
-2 & 3 & 40
\end{bmatrix} \text{ kg⋅m²}
\end{align*}\] angular velocity \(\vec{\omega} = [0.1, 0.2, 0.3]^T\) rad/s, and linear velocity \(\vec{v} = [1, -0.5, 0.2]^T\) m/s:
Derive the complete expression for \(\boldsymbol{C}_{RB}(\{v\})\{v\}\),
Show that \(\{v\}^T\boldsymbol{C}_{RB}(\{v\})\{v\} = 0\),
Explain the physical significance of this property.
Given a rigid body with mass \(m = 600\) kg, \(\vec{r}_G = [0.2, 0.1, -0.1]^T\) m, and inertia tensor \[\begin{align*}
\boldsymbol{I}_G = \begin{bmatrix}
500 & -25 & -30 \\
-25 & 600 & 35 \\
-30 & 35 & 650
\end{bmatrix} \text{ kg⋅m²}
\end{align*}\]
Derive the complete expression for the Coriolis matrix \(\boldsymbol{C}_{RB}(\{v\})\),
Show that \(\dot{\boldsymbol{M}}_{RB} - 2\boldsymbol{C}_{RB}(\{v\})\) is skew-symmetric,
Explain why this property is important for energy conservation.
Hint: What should be the rate of change of kinetic energy (\(\dot{T}\)) for a rigid body that is not subject to any external forces (i.e. \(\{\tau_{RB}\}=0\))? Remember that the total kinetic energy of a rigid body is given by \(T = \frac{1}{2} \vec{v}^T \boldsymbol{M}_{RB} \vec{v}\)
A marine vehicle has the following inertia tensors measured about its center of gravity (G) and about the origin of BCS (O): \[\begin{align*}
\boldsymbol{I}_G = \begin{bmatrix}
400 & -20 & -15 \\
-20 & 450 & 25 \\
-15 & 25 & 500
\end{bmatrix} \text{ kg⋅m²}, \quad
\boldsymbol{I}_O = \begin{bmatrix}
410 & -10 & -25 \\
-10 & 475 & 30 \\
-25 & 30 & 525
\end{bmatrix} \text{ kg⋅m²}
\end{align*}\] If the mass of the vehicle is \(m = 500\) kg:
Determine the coordinates of the center of gravity \(\vec{r}_G = [x_G, y_G, z_G]^T\)
Calculate \(\boldsymbol{S}(\vec{r}_G)\boldsymbol{S}(\vec{r}_G)\) and verify it matches with \(\frac{1}{m}(\boldsymbol{I}_O - \boldsymbol{I}_G)\)
Programming Exercise: A marine vehicle with mass \(m = 1000\) kg, center of gravity at \(\vec{r}_G = [0.2, 0, 0.1]^T\) m, and inertia tensor \[\begin{align*}
\boldsymbol{I}_G = \begin{bmatrix}
800 & 0 & -50 \\
0 & 1200 & 0 \\
-50 & 0 & 1000
\end{bmatrix} \text{ kg⋅m²}
\end{align*}\] starts from rest with initial orientation \(\phi = 20^{\circ}\), \(\theta = 30^{\circ}\), \(\psi = 45^{\circ}\) and experiences time-varying forces and moments for \(25\) seconds:
\(\vec{F}(t) = [500\cos\left(\frac{\pi t}{50}\right), 200\sin\left(\frac{\pi t}{50}\right), -100]^T\) N
The skew-symmetry of \(\dot{\boldsymbol{M}}_{RB} - 2\boldsymbol{C}_{RB}\) ensures that the rate of change of kinetic energy is zero when there are no external forces or moments.