# Copyright (C) 2025 Antoine ANCEAU
#
# This file is part of pynavaltoolbox.
#
# pynavaltoolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Hydrostatic calculations module.
This module provides the HydrostaticsCalculator class for computing
hydrostatic properties of hull geometries.
"""
import vtk
import numpy as np
from typing import Tuple, Optional, TYPE_CHECKING
from scipy.optimize import minimize
if TYPE_CHECKING:
from ..vessel import Vessel
from .dataclasses import HydrostaticState
from .tools import calculate_submerged_properties, calculate_shell_volume, calculate_section_area
[docs]
class HydrostaticsCalculator:
"""
Calculates hydrostatic properties for vessel geometries.
Supports both single-hull and multi-hull vessels. For multi-hull vessels,
hydrostatic properties are calculated for each hull individually and then
combined (by summation for volumes/areas, or weighted averages for centers).
It also correctly handles hulls that are fully submerged (e.g., submarine bodies)
or fully emerged (e.g., dry dock scenarios).
Attributes:
vessel: The Vessel object to analyze.
water_density: Density of water in kg/m³.
shell_thickness: Shell thickness in meters (adds to submerged volume).
"""
[docs]
def __init__(self, vessel: 'Vessel', water_density: float = 1025.0, shell_thickness: float = 0.0):
"""
Initialize the hydrostatics calculator.
Args:
vessel: The Vessel object.
water_density: Density of water in kg/m³ (default: 1025 for seawater).
shell_thickness: Shell thickness in meters (default: 0 for no shell).
Adds volume to account for water displaced by shell thickness.
"""
self.vessel = vessel
self.water_density = water_density
self.shell_thickness = shell_thickness
[docs]
def calculate_at_draft(self,
draft: float,
trim: float = 0,
heel: float = 0,
vcg: float = 0) -> HydrostaticState:
"""
Calculates hydrostatics for a fixed draft, trim, and heel.
Args:
draft: Draft at the reference point in meters.
trim: Trim angle in degrees (positive = bow down).
heel: Heel angle in degrees (positive = starboard down).
vcg: Vertical Center of Gravity (KG) in meters.
Returns:
HydrostaticState object with all calculated properties.
"""
# MULTI-HULL SUPPORT: Process each hull individually then combine
from ..hull import Hull
# Determine pivot point for transformations (vessel center)
vessel_bounds = self.vessel.get_bounds()
center_x = (vessel_bounds[0] + vessel_bounds[1]) / 2
center_y = (vessel_bounds[2] + vessel_bounds[3]) / 2
# Lists to accumulate properties from each hull
hull_volumes = []
hull_cb_coords = []
hull_wetted_areas = []
hull_waterplane_areas = []
hull_cf_coords = []
hull_i_t = []
hull_i_l = []
hull_midship_areas = []
all_sub_bounds = []
# Process each hull
for hull_geometry in self.vessel.hulls:
# Create a working copy
working_hull = Hull.__new__(Hull)
working_hull.polydata = vtk.vtkPolyData()
working_hull.polydata.DeepCopy(hull_geometry.polydata)
working_hull.file_path = hull_geometry.file_path
# Apply trim and heel rotations using vessel center as pivot
if trim != 0 or heel != 0:
working_hull.transform(
rotation=(heel, trim, 0), # Roll (Heel), Pitch (Trim), Yaw
pivot=(center_x, center_y, draft)
)
# Calculate submerged properties for this hull
submerged_volume, cb_coords, wetted_area, sub_bounds = calculate_submerged_properties(
working_hull.polydata,
surface_function=lambda x, y: draft,
surface_type='planar',
draft=draft
)
# Add shell volume if thickness is specified
if self.shell_thickness > 0:
shell_volume = calculate_shell_volume(wetted_area, self.shell_thickness)
submerged_volume += shell_volume
# Calculate waterplane properties for this hull
waterplane_area, cf_coords, i_t, i_l, _, _ = self._calculate_waterplane_properties(working_hull, draft)
# Calculate midship area for this hull
hull_bounds = working_hull.get_bounds()
midship_x = (hull_bounds[0] + hull_bounds[1]) / 2.0
midship_area = calculate_section_area(working_hull.polydata, midship_x, draft)
# Store properties
hull_volumes.append(submerged_volume)
hull_cb_coords.append(cb_coords)
hull_wetted_areas.append(wetted_area)
hull_waterplane_areas.append(waterplane_area)
hull_cf_coords.append(cf_coords)
hull_i_t.append(i_t)
hull_i_l.append(i_l)
hull_midship_areas.append(midship_area)
all_sub_bounds.append(sub_bounds)
# COMBINE PROPERTIES FROM ALL HULLS
# Total volumes and areas (sum)
total_volume = sum(hull_volumes)
total_wetted_area = sum(hull_wetted_areas)
total_waterplane_area = sum(hull_waterplane_areas)
total_midship_area = sum(hull_midship_areas)
# Combined center of buoyancy (volume-weighted average)
if total_volume > 0:
combined_lcb = sum(vol * cb[0] for vol, cb in zip(hull_volumes, hull_cb_coords)) / total_volume
combined_tcb = sum(vol * cb[1] for vol, cb in zip(hull_volumes, hull_cb_coords)) / total_volume
combined_vcb = sum(vol * cb[2] for vol, cb in zip(hull_volumes, hull_cb_coords)) / total_volume
else:
combined_lcb = combined_tcb = combined_vcb = 0.0
# Combined center of flotation (area-weighted average)
if total_waterplane_area > 0:
combined_lcf = sum(area * cf[0] for area, cf in zip(hull_waterplane_areas, hull_cf_coords)) / total_waterplane_area
else:
combined_lcf = 0.0
# Combined moments of inertia using parallel axis theorem
# I_total = sum(I_i + A_i * d_i^2) where d_i is distance from centroid to combined centroid
combined_i_t = 0.0
combined_i_l = 0.0
for i, (area, cf, i_t_local, i_l_local) in enumerate(zip(hull_waterplane_areas, hull_cf_coords, hull_i_t, hull_i_l)):
# Distance from this hull's centroid to combined centroid
dy = cf[1] - combined_tcb # Transverse distance
dx = cf[0] - combined_lcb # Longitudinal distance
# Parallel axis theorem
combined_i_t += i_t_local + area * dy**2
combined_i_l += i_l_local + area * dx**2
# Calculate combined waterline dimensions
# LWL and BWL from combined bounds
if all_sub_bounds:
combined_xmin = min(b[0] for b in all_sub_bounds)
combined_xmax = max(b[1] for b in all_sub_bounds)
combined_ymin = min(b[2] for b in all_sub_bounds)
combined_ymax = max(b[3] for b in all_sub_bounds)
lwl = combined_xmax - combined_xmin
bwl = combined_ymax - combined_ymin
los = lwl # Length overall submerged
else:
lwl = bwl = los = 0.0
# Calculate displacement mass
displacement_mass = total_volume * self.water_density
# Calculate AP/FP drafts (transform perpendicular points)
draft_ap = draft
draft_fp = draft
if trim != 0 or heel != 0:
t = vtk.vtkTransform()
t.Translate(center_x, center_y, draft)
t.RotateZ(0)
t.RotateY(trim)
t.RotateX(heel)
t.Translate(-center_x, -center_y, -draft)
ap_x = self.vessel.ap
fp_x = self.vessel.fp
p_ap = t.TransformDoublePoint(ap_x, 0.0, 0.0)
p_fp = t.TransformDoublePoint(fp_x, 0.0, 0.0)
draft_ap = draft - p_ap[2]
draft_fp = draft - p_fp[2]
# Calculate coefficients
cm = total_midship_area / (bwl * draft) if bwl > 0 and draft > 0 else 0.0
cb = total_volume / (lwl * bwl * draft) if lwl > 0 and bwl > 0 and draft > 0 else 0.0
cp = total_volume / (total_midship_area * lwl) if total_midship_area > 0 and lwl > 0 else 0.0
# Calculate metacentric radii
bmt = combined_i_t / total_volume if total_volume > 0 else 0.0
bml = combined_i_l / total_volume if total_volume > 0 else 0.0
kmt = combined_vcb + bmt
kml = combined_vcb + bml
# Calculate GM (dry) - without free surface correction
# GM = KM - KG = VCB + BM - VCG
gmt_dry = kmt - vcg
gml_dry = kml - vcg
# Calculate free surface corrections from tanks
# GG' = FSC / Δ where FSC is in m⁴ and Δ is displacement in kg
free_surface_correction_t = 0.0
free_surface_correction_l = 0.0
if displacement_mass > 0 and hasattr(self.vessel, 'tanks'):
fsc_t, fsc_l = self.vessel.get_total_free_surface_correction()
# FSC is in m⁴, need to convert to GG' in meters
# GG' = (FSC * ρ_water) / Δ = FSC * ρ_water / (Volume * ρ_water) = FSC / Volume
# But FSC already includes density ratio, so:
# GG' = FSC / Volume (since FSC = I * ρ_fluid / ρ_water)
if total_volume > 0:
free_surface_correction_t = fsc_t / total_volume
free_surface_correction_l = fsc_l / total_volume
# Calculate GM (wet) - with free surface correction
gmt_wet = gmt_dry - free_surface_correction_t
gml_wet = gml_dry - free_surface_correction_l
# Calculate stiffness matrix (using corrected GM)
stiffness = self._calculate_stiffness_matrix(
displacement_mass, (combined_lcb, combined_tcb, combined_vcb),
vcg, total_waterplane_area, bmt, bml,
free_surface_correction_t, free_surface_correction_l
)
# Create and return the state
return HydrostaticState(
draft=draft,
trim=trim,
heel=heel,
displacement_volume=total_volume,
displacement_mass=displacement_mass,
lcb=combined_lcb,
tcb=combined_tcb,
vcb=combined_vcb,
lcg=combined_lcb, # Assuming LCG = LCB at equilibrium
tcg=combined_tcb, # Assuming TCG = TCB at equilibrium
vcg=vcg,
waterplane_area=total_waterplane_area,
wetted_surface_area=total_wetted_area,
lcf=combined_lcf,
kmt=kmt,
kml=kml,
bmt=bmt,
bml=bml,
gmt_dry=gmt_dry,
gml_dry=gml_dry,
gmt_wet=gmt_wet,
gml_wet=gml_wet,
free_surface_correction_t=free_surface_correction_t,
free_surface_correction_l=free_surface_correction_l,
lwl=lwl,
bwl=bwl,
midship_area=total_midship_area,
cm=cm,
stiffness_matrix=stiffness,
los=los,
cb=cb,
cp=cp,
draft_ap=draft_ap,
draft_fp=draft_fp
)
[docs]
def find_equilibrium(self,
displacement_mass: float,
center_of_gravity: Tuple[float, float, float],
initial_draft: Optional[float] = None) -> HydrostaticState:
"""
Finds the equilibrium position for a given displacement and CoG.
Args:
displacement_mass: Target displacement mass in kg.
center_of_gravity: (LCG, TCG, VCG) in meters.
initial_draft: Initial guess for draft. If None, estimated from bounds.
Returns:
HydrostaticState at equilibrium.
"""
lcg, tcg, vcg = center_of_gravity
# Initial guess
if initial_draft is None:
bounds = self.vessel.get_bounds()
initial_draft = (bounds[4] + bounds[5]) / 2
# Better initial guess for trim?
# If LCB is far from LCG, we can estimate trim change needed using BML
# But we don't know BML yet. Start with 0.
x0 = np.array([initial_draft, 0.0, 0.0]) # [draft, trim, heel]
def objective(x):
"""Objective function to minimize."""
draft, trim, heel = x
try:
state = self.calculate_at_draft(draft, trim, heel, vcg)
# Residuals
# Displacement error
disp_error = (state.displacement_mass - displacement_mass) / displacement_mass
# Center of Buoyancy - Center of Gravity error (Moment equilibrium)
# For free floating body, B and G must align in X and Y (projected)
# Actually, this is simplified.
# Equilibrium condition:
# 1. Displacement = Weight
# 2. G matches B in x and y (if heel/trim account for shift)
# But state.lcb is in global coords. G is in global coords (rotated?)
# Wait, G is fixed on the body. We rotate the body.
# calculate_at_draft rotates the hull by (heel, trim).
# So the hull points change.
# G should also be rotated if we are checking vs global B.
# OR, calculate_at_draft returns B in global frame?
# working_hull.transform rotates the mesh. So B is in global frame.
# We need to rotate G into global frame to compare.
# BUT, `calculate_at_draft` input `vcg` is just used for stiffness matrix.
# It does NOT rotate G.
# We need to transform G inside objective function or compare properly.
# HOWEVER, standard approach:
# Equilibrium means Sum of Forces = 0 (Disp = Weight)
# Sum of Moments = 0
# We can just minimize distance between B and G in horizontal plane
# IF we rotate G correctly.
# Let's pivot G.
# Pivot used in calculate_at_draft was center of bounding box at draft.
# We don't have access to that dynamic pivot here easily without recalculating bounds.
# Using 0,0,0 as pivot for simple rotation? No, hull rotates around specific pivot.
# Easier approach used in naval arch:
# Longitudinal moment: Displacement * (LCB - LCG) = 0
# Transverse moment: Displacement * (TCB - TCG) = 0
# calculated in BODY frame or GLOBAL frame?
# Usually we align LCB to LCG in the frame where gravity is vertical.
# Our simulation rotates the HULL so gravity is effectively (0,0,-1).
# So we simply need B_x_global == G_x_global and B_y_global == G_y_global.
# We need to rotate input G (Body frame) to G' (Global frame)
# using the same rotation as hull.
# Reconstruct pivot used in calculate_at_draft:
# This is tricky because calculate_at_draft creates a new hull and gets bounds.
# We can't replicate that exactly without cost.
# Hack: Assume small angles? No, we need robust solving.
# Better: calculate_at_draft should probably accept G and return Righting Arms (GZ, etc).
# But we are finding trim/draft for UPRIGHT (or specified state).
# If we assume we are solving for FREE TRIM/DRAFT at Specific Heel (or free heel):
# We can minimize (LCB - LCG_rotated) and (TCB - TCG_rotated).
# dx, dy = state.lcb - lcg, state.tcb - tcg (Unused here)
# For now, let's stick to the existing simple diff, assuming small angles or
# that checking LCB=LCG is sufficient for trim.
# For trim optimization, rotating LCG is critical if finding large trim.
# LCG_x_global approx LCG - VCG*sin(trim) ?
lcb_error = (state.lcb - lcg)
tcb_error = (state.tcb - tcg)
# Weight errors
return (disp_error * 10.0)**2 + (lcb_error)**2 + (tcb_error)**2
except Exception:
return 1e10
# Optimize
# Use Nelder-Mead for robustness with non-smooth gradients
result = minimize(objective, x0, method='Nelder-Mead',
options={'maxiter': 200, 'xatol': 1e-4})
draft, trim, heel = result.x
return self.calculate_at_draft(draft, trim, heel, vcg)
def _calculate_waterplane_properties(self,
hull,
draft: float) -> Tuple[float, np.ndarray, float, float, float, float]:
"""
Calculate waterplane properties using exact polygon geometry.
Returns:
Area, Centroid [x,y,z], I_t, I_l, lwl, bwl
"""
# Create a plane at the waterline
plane = vtk.vtkPlane()
plane.SetOrigin(0, 0, draft)
plane.SetNormal(0, 0, 1)
# Cut the hull at the waterline
cutter = vtk.vtkCutter()
cutter.SetInputData(hull.polydata)
cutter.SetCutFunction(plane)
cutter.Update()
waterline = cutter.GetOutput()
# Use the tools function
from .tools import calculate_waterplane_geometric_properties
return calculate_waterplane_geometric_properties(waterline)
def _calculate_metacentric_radii(self, hull, draft: float, volume: float) -> Tuple[float, float]:
"""Deprecated: Internal calculation now handled in calculate_at_draft via waterplane props."""
pass
def _calculate_stiffness_matrix(self,
displacement: float,
cb: np.ndarray,
vcg: float,
waterplane_area: float,
bmt: float,
bml: float,
fsc_t: float = 0.0,
fsc_l: float = 0.0) -> np.ndarray:
"""
Calculate the 6x6 hydrostatic stiffness matrix.
Parameters
----------
displacement : float
Displacement mass in kg.
cb : tuple
Center of buoyancy (x, y, z) in meters.
vcg : float
Vertical center of gravity (KG) in meters.
waterplane_area : float
Waterplane area in m².
bmt : float
Transverse metacentric radius (BM_t) in meters.
bml : float
Longitudinal metacentric radius (BM_l) in meters.
fsc_t : float, optional
Transverse free surface correction (GG') in meters.
fsc_l : float, optional
Longitudinal free surface correction (GG') in meters.
Returns
-------
np.ndarray
6x6 stiffness matrix.
"""
K = np.zeros((6, 6))
g = 9.81 # gravity
# Heave stiffness (3,3)
K[2, 2] = self.water_density * g * waterplane_area
# Roll stiffness (4,4) - using corrected GM
if vcg > 0:
gmt_dry = cb[2] + bmt - vcg
gmt_wet = gmt_dry - fsc_t
K[3, 3] = displacement * g * gmt_wet
# Pitch stiffness (5,5) - using corrected GM
if vcg > 0:
gml_dry = cb[2] + bml - vcg
gml_wet = gml_dry - fsc_l
K[4, 4] = displacement * g * gml_wet
# Cross-coupling terms (simplified)
# K[2,4] = K[4,2] = -rho * g * waterplane_area * cb[0]
# K[2,3] = K[3,2] = -rho * g * waterplane_area * cb[1]
return K