This post shows how simulations with State Space Models is used in my research. This is not "rocket science", since it is not very complicated, but actually it is, since these methods were developed during the Apollo project! :)

Roll damping

A systems ability to damp reconance motions can be crucial as in the case of the Tacoma-narrows-bridge-collapse:

Resonance phenomena can also occur for ships, where for instance something called parametric roll can happen as can be seen in this video:

Having enough roll damping is one way to deal with this problem.

Ship roll-decay test

One way to determine the roll damping is to conduct a so called Roll-decay test. This test is conducted by giving the ship an initial roll angle and then observing the resulting oscillation. The oscillation motion will decrease with time due to energy losses from friction with the water and other things such as wave and eddy generation. Here is an example of a roll decay model test:

from scipy.integrate import solve_ivp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from ipywidgets import interactive
from IPython.core.display import HTML

Simulation

Drop simulation

Newtons 2nd law can be used to predict the resulting motion from forces: $$F=m \cdot a$$ so we want to determine acceleration $a$ which can then be integrated to get the velocity and position to describe the motion.

def calculate_acceleration(F,m):
    a = F/m
    return a

Lets simulate a 1 kg object falling without any drag. It can be done by:

  1. calculate acceleration
  2. calculate the velocity by adding acceleration
  3. calculate the position by adding the velocity
dt = 1/10
ts = np.arange(0,10,dt)

g=9.81
F=-g
m=1
velocity=0.0  # initial velocity
position=0.0  # initial position
positions=[]
for n in range(len(ts)):
    
    acceleration = calculate_acceleration(F=F, m=m)
    velocity = velocity + acceleration*dt
    position = position + velocity*dt
    positions.append(position)
positions = np.array(positions)
    
fig,ax=plt.subplots()
ax.plot(ts,positions);
ax.set_xlabel('time [s]')
ax.set_ylabel('position [m]');
ax.grid(True)

And here is an interactive animation of the same thing (use the Binder link at the top to get this running):

#collapse
from matplotlib.lines import Line2D
import matplotlib.animation as animation

class SubplotAnimation(animation.TimedAnimation):
    def __init__(self):
        fig,axes = plt.subplots(ncols=2)
        fig.set_size_inches(12,8)
        self.fig = fig
        self.axes=axes
        
        ax = axes[0]
        ax.set_ylim(-500,0)
        ax.set_xlim(0,0)
        ax.set_ylabel('Position [m]')
        ax.grid(True)
        
        ax = axes[1]
        ax.set_ylim(-500,0)
        ax.set_xlim(0,ts[-1])
        ax.set_xlabel('Time [s]')
        ax.set_ylabel('Position [m]')
        ax.grid(True)
        
        self.t = ts
        self.y = positions
        
        self.ball = Line2D([], [], linestyle='', marker='o',color='b', ms=20)
        self.ball_time = Line2D([], [], color='red', linewidth=2)
        
        axes[0].add_line(self.ball)
        axes[1].add_line(self.ball_time)
        
        return animation.TimedAnimation.__init__(self, fig, interval=50, blit=True)

    def _draw_frame(self, framedata):
        i = framedata
        x = 0
        position=positions[i]
        y = position
            
        self.ball.set_data(x, y)
        
        times = np.arange(i, dtype=int)
        self.ball_time.set_data(self.t[times],positions[times])  
        
        self._drawn_artists = [self.ball,self.ball_time]
                
    def new_frame_seq(self):
        return iter(range(self.t.size))

    def _init_draw(self):
        lines = [self.ball,self.ball_time]

        for l in lines:
            l.set_data([], [])

anim_ball = SubplotAnimation()
plt.close()

#collapse
plt.rcParams["animation.html"] = "jshtml"
HTML(anim_ball.to_jshtml())
</input>

Ship roll simulation

So now doing something similar for the ship:

$$ A\cdot\ddot{\phi} + C\cdot\phi = 0 $$ 

dt = 0.1
ts = np.arange(0,30,dt)

