21  Tutorial 09 - Full scale maneuvering tests

Question 01

Consider a ship of length \(100 ~m\) and design speed of \(10 ~m/s\) with \(T'=2\) and \(K'=-1\). Assuming that the ship starts from the nominal operating point and starts executing a \(20^{\circ{}}-10^{\circ{}}\) zig-zag maneuver. Determine:

  • Overshoot
  • Time of the overshoot

You may assume that the rudder angle can jump instantaneously (no steering gear inertia). Does the ship satisfy the standard maneuverability criteria as per MSC 137(76)?

Note: For finding the root of a function \(f(t)\), you may use the Newton-Raphson method

\[ t_{k+1} = t_k - \frac{f(t_k)}{f'(t_k)} \]

In your case if you choose initial guess as the time constant, you will converge to the root in 3 iterations.

Question 02

Consider a ship which is straight line unstable. Assuming that the yaw rate \(r(t)\) follows the dynamics:

\[ \dot{r'} - r' + 3r'^3 = -\delta(t) \]

Determine the area under the hysteresis loop observed in \(\delta\) vs \(r'\) plot of the spiral maneuver.

Question 01 Solution

Consider a ship of length \(100 ~m\) and design speed of \(10 ~m/s\) with \(T'=2\) and \(K'=-1\). Assuming that the ship starts from the nominal operating point and starts executing a \(20^{\circ{}}-10^{\circ{}}\) zig-zag maneuver. Determine:

  • Overshoot
  • Time of the overshoot

You may assume that the rudder angle can jump instantaneously (no steering gear inertia). Does the ship satisfy the standard maneuverability criteria as per MSC 137(76)?

Note: For finding the root of a function \(f(t)\), you may use the Newton-Raphson method

\[ t_{k+1} = t_k - \frac{f(t_k)}{f'(t_k)} \]

In your case if you choose initial guess as the time constant, you will converge to the root in 3 iterations.

Code 
import numpy as np

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def solve_problem_01():
    L = 100
    U = 10
    
    Tp = 2
    Kp = -1
    
    T = Tp * L / U
    K = Kp * U / L
    
    delta0 = 20 * np.pi / 180
    psi0 = -10 * np.pi / 180
    
    def psi_r(t, t0, T, K, delta, C1, C2):
        r = K*delta - C1*np.exp(-(t-t0)/T)
        psi = K*delta*(t-t0) + T*C1*np.exp(-(t-t0)/T) + C2
        
        return psi, r
    
    t0_1 = 0
    C1_1 = K*delta0
    C2_1 = -T*C1_1
    
    t_flip_1 = T
    for i in range(3):
        psi, r = psi_r(t_flip_1, 0, T, K, delta0, C1_1, C2_1)
        t_flip_1 = t_flip_1 - (psi - psi0)/r
    
    psi_flip_1, r_flip_1 = psi_r(t_flip_1, 0, T, K, delta0, C1_1, C2_1)
    
    C1_2 = -K*delta0 - r_flip_1
    C2_2 = -T*C1_2 + psi_flip_1
    
    psi_flip_2, r_flip_2 = psi_r(t_flip_1, t_flip_1, T, K, -delta0, C1_2, C2_2)
        
    t_peak = t_flip_1 - T*np.log(-K*delta0 / C1_2)
    psi_peak, r_peak = psi_r(t_peak, t_flip_1, T, K, -delta0, C1_2, C2_2)
    
    def nomoto(t, x, delta):
        psi = x[0]
        r = x[1]
        
        xd = np.zeros(2)
        xd[0] = r
        xd[1] = (K*delta - r)/T
        
        return xd
    
    t_uni = np.linspace(0, 30, 1000)
    
    psi1 = np.zeros_like(t_uni)
    psi2 = np.zeros_like(t_uni)
    r1 = np.zeros_like(t_uni)
    r2 = np.zeros_like(t_uni)
    
    for i in range(len(t_uni)):
        psi1[i], r1[i] = psi_r(t_uni[i], 0, T, K, delta0, C1_1, C2_1)
        psi2[i], r2[i] = psi_r(t_uni[i], t_flip_1, T, K, -delta0, C1_2, C2_2)
    
    sol1 = solve_ivp(nomoto, (t_uni[0], t_flip_1), np.zeros(2), args=(delta0,), t_eval=t_uni[t_uni < t_flip_1])
    sol2 = solve_ivp(nomoto, (t_flip_1, t_uni[-1]), sol1.y[:, -1], args=(-delta0,), t_eval = t_uni[t_uni > t_flip_1])
    
    plt.figure()
    plt.plot(t_uni, psi1*180/np.pi, 'r', label="$\\psi_1(t)$")
    plt.plot(t_uni, psi2*180/np.pi, 'b', label="$\\psi_2(t)$")
    plt.plot(sol1.t, sol1.y[0]*180/np.pi, 'g--', label="$\\psi_1(t)$ Numerical")
    plt.plot(sol2.t, sol2.y[0]*180/np.pi, 'g--', label="$\\psi_2(t)$ Numerical")
    plt.plot(t_peak, psi_peak*180/np.pi,'kp')
    plt.xlabel('Time (s)')
    plt.ylabel('$\\psi(t)$')
    plt.legend()
    plt.grid()
    
solve_problem_01()

Question 02 Solution

Consider a ship which is straight line unstable. Assuming that the yaw rate \(r(t)\) follows the dynamics:

\[ \dot{r'} - r' + 3r'^3 = -\delta(t) \]

Determine the area under the hysteresis loop observed in \(\delta\) vs \(r'\) plot of the spiral maneuver.

Code 
import numpy as np

def solve_problem_02():

    def delta(r):
        return r - 3*r**3
        
    r1_slope0 = 1/3
    r2_slope0 = -1/3

    d2 = delta(r2_slope0)
    r_max = np.max(np.real(np.roots([-3, 0, 1, 2/9])))
    r_min = np.min(np.real(np.roots([3, 0, -1, 2/9])))

    p = np.polyint([-3, 0, 1, 2/9])
    
    print(f"Area of hysteresis loop: {2*(np.polyval(p, r_max) - np.polyval(p, r2_slope0)):.2f}")

solve_problem_02()
Area of hysteresis loop: 0.50