15 Tutorial 03 - Inertia Matrix

Question 1

Consider a box barge with length \(40 ~m\), breadth \(10 ~m\) and draft \(5 ~m\) carrying cement in an inland waterway with freshwater. Let its BCS origin be located at the intersection of section \(5 ~m\) aft of the center of gravity, waterline and centerline and its x-axis pointed towards the bow and z-axis pointed towards the keel. Let the metacentric height be \(\overline{GM} = 2 ~m\). Calculate the inertia matrix about the BCS if the inertia matrix about center of gravity is given by:

\[\begin{align} I_G = \begin{bmatrix} 32 & 0 & 0 \\ 0 & 200 & 0 \\ 0 & 0 & 200 \end{bmatrix} \times 10^6 ~kg.m^2 \nonumber \end{align}\]

Question 2

Let a rigid body of mass \(2 ~kg\) have an inertia matrix about its BCS as shown below:

\[\begin{align} I_O = \begin{bmatrix} 20 & 5 & -5 \\ 5 & 25 & 5 \\ -5 & 5 & 30 \end{bmatrix} ~kg.m^2 \nonumber \end{align}\]

If \(r_G = [1, 2, 1]^T ~m\), calculate the moment of inertia about the center of gravity G.

Question 3

For the inertia matrix \(I_O\) specified in question 2, calculate the inertia matrix about the principal axes. Principal axes are the axes about which the inertia matrix is diagonal.

Hint: Determine the eigenvalues \((\lambda_1, \lambda_2, \lambda_3)\) and unit eigenvectors \((\vec{v}_1, \vec{v}_2, \vec{v}_3)\) of the inertia matrix. The matrix of eigenvectors defined as:

\[\begin{align} V = \begin{bmatrix} | & | & | \\ \vec{v}_1 & \vec{v}_2 & \vec{v}_3 \\ | & | & | \end{bmatrix} \nonumber \end{align}\]

is the rotation matrix required to transform from BCS to principal axes. By the definition of the eigenvalue problem we have:

\[\begin{align} I_O V = V \Lambda \nonumber \end{align}\]

where

\[\begin{align} \Lambda = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \end{bmatrix} \nonumber \end{align}\]

is the inertia matrix in the principal axes. The angular momentum in BCS is given by:

\[\begin{align} \vec{L}_O = I_O \vec{\omega} \nonumber \end{align}\]

The corresponding vector when expressed in principal axes is given by:

\[\begin{align} \vec{L}_P = V^T \vec{L}_O = V^T I_O \vec{\omega} = V^T I_O V \vec{\omega}_P \end{align}\]

where \(\vec{\omega}_P\) is the angular velocity expressed in principal axes. It can be seen that the inertia matrix in the principal axes is given by:

\[\begin{align} I_P = V^T I_O V = \Lambda \end{align}\]

Question 4

Calculate the euler angles to transform from the BCS frame to the principal axes frame for the inertia matrix shown in question 2.

Solution 1

\[\begin{align} I_G = \begin{bmatrix} 32 & 0 & 0 \\ 0 & 200 & 0 \\ 0 & 0 & 200 \end{bmatrix} \times 10^6 ~kg.m^2 \nonumber \end{align}\]

Code

import numpy as np

def Smat(v):
    return np.array([
        [0, -v[2], v[1]],
        [v[2], 0, -v[0]],
        [-v[1], v[0], 0]
    ])

L = 40
B = 10
T = 5

xG = 5
GM = 2 

I_G = np.array([
    [32, 0, 0],
    [0, 200, 0],
    [0, 0, 200]
]) * 1e6

I_wp = L*B**2/12
Disp = L*B*T
BM = I_wp/Disp
KB = T/2
KG = KB + BM - GM

m = Disp * 1000

zG = T - KG
rG = np.array([xG, 0, zG])

print("The center of gravity is located at")
print(f"x_G = {xG:.2f} m, y_G = {0.0:.2f} m, z_G = {zG:.2f} m\n")

S = Smat(rG)

I_O = I_G - m * (S @ S)

print("The inertia matrix about O is given by:")
print(f"{I_O[0, 0]:.3e}, {I_O[0, 1]:.3e}, {I_O[0, 2]:.3e}")
print(f"{I_O[1, 0]:.3e}, {I_O[1, 1]:.3e}, {I_O[1, 2]:.3e}")
print(f"{I_O[2, 0]:.3e}, {I_O[2, 1]:.3e}, {I_O[2, 2]:.3e}")

The center of gravity is located at
x_G = 5.00 m, y_G = 0.00 m, z_G = 4.33 m

The inertia matrix about O is given by:
6.956e+07, 0.000e+00, -4.333e+07
0.000e+00, 2.876e+08, 0.000e+00
-4.333e+07, 0.000e+00, 2.500e+08

