1  Rigid Body Kinematics

Motion of a marine vehicle (assumed to be rigid) requires the knowledge of the pose and velocities of the vehicle as a function of time. Vehicle motion is governed by both kinematics and dynamics. Kinematics refers to the study of the motion of a rigid body without considering the forces and moments that cause the motion. The study of the forces and moments acting on the vehicle that cause its motion is known as dynamics. This first chapter will focus on the kinematics of a rigid body.

1.1 Coordinate Systems

The pose and velocity of a rigid body is always defined with respect a coordinate frame of reference. Therefore, coordinate systems need to be defined before the description of motion of a rigid body is possible. Two right hand coordinate systems are defined as follows:

  • Global Coordinate System (GCS): This is an inertial coordinate system fixed to the Earth. The usual convention when considering vehicle’s maneuvering motions is to consider a North-East-Down (NED) frame of reference with origin on the surface of the earth.
  • Body Coordinate System (BCS): This is a coordinate system fixed with the vehicle. As the body undergoes motion the BCS also moves with it. Thus BCS is not an inertial frame of reference. Usual convention is to choose the x-axis of BCS to be along the length of the vehicle, z-axis pointing downward and y-axis following the right hand rule. For a ship floating in even keel, this will result in x-axis along the longitudinal direction (from aft to forward) of the ship, y-axis pointing towards starboard direction and z-axis pointing vertically down. The origin of the BCS may be chosen either to be at the center of gravity or a geometric point on the vessel. In general, the center of gravity of the vehicle would change depending on the vessel’s loading condition (fully loaded or ballast loaded), consumption of fuel during the voyage and due to loading or unloading of cargo at ports. Thus, it is better to choose a geometric point on the vehicle that does not change with time as the origin of BCS. The intersection of midship section, centerline and fully loaded waterline is a good choice for the origin of BCS. For underwater vehicles the geometric center of the vehicle may be chosen as the BCS origin.

The above described choice of BCS and GCS origins and their orientations are the generally used convention. One can choose the origins and frame orientations arbitrarily as long as they are maintained consistently. For ease of understanding, the convention described above will be followed in this text.

Both coordinate systems are illustrated in Figure 1.1. The origin of GCS is denoted as \(O_0\) and its axes are denoted as \(X_0Y_0Z_0\). The coordinates of a point in GCS frame will be denoted by \((x_0, y_0, z_0)\). The origin of the BCS frame is denoted by \(O'\) and its axes are denoted as \(X'Y'Z'\). The coordinates of a point expressed in BCS frame will be denoted by \((x, y, z)\).

Figure 1.1: Global Coordinate System (GCS) and Body Coordinate System (BCS)

The unit vectors along the x, y and z axes of GCS are denoted as \(\hat{i}_0\), \(\hat{j}_0\) and \(\hat{k}_0\) respectively. Similarly the unit vectors along the x, y and z axes of the BCS are denoted as \(\hat{i}\), \(\hat{j}\) and \(\hat{k}\) respectively.

1.2 Pose of a Rigid Body

The pose of a rigid body can be described in terms of the above defined coordinate systems. The translation of the body is described as the position vector from the GCS origin to the BCS origin expressed in the GCS frame. In Figure 1.1 we can see that position vector of the origin of BCS in the GCS frame is given by \(\eqref{eq-gcsbcs-posvec}\).

\[\begin{align} \vec{r}_0 = x_0 \hat{i}_0 + y_0 \hat{j}_0 + z_0 \hat{k}_0 \label{eq-gcsbcs-posvec} \end{align}\]

The three parameters \(x_0\), \(y_0\) and \(z_0\) describe the translation of the body. However, as the BCS moves with the body, it may be in a different orientation as compared to the GCS at any particular instant. This means that the axes of BCS need not be parallel to the axes of GCS. Thus, the orientation also needs to be known to characterize the motion of the vehicle completely. The orientation can be described using one of the following ways:

Euler angle representation: The orientation between two coordinate frames is expressed as a sequence of three rotations in a specified order. Unlike translations, rotations are not commutative and thus a specific order needs to be maintained. The non-cummutative property of rotations can be seen from Figure 1.2 where two \(90^\circ\) rotations are applied in two different orders. The first approach applies a \(90^\circ\) rotation about the x-axis followed by a \(90^\circ\) rotation about the new z-axis. The second approach interchanges the rotation order - a \(90^\circ\) rotation about z-axis is applied first and then a \(90^\circ\) rotation is applied about the new x-axis. The orientations of the new coordinate frames are completely different when the order of rotations is swapped. The Python code used to generate the results showns in Figure 1.2 is shown below. Euler angles will be described in detail in Section 1.3.

Code 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x_vec = np.array([
    [0.0, 0.0, 0.0],
    [1.0, 0.0, 0.0]
]).T

y_vec = np.array([
    [0.0, 0.0, 0.0],
    [0.0, 1.0, 0.0]
]).T

z_vec = np.array([
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 1.0]
]).T


psi = np.pi/2
phi = np.pi/2

Rz = np.array([
    [np.cos(psi), -np.sin(psi), 0.0],
    [np.sin(psi), np.cos(psi), 0.0],
    [0.0, 0.0, 1.0]
])

Rx = np.array([
    [1.0, 0.0, 0.0],
    [0.0, np.cos(phi), -np.sin(phi)],
    [0.0, np.sin(phi), np.cos(phi)]
])

R1 = Rz @ Rx
R2 = Rx @ Rz

x_vec1 = R1 @ x_vec
y_vec1 = R1 @ y_vec
z_vec1 = R1 @ z_vec

x_vec2 = R2 @ x_vec
y_vec2 = R2 @ y_vec
z_vec2 = R2 @ z_vec

fig = plt.figure()
ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122, projection='3d')

ax1.plot(x_vec[0], x_vec[1], x_vec[2], 'r', linewidth=5)
ax1.plot(y_vec[0], y_vec[1], y_vec[2], 'g', linewidth=5)
ax1.plot(z_vec[0], z_vec[1], z_vec[2], 'b', linewidth=5)

ax1.plot(x_vec1[0], x_vec1[1], x_vec1[2], 'c:', linewidth=5)
ax1.plot(y_vec1[0], y_vec1[1], y_vec1[2], 'm:', linewidth=5)
ax1.plot(z_vec1[0], z_vec1[1], z_vec1[2], 'y:', linewidth=5)

ax2.plot(x_vec[0], x_vec[1], x_vec[2], 'r', linewidth=5)
ax2.plot(y_vec[0], y_vec[1], y_vec[2], 'g', linewidth=5)
ax2.plot(z_vec[0], z_vec[1], z_vec[2], 'b', linewidth=5)

ax2.plot(x_vec2[0], x_vec2[1], x_vec2[2], 'c:', linewidth=5)
ax2.plot(y_vec2[0], y_vec2[1], y_vec2[2], 'm:', linewidth=5)
ax2.plot(z_vec2[0], z_vec2[1], z_vec2[2], 'y:', linewidth=5)

ax1.set_title("Rotation Order Z-X")
ax1.set_xlabel("x-axis")
ax1.set_ylabel("y-axis")
ax1.set_zlabel("z-axis")

ax2.set_title("Rotation Order X-Z")
ax2.set_xlabel("x-axis")
ax2.set_ylabel("y-axis")
ax2.set_zlabel("z-axis")

# plt.tight_layout()
plt.show()
Figure 1.2: The solid line represents the GCS orientation and the dotted line represents the BCS oreintation. For GCS, x-axis is shown in red, y-axis in green and z-axis in blue. For BCS, x-axis is shown in cyan, y-axis in magenta and z-axis in yellow. Applying 90 degrees rotation about z-axis followed by 90 degrees rotation about x-axis leads the orientation shown in the right subplot. Interchanging the rotation order lead to the orientation shown in left subplot.

Axis-angle representation: The orientation between the two coordinate frames (GCS and BCS) can be described as a single rotation about an axis of rotation that results in GCS being oriented along the BCS. The rotation is parametrized by the unit vector \(\hat{a}\) along the axis about which the rotation is performed and the magnitude of the angle of rotation \(\Theta\). This parametrization has three independent parameters. The axis being a unit vector has two independent parameters with the third parameter being constrained by the property that the vector \(\hat{a}\) is a unit vector. The third independent parameter is the angle of rotation \(\Theta\).

