| name | ship-dynamics-6dof |
| version | 1.0.0 |
| description | 6DOF ship dynamics, equations of motion, seakeeping analysis, and natural frequency calculations |
| author | workspace-hub |
| category | subject-matter-expert |
| tags | ["6dof","ship-dynamics","seakeeping","equations-of-motion","natural-frequency","vessel-motions"] |
| platforms | ["engineering"] |
Ship Dynamics (6DOF) SME Skill
Comprehensive 6 degrees of freedom ship dynamics expertise including equations of motion, seakeeping analysis, natural frequencies, and coupled motion analysis.
When to Use This Skill
Use 6DOF ship dynamics when:
- Equations of motion - Set up and solve 6DOF coupled equations
- Natural frequencies - Calculate natural periods for all DOFs
- Seakeeping analysis - Predict vessel motions in waves
- Coupled dynamics - Analyze roll-heave-pitch coupling
- Time-domain simulation - Integrate equations of motion
- Frequency-domain analysis - RAO-based motion prediction
Core Knowledge Areas
1. 6 Degrees of Freedom
DOF Definition:
translations:
surge:
positive: "Forward"
typical_natural_period: "50-150 seconds"
sway:
positive: "Port"
typical_natural_period: "50-150 seconds"
heave:
positive: "Upward"
typical_natural_period: "6-15 seconds"
rotations:
roll:
positive: "Starboard down"
typical_natural_period: "15-30 seconds"
pitch:
positive: "Bow up"
typical_natural_period: "6-12 seconds"
yaw:
positive: "Bow to starboard"
typical_natural_period: "60-200 seconds"
2. Equations of Motion
General Form:
[M + A(ω)]{ẍ} + [B(ω)]{ẋ} + [C]{x} = {F(t)}
Where:
- [M] = Mass/inertia matrix (6x6)
- [A] = Added mass matrix (6x6, frequency-dependent)
- [B] = Damping matrix (6x6, frequency-dependent)
- [C] = Hydrostatic restoring matrix (6x6)
- {F} = External force vector (6x1)
- {x} = Displacement vector [surge, sway, heave, roll, pitch, yaw]
Mass Matrix:
import numpy as np
def create_mass_matrix(
mass: float,
radii_of_gyration: dict,
center_of_gravity: np.ndarray = None
) -> np.ndarray:
"""
Create 6x6 mass matrix for vessel.
Args:
mass: Vessel mass (tonnes)
radii_of_gyration: {'Rxx': roll, 'Ryy': pitch, 'Rzz': yaw} (m)
center_of_gravity: [x, y, z] from origin (m)
Returns:
6x6 mass matrix
"""
if center_of_gravity is None:
center_of_gravity = np.zeros(3)
xg, yg, zg = center_of_gravity
m = mass * 1000
Ixx = m * radii_of_gyration['Rxx']**2
Iyy = m * radii_of_gyration['Ryy']**2
Izz = m * radii_of_gyration['Rzz']**2
M = np.array([
[m, 0, 0, 0, m*zg, -m*yg],
[0, m, 0, -m*zg, 0, m*xg],
[0, 0, m, m*yg, -m*xg, 0 ],
[0, -m*zg, m*yg, Ixx, 0, 0 ],
[m*zg, 0, -m*xg, 0, Iyy, 0 ],
[-m*yg, m*xg, 0, 0, 0, Izz ]
])
return M
M_fpso = create_mass_matrix(
mass=150000,
radii_of_gyration={
'Rxx': 22,
'Ryy': 95,
'Rzz': 95
},
center_of_gravity=np.array([160, 0, 15])
)
print("Mass Matrix:")
print(M_fpso)
3. Natural Frequencies and Periods
Uncoupled Natural Frequency:
def calculate_natural_frequency_uncoupled(
mass: float,
stiffness: float
) -> dict:
"""
Calculate natural frequency for single DOF.
ω_n = sqrt(K / M)
T_n = 2π / ω_n
Args:
mass: Mass or moment of inertia
stiffness: Stiffness or restoring coefficient
Returns:
Natural frequency and period
"""
omega_n = np.sqrt(stiffness / mass)
period_n = 2 * np.pi / omega_n
frequency_hz = omega_n / (2 * np.pi)
return {
'omega_rad_s': omega_n,
'frequency_hz': frequency_hz,
'period_s': period_n
}
m = 150000 * 1000
A33 = 50000 * 1000
K33 = 1025 * 9.81 * 15000
heave_freq = calculate_natural_frequency_uncoupled(
mass=m + A33,
stiffness=K33
)
print(f"Heave natural period: {heave_freq['period_s']:.2f} seconds")
Coupled Natural Frequencies:
def calculate_coupled_natural_frequencies(
mass_matrix: np.ndarray,
stiffness_matrix: np.ndarray
) -> dict:
"""
Calculate coupled natural frequencies from eigenvalue problem.
det([K] - ω²[M]) = 0
Args:
mass_matrix: 6x6 mass matrix (including added mass)
stiffness_matrix: 6x6 stiffness matrix
Returns:
Natural frequencies for all modes
"""
eigenvalues, eigenvectors = np.linalg.eig(
np.linalg.solve(mass_matrix, stiffness_matrix)
)
omega_n = np.sqrt(eigenvalues.real)
periods = 2 * np.pi / omega_n
sort_idx = np.argsort(periods)
periods = periods[sort_idx]
omega_n = omega_n[sort_idx]
eigenvectors = eigenvectors[:, sort_idx]
dof_names = ['Surge', 'Sway', 'Heave', 'Roll', 'Pitch', 'Yaw']
return {
'periods_s': periods,
'frequencies_rad_s': omega_n,
'frequencies_hz': omega_n / (2*np.pi),
'mode_shapes': eigenvectors,
'dof_names': dof_names
}
M_total = M_fpso + np.diag([15000e3, 15000e3, 50000e3, 1e9, 1e9, 5e8])
K = np.diag([0, 0, 150e6, 5e9, 8e9, 0])
natural_freq = calculate_coupled_natural_frequencies(M_total, K)
print("Natural Periods:")
for i, (dof, T) in enumerate(zip(natural_freq['dof_names'], natural_freq['periods_s'])):
print(f" {dof}: {T:.2f} seconds")
4. Hydrostatic Restoring
Complete Stiffness Matrix:
def calculate_complete_hydrostatic_stiffness(
rho: float,
g: float,
displacement: float,
waterplane_area: float,
waterplane_inertia: dict,
center_of_buoyancy: np.ndarray,
center_of_gravity: np.ndarray,
metacentric_height: dict
) -> np.ndarray:
"""
Calculate complete 6x6 hydrostatic stiffness matrix.
Args:
rho: Water density (kg/m³)
g: Gravity (m/s²)
displacement: Volume displacement (m³)
waterplane_area: Waterplane area (m²)
waterplane_inertia: {'Ixx': Ixx, 'Iyy': Iyy} second moments (m⁴)
center_of_buoyancy: [xb, yb, zb] (m)
center_of_gravity: [xg, yg, zg] (m)
metacentric_height: {'GMT': transverse, 'GML': longitudinal} (m)
Returns:
6x6 hydrostatic stiffness matrix
"""
xb, yb, zb = center_of_buoyancy
xg, yg, zg = center_of_gravity
C = np.zeros((6, 6))
C[2, 2] = rho * g * waterplane_area
C[3, 3] = rho * g * displacement * metacentric_height['GMT']
C[4, 4] = rho * g * displacement * metacentric_height['GML']
C[2, 4] = -rho * g * waterplane_area * xb
C[4, 2] = C[2, 4]
C[2, 3] = -rho * g * waterplane_area * yb
C[3, 2] = C[2, 3]
C[3, 4] = -rho * g * displacement * (zg - zb)
C[4, 3] = C[3, 4]
return C
C_hydro = calculate_complete_hydrostatic_stiffness(
rho=1025,
g=9.81,
displacement=150000,
waterplane_area=15000,
waterplane_inertia={'Ixx': 5e5, 'Iyy': 3e7},
center_of_buoyancy=np.array([160, 0, -10]),
center_of_gravity=np.array([160, 0, 15]),
metacentric_height={'GMT': 3.0, 'GML': 5.0}
)
print("Hydrostatic Stiffness Matrix (diagonal terms):")
print(np.diag(C_hydro))
5. Time-Domain Simulation
Newmark-Beta Integration:
def newmark_beta_integration(
M: np.ndarray,
C: np.ndarray,
K: np.ndarray,
F_t: np.ndarray,
x0: np.ndarray,
v0: np.ndarray,
t: np.ndarray,
beta: float = 0.25,
gamma: float = 0.5
) -> dict:
"""
Newmark-Beta time integration for 6DOF dynamics.
Args:
M: Mass matrix (6x6)
C: Damping matrix (6x6)
K: Stiffness matrix (6x6)
F_t: Force time series (n_steps x 6)
x0: Initial displacement (6,)
v0: Initial velocity (6,)
t: Time array
beta: Newmark beta parameter (0.25 = const accel)
gamma: Newmark gamma parameter (0.5)
Returns:
Dictionary with motion time series
"""
n_steps = len(t)
dt = t[1] - t[0]
x = np.zeros((n_steps, 6))
v = np.zeros((n_steps, 6))
a = np.zeros((n_steps, 6))
x[0] = x0
v[0] = v0
a[0] = np.linalg.solve(M, F_t[0] - C @ v[0] - K @ x[0])
K_eff = K + gamma/(beta*dt) * C + 1/(beta*dt**2) * M
for i in range(n_steps - 1):
F_eff = (
F_t[i+1] +
M @ (x[i]/(beta*dt**2) + v[i]/(beta*dt) + (0.5/beta - 1)*a[i]) +
C @ (gamma/(beta*dt)*x[i] + (gamma/beta - 1)*v[i] + dt*(gamma/(2*beta) - 1)*a[i])
)
x[i+1] = np.linalg.solve(K_eff, F_eff)
a[i+1] = (x[i+1] - x[i])/(beta*dt**2) - v[i]/(beta*dt) - (0.5/beta - 1)*a[i]
v[i+1] = v[i] + dt*((1-gamma)*a[i] + gamma*a[i+1])
return {
'time': t,
'displacement': x,
'velocity': v,
'acceleration': a
}
t = np.linspace(0, 100, 10000)
dt = t[1] - t[0]
M_simple = np.diag([0, 0, 200e6, 0, 0, 0])
C_simple = np.diag([0, 0, 100e6, 0, 0, 0])
K_simple = np.diag([0, 0, 150e6, 0, 0, 0])
F = np.zeros((len(t), 6))
F[:, 2] = 1e6 * np.sin(2*np.pi*t / 10)
result = newmark_beta_integration(
M_simple, C_simple, K_simple, F,
x0=np.zeros(6), v0=np.zeros(6), t=t
)
print(f"Max heave: {np.max(np.abs(result['displacement'][:, 2])):.2f} m")
6. Seakeeping Analysis
Motion Statistics:
def calculate_seakeeping_statistics(
motion_time_series: np.ndarray,
dt: float,
dof_name: str = "Motion"
) -> dict:
"""
Calculate seakeeping statistics from motion time series.
Args:
motion_time_series: Time series of motion
dt: Time step
dof_name: Name of DOF
Returns:
Statistical parameters
"""
mean = np.mean(motion_time_series)
std = np.std(motion_time_series)
sorted_amplitudes = np.sort(np.abs(motion_time_series))
n_third = len(sorted_amplitudes) // 3
significant_amplitude = np.mean(sorted_amplitudes[-n_third:])
max_amplitude = np.max(np.abs(motion_time_series))
zero_crossings = np.where(np.diff(np.sign(motion_time_series)))[0]
if len(zero_crossings) > 1:
Tz = np.mean(np.diff(zero_crossings)) * dt * 2
else:
Tz = np.nan
rms = np.sqrt(np.mean(motion_time_series**2))
return {
'dof': dof_name,
'mean': mean,
'std_dev': std,
'rms': rms,
'significant_amplitude': significant_amplitude,
'max_amplitude': max_amplitude,
'zero_crossing_period': Tz
}
heave_motion = result['displacement'][:, 2]
heave_stats = calculate_seakeeping_statistics(heave_motion, dt, 'Heave')
print(f"Heave Statistics:")
print(f" Significant amplitude: {heave_stats['significant_amplitude']:.2f} m")
print(f" Max amplitude: {heave_stats['max_amplitude']:.2f} m")
print(f" Zero-crossing period: {heave_stats['zero_crossing_period']:.2f} s")
Motion Sickness Incidence (MSI):
def calculate_motion_sickness_incidence(
acceleration_rms: float,
frequency_hz: float,
duration_hours: float = 2
) -> float:
"""
Calculate Motion Sickness Incidence (MSI) using ISO 2631-1.
MSI = % of people experiencing motion sickness
Args:
acceleration_rms: RMS vertical acceleration (m/s²)
frequency_hz: Dominant frequency (Hz)
duration_hours: Exposure duration (hours)
Returns:
MSI percentage
"""
if 0.1 <= frequency_hz <= 0.5:
weighting = 1.0
else:
weighting = 0.5
a_w = acceleration_rms * weighting
time_factor = (duration_hours / 2) ** 0.5
msdv = a_w * time_factor
MSI = 100 * (1 / (1 + np.exp(-(msdv - 3.5) / 0.7)))
return MSI
a_heave = np.std(np.diff(result['velocity'][:, 2]) / dt)
freq_heave = 1 / heave_stats['zero_crossing_period']
msi = calculate_motion_sickness_incidence(a_heave, freq_heave, duration_hours=2)
print(f"Motion Sickness Incidence (2 hours): {msi:.1f}%")
Complete Examples
Example 1: Full 6DOF Simulation
def simulate_vessel_6dof_in_waves(
vessel_properties: dict,
wave_conditions: dict,
duration: float = 3600,
dt: float = 0.1
) -> dict:
"""
Complete 6DOF vessel simulation in irregular waves.
Args:
vessel_properties: Vessel mass, stiffness, damping
wave_conditions: Hs, Tp, heading
duration: Simulation duration (s)
dt: Time step (s)
Returns:
Complete motion results
"""
import plotly.graph_objects as go
from plotly.subplots import make_subplots
t = np.arange(0, duration, dt)
n_steps = len(t)
M = vessel_properties['mass_matrix']
C_damp = vessel_properties['damping_matrix']
K = vessel_properties['stiffness_matrix']
Hs = wave_conditions['Hs']
Tp = wave_conditions['Tp']
heading = wave_conditions['heading']
omega_p = 2 * np.pi / Tp
F_wave = np.zeros((n_steps, 6))
if heading == 0:
F_wave[:, 0] = Hs * 1e5 * np.sin(omega_p * t)
F_wave[:, 2] = Hs * 5e5 * np.sin(omega_p * t)
F_wave[:, 4] = Hs * 1e6 * np.sin(omega_p * t)
elif heading == 90:
F_wave[:, 1] = Hs * 1e5 * np.sin(omega_p * t)
F_wave[:, 3] = Hs * 2e6 * np.sin(omega_p * t)
for i in range(6):
F_wave[:, i] += np.random.randn(n_steps) * 0.2 * np.std(F_wave[:, i])
result = newmark_beta_integration(
M, C_damp, K, F_wave,
x0=np.zeros(6), v0=np.zeros(6), t=t
)
dof_names = ['Surge', 'Sway', 'Heave', 'Roll', 'Pitch', 'Yaw']
statistics = {}
for i, dof in enumerate(dof_names):
statistics[dof] = calculate_seakeeping_statistics(
result['displacement'][:, i], dt, dof
)
fig = make_subplots(
rows=2, cols=3,
subplot_titles=dof_names
)
for i, dof in enumerate(dof_names):
row = i // 3 + 1
col = i % 3 + 1
fig.add_trace(
go.Scatter(
x=t,
y=result['displacement'][:, i],
name=dof,
showlegend=False
),
row=row, col=col
)
fig.update_layout(
title=f'6DOF Vessel Motions (Hs={Hs}m, Tp={Tp}s, Heading={heading}°)',
height=800
)
fig.write_html('reports/6dof_simulation.html')
return {
'time': t,
'motions': result,
'statistics': statistics,
'wave_conditions': wave_conditions
}
vessel = {
'mass_matrix': M_fpso,
'damping_matrix': np.diag([50e6, 50e6, 100e6, 5e8, 5e8, 2e8]),
'stiffness_matrix': C_hydro
}
waves = {
'Hs': 8.5,
'Tp': 12.0,
'heading': 0
}
results = simulate_vessel_6dof_in_waves(vessel, waves, duration=600, dt=0.1)
print("Motion Statistics:")
for dof, stats in results['statistics'].items():
print(f"{dof}: Sig = {stats['significant_amplitude']:.2f}, Max = {stats['max_amplitude']:.2f}")
Example 2: Natural Frequency Sensitivity Study
def natural_frequency_sensitivity_study(
base_properties: dict,
parameter_ranges: dict
) -> dict:
"""
Sensitivity study of natural frequencies to design parameters.
Args:
base_properties: Base vessel properties
parameter_ranges: Parameters to vary
Returns:
Sensitivity results
"""
import plotly.graph_objects as go
results = {}
for param_name, param_values in parameter_ranges.items():
natural_periods = []
for value in param_values:
props = base_properties.copy()
if param_name == 'GMT':
props['K'][3, 3] *= value / base_properties['GMT']
elif param_name == 'Rxx':
m = props['M'][0, 0]
props['M'][3, 3] = m * value**2
freq_result = calculate_coupled_natural_frequencies(
props['M'], props['K']
)
natural_periods.append(freq_result['periods_s'][3])
results[param_name] = {
'values': param_values,
'roll_periods': natural_periods
}
fig = go.Figure()
for param_name, data in results.items():
fig.add_trace(go.Scatter(
x=data['values'],
y=data['roll_periods'],
name=param_name,
mode='lines+markers'
))
fig.update_layout(
title='Roll Natural Period Sensitivity',
xaxis_title='Parameter Value',
yaxis_title='Roll Natural Period (s)'
)
fig.write_html('reports/sensitivity_analysis.html')
return results
base = {
'M': M_fpso,
'K': C_hydro,
'GMT': 3.0
}
param_ranges = {
'GMT': np.linspace(1.0, 5.0, 20),
'Rxx': np.linspace(18, 26, 20)
}
sensitivity = natural_frequency_sensitivity_study(base, param_ranges)
Resources
- Principles of Naval Architecture: SNAME
- Ship Hydrostatics and Stability: Adrian Biran
- Seakeeping: Ship Behaviour in Rough Weather: A.R.J.M. Lloyd
- DNV-RP-C205: Environmental Conditions and Environmental Loads
- ISO 2631-1: Mechanical vibration and shock
Use this skill for all 6DOF dynamics analysis in DigitalModel!