g=9.81
A=1
C=1
phi1d=0.0          # initial velocity
phi=np.deg2rad(10) # initial roll angle

phis=[]
for n in range(len(ts)):
    
    phi2d = -C*phi/A
    phi1d = phi1d + phi2d*dt
    phi = phi + phi1d*dt
    phis.append(phi)
phis = np.array(phis)
    
fig,ax=plt.subplots()
ax.plot(ts,np.rad2deg(phis));
ax.set_xlabel('time [s]')
ax.set_ylabel('phi [deg]');
ax.grid(True)


#collapse
from matplotlib.lines import Line2D
import matplotlib.animation as animation
from matplotlib.patches import Rectangle
from matplotlib import transforms

def add_ship(ax):
    beam=1
    height=0.5
    ship = Rectangle(xy=(-beam/2,-height/2),width=1,height=0.5)
    ax.add_patch(ship)
    
    ax.axis('equal')
    xlim=1.3*beam/2*np.array([-1,1])
    ax.set_xlim(xlim)
    ax.set_ylim(xlim)
    
    ax.plot(xlim,[0,0],'b-');
    
    return ship
    
def rotate_ship(ship,ax,angle):
    t2 = transforms.Affine2D().rotate_deg(-np.rad2deg(angle)) + ax.transData
    ship.set_transform(t2)

class RollDecayAnimation(animation.TimedAnimation):
    def __init__(self, ts, phis):
        
        fig,axes = plt.subplots(ncols=2)
        fig.set_size_inches(12,8)
        self.fig = fig
        self.axes=axes
        
        self.t = ts
        self.phis = phis
        
        ax = axes[0]
        self.ship=add_ship(ax=ax)
        
        ax = axes[1]
        ax.set_ylim(np.min(np.rad2deg(self.phis)),np.rad2deg(np.max(self.phis)))
        ax.set_xlim(0,self.t[-1])
        ax.set_xlabel('Time [s]')
        ax.set_ylabel('roll [deg]')
        ax.grid(True)
    
        
        self.phi_time = Line2D([], [], color='red', linewidth=2)
        
        #axes[0].add_line(self.ball)
        axes[1].add_line(self.phi_time)
        
        return animation.TimedAnimation.__init__(self, fig, interval=30, blit=True)

    def _draw_frame(self, framedata):
        i = framedata
        
        phi=self.phis[i]
        times = np.arange(i, dtype=int)
        self.phi_time.set_data(self.t[times],np.rad2deg(self.phis[times]))  
        
        rotate_ship(ship=self.ship, ax=self.axes[0], angle=phi)
        
        self._drawn_artists = [self.phi_time, self.ship]
                
    def new_frame_seq(self):
        return iter(range(self.t.size))

    def _init_draw(self):
        lines = [self.phi_time]

        for l in lines:
            l.set_data([], [])

anim_roll = RollDecayAnimation(ts=ts, phis=phis)
plt.close()

plt.rcParams["animation.html"] = "jshtml"
HTML(anim_roll.to_jshtml())
</input>

This is however not a damped motions as the motion does not decay. In order to get the decay one must also include a damping term $B$: $$ A\cdot\ddot{\phi} + B\cdot\dot{\phi} + C\cdot\phi = 0 $$

B=0.2
phi1d=0.0          # initial velocity
phi=np.deg2rad(10) # initial roll angle

phis=[]
for n in range(len(ts)):
    
    phi2d = (-C*phi - B*phi1d)/A
    phi1d = phi1d + phi2d*dt
    phi = phi + phi1d*dt
    phis.append(phi)
    
phis = np.array(phis)
    
fig,ax=plt.subplots()
ax.plot(ts,np.rad2deg(phis));
ax.set_xlabel('time [s]')
ax.set_ylabel('phi [deg]');
ax.grid(True)


#collapse
anim_rolldecay = RollDecayAnimation(ts=ts, phis=phis)
plt.close()

plt.rcParams["animation.html"] = "jshtml"
HTML(anim_rolldecay.to_jshtml())
</input>

Think of something that could be simulated in your field of work or research?