Quaternion: In this parametrization, the orientation between two coordinate frames (GCS and BCS) is characterized by four parameters. The parameters of a quaternion are closely related to the axis-angle representation parameters \(\Theta\) and \(\hat{a}\). The first parameter of a quaternion is \(\cos(\frac{\Theta}{2})\). The other three parameters are the components of the vector \(\sin(\frac{\Theta}{2}) \hat{a}\). The Quaternions have the primary advantage of being numerically stable, compact and efficient as compared to the other rotation formalisms. However, in contrast to the Euler angles or axis-angle representation they have one extra parameter. Quaternions have found applications in several fields including robotics, space navigation, computer vision, flight dynamics etc.

Rotation Matrix: The orientation between the two coordinate frames can also be characterized by a \(3 \times 3\) rotation matrix that has \(9\) parameters. The rotation matrix projects the vectors from one frame into another. The rotation matrix \(\boldsymbol{R}\) to project vectors from BCS to GCS is given by \(\eqref{eq-rotm-def}\). The rotation matrix to project vectors from GCS to BCS is given by the transpose of the rotation matrix denoted by \(\boldsymbol{R}^T\).

\[\begin{align} \boldsymbol{R} = \begin{bmatrix} \hat{i}.\hat{i}_0 & \hat{j}.\hat{i}_0 & \hat{k}.\hat{i}_0 \\ \hat{i}.\hat{j}_0 & \hat{j}.\hat{j}_0 & \hat{k}.\hat{j}_0 \\ \hat{i}.\hat{k}_0 & \hat{j}.\hat{k}_0 & \hat{k}.\hat{k}_0 \end{bmatrix} \label{eq-rotm-def} \end{align}\]

Consider a vector \(\vec{p}\) expressed in BCS. The same vector expressed in GCS \(\vec{p}_0\) is then given by premultiplying the vector \(\vec{p}\) by the rotation matrix \(\boldsymbol{R}\) as shown in \(\eqref{eq-rotm-bcs-to-gcs}\).

\[\begin{align} \vec{p}_0 = \boldsymbol{R} \vec{p} \label{eq-rotm-bcs-to-gcs} \end{align}\]

Some of the other formalisms of rotation in 3D space can be found here.

Each of the rotation formalisms have their own advantages and disadvantages. The rotation matrix has more parameters than any of the other formulations. Quaternions, on the other hand, represent the most compact formulation and are also numerically stable. However, they are difficult to intuitively appreciate. On the contrary, the Euler angles are most intuitive to understand and only have \(3\) parameters. However, the Euler angle formulation suffers from singularity at two specific orientations. This phenomenon is known as gimbal lock and will be explored more in Section 1.3.1. For marine vehicles like ships, if the order of rotation of the Euler angles is chosen well, the gimbal lock orientation is not encountered in regular operations. Thus for ships, Euler angles have become a popular approach as they are intuitive to understand.

1.3 Euler Angle Representation