Solution 2

Let a rigid body of mass \(2 ~kg\) have an inertia matrix about its BCS as shown below:

\[\begin{align} I_O = \begin{bmatrix} 20 & 5 & -5 \\ 5 & 25 & 5 \\ -5 & 5 & 30 \end{bmatrix} ~kg.m^2 \nonumber \end{align}\]

If \(r_G = [1, 2, 1]^T ~m\), calculate the moment of inertia about the center of gravity G.

Code

import numpy as np

def Smat(v):
    return np.array([
        [0, -v[2], v[1]],
        [v[2], 0, -v[0]],
        [-v[1], v[0], 0]
    ])

m = 2
rG = np.array([1, 2, 1])
IO = np.array([
    [20, 5, -5],
    [5, 25, 5],
    [-5, 5, 30]
])

S = Smat(rG)
IG = IO + m * (S @ S)

print("The inertia matrix about G is given by:")
print(f"{IG[0, 0]:.2f}, {IG[0, 1]:.2f}, {IG[0, 2]:.2f}")
print(f"{IG[1, 0]:.2f}, {IG[1, 1]:.2f}, {IG[1, 2]:.2f}")
print(f"{IG[2, 0]:.2f}, {IG[2, 1]:.2f}, {IG[2, 2]:.2f}")

The inertia matrix about G is given by:
10.00, 9.00, -3.00
9.00, 21.00, 9.00
-3.00, 9.00, 20.00

Solution 3

For the inertia matrix \(I_O\) specified in question 2, calculate the inertia matrix about the principal axes. Principal axes are the axes about which the inertia matrix is diagonal.

Code

import numpy as np

IO = np.array([
    [20, 5, -5],
    [5, 25, 5],
    [-5, 5, 30]
])

EO, VO = np.linalg.eig(IO)

EO_mat = np.diag(EO)

print("Inertia matrix in the principal axes:")
print(f"{EO_mat[0, 0]:.2f}, {EO_mat[0, 1]:.2f}, {EO_mat[0, 2]:.2f}")
print(f"{EO_mat[1, 0]:.2f}, {EO_mat[1, 1]:.2f}, {EO_mat[1, 2]:.2f}")
print(f"{EO_mat[2, 0]:.2f}, {EO_mat[2, 1]:.2f}, {EO_mat[2, 2]:.2f}")

Inertia matrix in the principal axes:
13.93, 0.00, 0.00
0.00, 27.70, 0.00
0.00, 0.00, 33.38

Solution 4

Calculate the euler angles to transform from the BCS frame to the principal axes frame for the inertia matrix shown in question 2.

Code

import numpy as np

def rotm_to_eul(R, alternate=False):
    eul = np.zeros(3)
    
    eul[1] = np.arcsin(-R[2, 0])
    eul[2] = np.arctan2(R[1, 0] / np.cos(eul[1]), R[0, 0] / np.cos(eul[1]))    
    eul[0] = np.arctan2(R[2, 1] / np.cos(eul[1]), R[2, 2] / np.cos(eul[1]))
    
    # When theta = 0, the alternate is same as original Euler angles
    if eul[1] != 0:
        eul_alt = np.zeros(3)
    else:
        eul_alt = eul
    
    if eul[1] > 0:
        eul_alt[1] = np.pi - eul[1]
        eul_alt[2] = np.arctan2(R[1, 0] / np.cos(eul_alt[1]), R[0, 0] / np.cos(eul_alt[1]))
        eul_alt[0] = np.arctan2(R[2, 1] / np.cos(eul_alt[1]), R[2, 2] / np.cos(eul_alt[1]))
    else:
        eul_alt[1] = -np.pi - eul[1]
        eul_alt[2] = np.arctan2(R[2, 1] / np.cos(eul_alt[1]), R[2, 2] / np.cos(eul_alt[1]))
        eul_alt[0] = np.arctan2(R[1, 0] / np.cos(eul_alt[1]), R[0, 0] / np.cos(eul_alt[1]))
        
    
    if alternate:
        return eul_alt
    else:
        return eul  

IO = np.array([
    [20, 5, -5],
    [5, 25, 5],
    [-5, 5, 30]
])

EO, VO = np.linalg.eig(IO)

eul = rotm_to_eul(VO)

print("Euler angles to transform from BCS to principal axes are:")
print(f"phi: {np.degrees(eul[0]):.2f} deg, theta: {np.degrees(eul[1]):.2f} deg, psi: {np.degrees(eul[2]):.2f} deg")

Euler angles to transform from BCS to principal axes are:
phi: -14.72 deg, theta: 23.40 deg, psi: 145.44 deg