Let us consider the GCS and BCS origins to be non-coincident. In order to transform the GCS into BCS, we first translate the GCS origin to the BCS origin. This intermediate coordinate system is labelled as \(X_3Y_3Z_3\) frame. Note that \(X_3Y_3Z_3\) is parallel to GCS (\(X_0Y_0Z_0\) frame). We will assume a ZYX order of rotation from \(X_3Y_3Z_3\) frame to reach BCS (\(X'Y'Z'\) frame).

\(X_3Y_3Z_3\) is now rotated by an angle \(\psi\) about the \(Z_3\) axis to yield a new coordinate system \(X_2Y_2Z_2\). Let us consider a vector \(\vec{r}_3 = [x_3 \quad y_3 \quad z_3]^T\) expressed in \(X_3Y_3Z_3\) frame. The same vector when expressed in \(X_2Y_2Z_2\) frame is denoted by \(\vec{r}_2 = [x_2 \quad y_2 \quad z_2]^T\). Since the rotation has happened along the \(Z_3\) axis, \(Z_2\) axis is coincident with \(Z_3\) axis. Therefore, the projection of the vector \(\vec{r}_3\) along \(Z_3\) axis and the projection of vector \(\vec{r}_2\) along \(Z_2\) axis is the same. Hence, \(z_2\) must be equal to \(z_3\).

Figure 1.3: Rotation from \(X_3Y_3Z_3\) to \(X_2Y_2Z_2\)

Figure 1.3 shows the rotation from \(X_3Y_3Z_3\) frame to \(X_2Y_2Z_2\) frame. Note that the \(Z_3\) and \(Z_2\) axes are coincident and are pointed into the plane. The positive sign convention for \(\psi\) follows the right hand rule about the \(Z_3\) (or \(Z_2\)) axis and is thus denoted in the clockwise direction in Figure 1.3. From Figure 1.3 it can be seen that the \(x_3\) and \(y_3\) can be expressed in terms of \(x_2\), \(y_2\) and \(\psi\) as shown in \(\eqref{eq-z3-rotm}\) and \(\eqref{eq-z3-rotm-expanded}\).

\[\begin{align} \vec{r}_3 = \boldsymbol{R}_{z,\psi} \vec{r}_2 \label{eq-z3-rotm} \end{align}\]

\[\begin{align} \begin{bmatrix} x_3 \\ y_3 \\ z_3 \end{bmatrix} = \begin{bmatrix} \cos(\psi) & -\sin(\psi) & 0 \\ \sin(\psi) & \cos(\psi) & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_2 \\ y_2 \\ z_2 \end{bmatrix} \label{eq-z3-rotm-expanded} \end{align}\]

Now \(X_2Y_2Z_2\) frame is further rotated by angle \(\theta\) about the \(Y_2\) axis to obtain the intermediate frame \(X_1Y_1Z_1\) frame. Let the \(\vec{r}_2\) vector when expressed in \(X_1Y_1Z_1\) frame be denoted by \(\vec{r}_1 = [x_1 \quad y_1 \quad z_1]^T\). Since the rotation has happened along the \(Y_2\) axis, \(Y_1\) axis is coincident with \(Y_2\) axis. Therefore, the projection of the vector \(\vec{r}_2\) along \(Y_2\) axis and the projection of vector \(\vec{r}_1\) along \(Y_1\) axis is the same. Hence, \(y_1\) must be equal to \(y_2\).

Figure 1.4: Rotation from \(X_2Y_2Z_2\) to \(X_1Y_1Z_1\)

Figure 1.4 shows the rotation from \(X_2Y_2Z_2\) frame to \(X_1Y_1Z_1\) frame. Note that the \(Y_2\) and \(Y_1\) axes are coincident and are pointed out of the plane. The positive sign convention for \(\theta\) follows the right hand rule about the \(Y_2\) (or \(Y_1\)) axis and is thus denoted in the anti-clockwise direction in Figure 1.4. From Figure 1.4 it can be seen that the \(x_2\) and \(z_2\) can be expressed in terms of \(x_1\), \(z_1\) and \(\theta\) as shown in \(\eqref{eq-y2-rotm}\) and \(\eqref{eq-y2-rotm-expanded}\).

\[\begin{align} \vec{r}_2 = \boldsymbol{R}_{y,\theta} \vec{r}_1 \label{eq-y2-rotm} \end{align}\]

\[\begin{align} \begin{bmatrix} x_2 \\ y_2 \\ z_2 \end{bmatrix} = \begin{bmatrix} \cos(\theta) & 0 & \sin(\theta) \\ 0 & 1 & 0 \\ -\sin(\theta) & 0 & \cos(\theta) \end{bmatrix} \begin{bmatrix} x_1 \\ y_1 \\ z_1 \end{bmatrix} \label{eq-y2-rotm-expanded} \end{align}\]

Now \(X_1Y_1Z_1\) frame is finally rotated by angle \(\phi\) about the \(X_1\) axis to obtain the BCS frame \(X'Y'Z'\). Let the \(\vec{r}_1\) vector when expressed in BCS frame be denoted by \(\vec{r} = [x \quad y \quad z]^T\). Since the rotation has happened along the \(X_1\) axis, \(X'\) axis is coincident with \(X_1\) axis. Therefore, the projection of the vector \(\vec{r}_1\) along \(X_1\) axis and the projection of vector \(\vec{r}\) along \(X'\) axes is the same. Hence, \(x\) must be equal to \(x_1\).

Figure 1.5: Rotation from \(X_1Y_1Z_1\) to \(X'Y'Z' (BCS)\)

Figure 1.5 shows the rotation from \(X_1Y_1Z_1\) frame to \(X'Y'Z'\) frame (BCS). Note that the \(X_1\) and \(X'\) axes are coincident and are pointed into the plane. The positive sign convention for \(\phi\) follows the right hand rule about the \(X_1\) (or \(X'\)) axis and is thus clockwise in Figure 1.5. From Figure 1.5 it can be seen that \(y_1\) and \(z_1\) can be expressed in terms of \(y\), \(z\) and \(\phi\) as shown in \(\eqref{eq-x1-rotm}\) and \(\eqref{eq-x1-rotm-expanded}\).

\[\begin{align} \vec{r}_1 = \boldsymbol{R}_{x,\phi} \vec{r} \label{eq-x1-rotm} \end{align}\]

\[\begin{align} \begin{bmatrix} x_1 \\ y_1 \\ z_1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\phi) & -\sin(\phi) \\ 0 & \sin(\phi) & \cos(\phi) \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} \label{eq-x1-rotm-expanded} \end{align}\]

By combining \(\eqref{eq-z3-rotm}\), \(\eqref{eq-y2-rotm}\) and \(\eqref{eq-x1-rotm}\) the vector \(\vec{r}_3\) can be expressed in terms of the vector \(\vec{r}\) as shown in \(\eqref{eq-rotm-fomulation}\).

\[\begin{align} \vec{r}_3 = \boldsymbol{R}_{z,\psi} \vec{r}_2 = \boldsymbol{R}_{z,\psi} \boldsymbol{R}_{y,\theta} \vec{r}_1 = \boldsymbol{R}_{z,\psi} \boldsymbol{R}_{y,\theta} \boldsymbol{R}_{x,\phi} \vec{r} = \boldsymbol{R} \vec{r} \label{eq-rotm-fomulation} \end{align}\]

The combined rotation matrix \(\boldsymbol{R}\) is given by \(\eqref{eq-rotm-expansion}\)

\[\begin{align} \boldsymbol{R} = \begin{bmatrix} c_2c_3 & -c_1s_3 + s_1s_2c_3 & s_1s_3 + c_1s_2c_3 \\ c_2s_3 & c_1c_3 + s_1s_2s_3 & -s_1c_3 + c_1s_2s_3 \\ -s_2 & s_1c_2 & c_1c_2 \end{bmatrix} \label{eq-rotm-expansion} \end{align}\]

where \(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)\). The rotation matrix \(\boldsymbol{R}\) transforms any vector in the BCS to the GCS. The rotation matrix is defined by the Euler angles \(\phi\), \(\theta\) and \(\psi\). Rotation matrices are orthogonal matrices, which means that the transpose of the matrix is its inverse. Mathematically it can be expressed as shown in \(\eqref{eq-rotm-orthogonal}\).

\[\begin{align} \boldsymbol{R}^{-1} = \boldsymbol{R}^T \implies \boldsymbol{R} \boldsymbol{R}^T = \boldsymbol{R}^T\boldsymbol{R} = \boldsymbol{I} \label{eq-rotm-orthogonal} \end{align}\]

1.3.1 Gimbal Lock

It can be seen that when \(\theta = \frac{\pi}{2}\) the rotation matrix \(\boldsymbol{R}\) reduces to as shown in \(\eqref{eq-rotm-gimbal-lock}\).

\[\begin{align} \boldsymbol{R} = \begin{bmatrix} 0 & -c_1s_3 + s_1c_3 & s_1s_3 + c_1c_3 \\ 0 & c_1c_3 + s_1s_3 & -s_1c_3 + c_1s_3 \\ -1 & 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & \sin(\phi - \psi) & \cos(\phi - \psi) \\ 0 & \cos(\phi - \psi) & -\sin(\phi - \psi) \\ -1 & 0 & 0 \end{bmatrix} \label{eq-rotm-gimbal-lock} \end{align}\]

In this case the rotation matrix becomes degenerate. This means that unique mapping is lost at this point. Several combinations of \(\phi\) and \(\psi\) that have the same value of \(\phi - \psi\) will lead to the same rotation matrix. For example, \((\phi, \theta, \psi) = (10^\circ, 90^\circ, 5^\circ)\) and \((\phi, \theta, \psi) = (5^\circ, 90^\circ, 0^\circ)\) will lead to the same rotation matrix. This phenomenon is called gimbal lock. It is called a gimbal lock as in this representation there is no way to roll the vehicle (for the ZYX rotation order) without pitching it out of the gimbal lock position.

For Euler angle formulation based on ZYX order of rotation, gimbal lock will happen when \(\theta = \frac{n\pi}{2}\) for \(n=\pm1, \pm2, ...\) Note that this phenomenon happens for all Euler angle formulations regardless of the order of rotation. Whenever the second rotation results in the first and the third axes of rotations being coincident, this degeneracy will manifest.

Common Confusion

The common confusion is - does physically the rigid body experience any problem when passing through this degenerate point of \(\theta = \pm \frac{\pi}{2}\)?

Remember that the Euler angle notation is our way of representing the pose of the rigid body. The degenracy is in our way of representing the orientation and not with the physical world. So when physically the rigid body moves through a gimbal lock orientation, our representation of orientation will result in a funny looking motion. Quaternions based representation of oreintation does not suffer from the problem of gimbal lock and hence has been widely adopted in the fields of robotics, computer vision, computer graphics, flight dynamics etc.

Check out this YouTube video on Gimbal lock for better visualization.

There exist 12 different possible rotation orders when defining Euler angles. These can be separated into two groups:

  • Proper Euler Angles (ZXZ, XYX, YZY, ZYZ, XZX and YXY)
  • Tait-Bryan Angles (XYZ, ZXY, YZX, XZY, ZYX and YZX)

The Tait-Bryan angles represent rotations about three distinct axes while proper Euler angles use the same axes for the first and third elemental rotations. For marine vehicles, the common rotation order used is ZYX order. Notice that for this choice of rotation order, the gimbal lock orientations occur when the vehicle pitches up or down by \(90^{\circ}\). As marine vehicles like ships and underwater vehicles do not pitch by large angles during regular operations, this representation serves well to describe orientations uniquely. Unless specified explicitly, this text will assume a ZYX order of rotation when referring to Euler angles.

1.4 Kinematic Equations

Let the position vector from the origin of GCS to the origin of the BCS expressed in GCS be defined as shown in \(\eqref{eq-position}\).

\[\begin{align} \vec{r}_0 = \{\eta_1\} = [x_0 \quad y_0 \quad z_0]^T \label{eq-position} \end{align}\]

Let the Euler angles be defined as a vector as shown in \(\eqref{eq-euler-angles}\). Note that this vector does not physically represent any geometric quantity. Remember that each of the angles sequentially denote the rotations from GCS through the intermediate frames to the BCS.

\[\begin{align} \{\eta_2\} = [\phi \quad \theta \quad \psi]^T \label{eq-euler-angles} \end{align}\]

The total pose vector \(\{\eta\}\) is defined as the concatenation of \(\{\eta_1\}\) and \(\{\eta_2\}\) as shown in \(\eqref{eq-pose-vec}\).

\[\begin{align} \{\eta\} = \begin{bmatrix}\{\eta_1\}^T & \{\eta_2\}^T \end{bmatrix}^T = \begin{bmatrix} x_0 & y_0 & z_0 & \phi & \theta & \psi\end{bmatrix}^T \label{eq-pose-vec} \end{align}\]

Let the translational velocity of the vehicle expressed in BCS be defined as shown in \(\eqref{eq-bcs-vel}\). Note that the vector notation and matrix notations are used interchangably. This interchange of notations will be useful later when we would want to express cross product of vectors in terms of matrix muliplications.

\[\begin{align} \vec{v} = \{v_1\} = \begin{bmatrix} u & v & w\end{bmatrix}^T = u \hat{i} + v\hat{j} + w\hat{k} \label{eq-bcs-vel} \end{align}\]

Similarly the angular velocity of the vehicle in BCS is defined as shown in \(\eqref{eq-bcs-ang-vel}\).

\[\begin{align} \vec{\omega} = \{v_2\} = \begin{bmatrix} p & q & r\end{bmatrix}^T = p \hat{i} + q\hat{j} + r\hat{k} \label{eq-bcs-ang-vel} \end{align}\]

The total velocity vector \(\{v\}\) is defined as the concatenation of \(\{v_1\}\) and \(\{v_2\}\) as shown in \(\eqref{eq-vel-vec}\).

\[\begin{align} \{v\} = \begin{bmatrix}\{v_1\}^T & \{v_2\}^T \end{bmatrix}^T = \begin{bmatrix} u & v & w & p & q & r \end{bmatrix}^T \label{eq-vel-vec} \end{align}\]

Now \(\frac{d}{dt}\{\eta\} = \{\dot{\eta}_1\}\) represents the time derivative of the position vector of BCS origin expressed in the GCS. Since GCS is an inertial frame of reference, \(\{\dot{\eta}_1\}\) is the translational velocity of the BCS origin with respect to GCS. The translational velocity of BCS origin expressed in BCS is given by \(\{v_1\}\). These two vectors are related through \(\eqref{eq-kinematics-translation}\).

\[\begin{align} \{\dot{\eta}_1\} = \boldsymbol{R} \{v_1\} \label{eq-kinematics-translation} \end{align}\]

Common Confusion

Why \(\{v_1\} \neq 0\)? \(\{v_1\}\) represents the velocity of the body in BCS frame and if the BCS is moving with the body, then the relative velocity of the body with respect to BCS should be zero?

Note that \(\{v_1\}\) represents the velocity of the BCS orgin with respect to GCS origin expressed in BCS. It is NOT the velocity of the BCS origin with respect to BCS origin expressed in BCS. The velocity of the BCS origin with respect to BCS origin expressed in BCS will definitely be zero!

The angular velocity \(\{v_2\}\) and the derivative of the Euler angles \(\{\dot{\eta}_2\}\) can also be related as was seen for the translational velocity in \(\eqref{eq-kinematics-translation}\). However, the relation is not simply through a rotation matrix!

The Euler angles sequentially denote the rotations from GCS through the intermediate frames to the BCS. Thus each Euler angle is defined with respect to different frames of reference and the contribution of each Euler rate to the total angular velocity needs to be computed separately. The total angular velocity of the body is the sum of angular velocities due to each of the Euler rates.

Let us consider the contribution to angular velocity due to \(\dot{\phi}\). The angle \(\phi\) is defined about the \(X'\) axis (or equivalently \(X_1\) axis). Since \(X'\) axis is one of the cardinal axes of the BCS, the contribution to angular velocity is given by \(\eqref{eq-ang-vel-comp-phi}\).

\[\begin{align} \{\omega_{\dot{\phi}}\} = \begin{bmatrix} \dot{\phi} \\ 0 \\ 0 \end{bmatrix} \label{eq-ang-vel-comp-phi} \end{align}\]

As \(\theta\) is defined about the \(Y_1\) axis (or equivalently \(Y_2\) axis), the contribution to angular velocity due to \(\dot{\theta}\) is given by \(\eqref{eq-ang-vel-comp-theta}\). Note that Euler rate \(\dot{\theta}\) is the angular velocity in the \(X_1Y_1Z_1\) frame and needs to be transformed into BCS. Therefore, it is premultiplied by \(\boldsymbol{R}_{x,\phi}^T\).

\[\begin{align} \{\omega_{\dot{\theta}}\} = \boldsymbol{R}_{x,\phi}^T\begin{bmatrix} 0 \\ \dot{\theta} \\ 0 \end{bmatrix} \label{eq-ang-vel-comp-theta} \end{align}\]

As \(\psi\) is defined about the \(Z_2\) axis (or equivalently \(Z_3\) axis), the contribution to angular velocity due to \(\dot{\psi}\) is given by \(\eqref{eq-ang-vel-comp-psi}\). Note that Euler rate \(\dot{\psi}\) is the angular velocity in the \(X_2Y_2Z_2\) frame and needs to be transformed into BCS. Therefore, it is premultiplied by \(\boldsymbol{R}_{y,\theta}^T\boldsymbol{R}_{x,\phi}^T\).

\[\begin{align} \{\omega_{\dot{\psi}}\} = \boldsymbol{R}_{y,\theta}^T\boldsymbol{R}_{x,\phi}^T\begin{bmatrix} 0 \\ 0 \\ \dot{\psi} \end{bmatrix} \label{eq-ang-vel-comp-psi} \end{align}\]

Thus combining \(\eqref{eq-ang-vel-comp-phi}\), \(\eqref{eq-ang-vel-comp-theta}\) and \(\eqref{eq-ang-vel-comp-psi}\) the total angular velocity is given by \(\eqref{eq-ang-vel-comp}\) and \(\eqref{eq-ang-vel}\).

\[\begin{align} \{\omega\} = \{\omega_{\dot{\phi}}\} + \{\omega_{\dot{\theta}}\} + \{\omega_{\dot{\psi}}\} \label{eq-ang-vel-comp} \end{align}\]

\[\begin{align} \{v_2\} = \{\omega\} = \begin{bmatrix} \dot{\phi} - s_2 \dot{\psi} \\ c_1 \dot{\theta} + s_1c_2 \dot{\psi} \\ -s_1 \dot{\theta} + c_1c_2 \dot{\psi} \end{bmatrix} = \begin{bmatrix} 1 & 0 & -s_2 \\ 0 & c_1 & s_1c_2 \\ 0 & -s_1 & c_1c_2 \end{bmatrix} \begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix} = \boldsymbol{J}_2^{-1} \{\dot{\eta}_2\} \label{eq-ang-vel} \end{align}\]

Thus the Euler rates are given by \(\eqref{eq-kinematics-rotation}\). Note that \(\boldsymbol{J}_2\) matrix becomes singular when gimbal lock occurs (\(\theta=\frac{n\pi}{2}\) for \(n=\pm 1, \pm 2, ...\)).

\[\begin{align} \{\dot{\eta}_2\} = \begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix} = \begin{bmatrix} 1 & \frac{s_1s_2}{c_2} & \frac{c_1s_2}{c_2} \\ 0 & c_1 & -s_1 \\ 0 & \frac{s_1}{c_2} & \frac{c_1}{c_2} \end{bmatrix} \begin{bmatrix} p \\ q \\ r \end{bmatrix} = \boldsymbol{J}_2 \{v_2\} \label{eq-kinematics-rotation} \end{align}\]

Combining \(\eqref{eq-kinematics-translation}\) and \(\eqref{eq-kinematics-rotation}\) yields \(\eqref{eq-kinematics}\).

\[\begin{align} \{\dot{\eta}\} = \begin{bmatrix} \dot{\eta}_1 \\ \dot{\eta}_2 \end{bmatrix} = \begin{bmatrix} \boldsymbol{R}(\vec{\eta}_2) & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{J}_2(\vec{\eta}_2) \end{bmatrix} \{v\} \label{eq-kinematics} \end{align}\]

In a compact form this can be written as shown in \(\eqref{eq-kinematics-compact}\). This matrix differential equation represents the kinematics of a rigid body.

\[\begin{align} \{\dot{\eta}\} = \boldsymbol{J}(\vec{\eta}_2) \{v\} \label{eq-kinematics-compact} \end{align}\]

Note that so far the discussion has been centered around the motion itself and not on the cause of the motion. In the next chapter the dynamic equation of motion that describe the cause of the motion will be discussed.