feat(stability): Módulo 3 — Curva GZ + criterios IMO IS Code 2008
- gz_integrator.py: GZCurve, GZPoint, compute_gz_wall_sided (fórmula pared lateral), compute_gz_direct (integración Sutherland-Hodgman) - imo_is2008.py: IMOCriterion, IMOResult, check_imo_is2008 — 6 criterios A.2.1.1–A.2.1.6 del IS Code 2008 Cap.2 - gz_curve_widget.py: GZCurveWidget QPainter — curva cian, áreas sombreadas, líneas IMO, marcador AVS, tabla PASS/FAIL integrada - main_window.py: GZCurveWidget en MOD_STABILITY, _compute_and_show_gz, _on_show_stability conectado al ribbon - dark.qss: estilos GZCurveWidget - test_module3_stability.py: 33 tests S-01..S-28 (315 total, todos pasan) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
@@ -1 +1,31 @@
|
||||
# arshipdesign/stability
|
||||
"""
|
||||
arshipdesign.stability — estabilidad estática e intacta.
|
||||
|
||||
Exporta:
|
||||
- GZPoint, GZCurve — estructuras de datos
|
||||
- compute_gz_wall_sided — fórmula de pared lateral (rápida)
|
||||
- compute_gz_direct — integración directa 3-D (precisa)
|
||||
- IMOCriterion, IMOResult — criterios IMO IS Code 2008
|
||||
- check_imo_is2008 — verificación IS Code 2008 Cap.2
|
||||
"""
|
||||
from arshipdesign.stability.gz_integrator import (
|
||||
GZPoint,
|
||||
GZCurve,
|
||||
compute_gz_wall_sided,
|
||||
compute_gz_direct,
|
||||
)
|
||||
from arshipdesign.stability.imo_is2008 import (
|
||||
IMOCriterion,
|
||||
IMOResult,
|
||||
check_imo_is2008,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"GZPoint",
|
||||
"GZCurve",
|
||||
"compute_gz_wall_sided",
|
||||
"compute_gz_direct",
|
||||
"IMOCriterion",
|
||||
"IMOResult",
|
||||
"check_imo_is2008",
|
||||
]
|
||||
|
||||
@@ -0,0 +1,579 @@
|
||||
"""
|
||||
gz_integrator.py — Motor de cálculo de la curva GZ de estabilidad estática.
|
||||
|
||||
Dos métodos:
|
||||
1. compute_gz_wall_sided — Fórmula de pared lateral (rápida, exacta hasta ~30°).
|
||||
2. compute_gz_direct — Integración directa 3-D con recorte de polígono
|
||||
por Sutherland-Hodgman (precisa a grandes ángulos).
|
||||
|
||||
Convenciones
|
||||
------------
|
||||
* Todos los ángulos de escora φ en grados en la API pública;
|
||||
en radianes internamente.
|
||||
* Sistema de coordenadas del casco (upright):
|
||||
y = manga (+ estribor), z = altura (0 = quilla, + hacia arriba).
|
||||
* Sistema de coordenadas del mundo (heeled):
|
||||
Eje z_w vertical apuntando arriba; eje y_w horizontal apuntando
|
||||
a estribor mundo.
|
||||
* La escora φ es positiva a estribor.
|
||||
|
||||
Fórmula GZ
|
||||
----------
|
||||
GZ = y_B_world − KG · sin(φ)
|
||||
donde y_B_world es la coordenada horizontal (eje y_w) del centro de
|
||||
carena B en el sistema del mundo.
|
||||
|
||||
Autor: Álvaro Romero
|
||||
Módulo 3 — AR-ShipDesign
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
from scipy.integrate import simpson as _scipy_simpson
|
||||
from scipy.optimize import brentq
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Dataclasses públicos
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
@dataclass
|
||||
class GZPoint:
|
||||
"""Un punto de la curva GZ."""
|
||||
phi_deg: float
|
||||
gz: float
|
||||
|
||||
|
||||
@dataclass
|
||||
class GZCurve:
|
||||
"""Curva GZ completa con estadísticas derivadas.
|
||||
|
||||
Atributos
|
||||
---------
|
||||
hull_name : str
|
||||
draft : float [m]
|
||||
kg : float [m] — Altura del centro de gravedad sobre la quilla.
|
||||
gm : float [m] — Altura metacéntrica transversal GM = KMT - KG.
|
||||
bmt : float [m] — Radio metacéntrico transversal.
|
||||
points : list[GZPoint]
|
||||
area_0_30 : float [m·rad]
|
||||
area_30_40 : float [m·rad]
|
||||
area_0_40 : float [m·rad]
|
||||
gz_max : float [m]
|
||||
phi_gz_max : float [°]
|
||||
avs : float [°] — Ángulo de estabilidad nula (Angle of Vanishing Stability).
|
||||
"""
|
||||
hull_name: str
|
||||
draft: float
|
||||
kg: float
|
||||
gm: float
|
||||
bmt: float
|
||||
points: list[GZPoint] = field(default_factory=list)
|
||||
area_0_30: float = 0.0
|
||||
area_30_40: float = 0.0
|
||||
area_0_40: float = 0.0
|
||||
gz_max: float = 0.0
|
||||
phi_gz_max: float = 0.0
|
||||
avs: float = 90.0
|
||||
|
||||
# ------------------------------------------------------------------
|
||||
# Propiedades de acceso rápido a arrays
|
||||
# ------------------------------------------------------------------
|
||||
|
||||
@property
|
||||
def angles_deg(self) -> np.ndarray:
|
||||
return np.array([p.phi_deg for p in self.points], dtype=float)
|
||||
|
||||
@property
|
||||
def gz_values(self) -> np.ndarray:
|
||||
return np.array([p.gz for p in self.points], dtype=float)
|
||||
|
||||
# ------------------------------------------------------------------
|
||||
# Cálculo de áreas, gz_max y AVS
|
||||
# ------------------------------------------------------------------
|
||||
|
||||
def _compute_areas(self) -> None:
|
||||
"""Integra las áreas bajo la curva GZ y calcula gz_max, phi_gz_max, avs.
|
||||
|
||||
Usa scipy.integrate.simpson sobre los puntos disponibles.
|
||||
Las áreas se calculan en [m·rad] (ángulos convertidos a radianes).
|
||||
"""
|
||||
phi = self.angles_deg
|
||||
gz = self.gz_values
|
||||
|
||||
if len(phi) < 2:
|
||||
return
|
||||
|
||||
# Convertir ángulos a radianes para la integración
|
||||
phi_rad = np.deg2rad(phi)
|
||||
|
||||
# ── Área 0–30° ──────────────────────────────────────────────
|
||||
mask_030 = phi <= 30.0
|
||||
if mask_030.sum() >= 2:
|
||||
self.area_0_30 = float(abs(_scipy_simpson(
|
||||
gz[mask_030], x=phi_rad[mask_030]
|
||||
)))
|
||||
else:
|
||||
self.area_0_30 = 0.0
|
||||
|
||||
# ── Área 30–40° ─────────────────────────────────────────────
|
||||
# Incluir el punto exacto en 30° (interpolado si no existe)
|
||||
gz30 = float(np.interp(30.0, phi, gz))
|
||||
gz40 = float(np.interp(40.0, phi, gz))
|
||||
|
||||
mask_3040 = (phi >= 30.0) & (phi <= 40.0)
|
||||
pts_3040_phi = np.concatenate([[30.0], phi[mask_3040], [40.0]])
|
||||
pts_3040_gz = np.concatenate([[gz30], gz[mask_3040], [gz40]])
|
||||
|
||||
# Eliminar duplicados
|
||||
_, unique_idx = np.unique(pts_3040_phi, return_index=True)
|
||||
pts_3040_phi = pts_3040_phi[unique_idx]
|
||||
pts_3040_gz = pts_3040_gz[unique_idx]
|
||||
|
||||
if len(pts_3040_phi) >= 2:
|
||||
self.area_30_40 = float(abs(_scipy_simpson(
|
||||
pts_3040_gz, x=np.deg2rad(pts_3040_phi)
|
||||
)))
|
||||
else:
|
||||
self.area_30_40 = 0.0
|
||||
|
||||
# ── Área 0–40° ──────────────────────────────────────────────
|
||||
mask_040 = phi <= 40.0
|
||||
pts_040_phi = np.concatenate([phi[mask_040], [40.0]])
|
||||
pts_040_gz = np.concatenate([gz[mask_040], [gz40]])
|
||||
_, unique_idx2 = np.unique(pts_040_phi, return_index=True)
|
||||
pts_040_phi = pts_040_phi[unique_idx2]
|
||||
pts_040_gz = pts_040_gz[unique_idx2]
|
||||
|
||||
if len(pts_040_phi) >= 2:
|
||||
self.area_0_40 = float(abs(_scipy_simpson(
|
||||
pts_040_gz, x=np.deg2rad(pts_040_phi)
|
||||
)))
|
||||
else:
|
||||
self.area_0_40 = 0.0
|
||||
|
||||
# ── GZ máximo y su ángulo ────────────────────────────────────
|
||||
idx_max = int(np.argmax(gz))
|
||||
self.gz_max = float(gz[idx_max])
|
||||
self.phi_gz_max = float(phi[idx_max])
|
||||
|
||||
# ── AVS: primer cruce de GZ = 0 después del pico ────────────
|
||||
# Buscar a partir del índice del máximo
|
||||
avs = 90.0
|
||||
for i in range(idx_max, len(gz) - 1):
|
||||
if gz[i] > 0.0 and gz[i + 1] <= 0.0:
|
||||
# Interpolación lineal para encontrar el cruce exacto
|
||||
phi1, gz1 = phi[i], gz[i]
|
||||
phi2, gz2 = phi[i + 1], gz[i + 1]
|
||||
if abs(gz2 - gz1) > 1e-12:
|
||||
avs = phi1 + (phi2 - phi1) * (-gz1) / (gz2 - gz1)
|
||||
else:
|
||||
avs = float(phi1)
|
||||
break
|
||||
self.avs = avs
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Método 1: Fórmula de pared lateral (Wall-Sided)
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
def compute_gz_wall_sided(
|
||||
hull,
|
||||
draft: float,
|
||||
kg: Optional[float] = None,
|
||||
angles_deg: Optional[np.ndarray] = None,
|
||||
rho: float = 1025.0,
|
||||
) -> GZCurve:
|
||||
"""Calcula la curva GZ por la fórmula de pared lateral.
|
||||
|
||||
GZ = sin(φ) · (GM + 0.5 · BM · tan²(φ))
|
||||
|
||||
Esta fórmula es exacta para cascos de paredes verticales (wall-sided)
|
||||
a cualquier ángulo; para cascos con flare proporciona buena aproximación
|
||||
hasta ~30° y subestima a ángulos mayores.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
hull : Hull
|
||||
draft : float [m]
|
||||
kg : float, optional — Si None, usa hull.depth * 0.55
|
||||
angles_deg : array_like, optional — Default: 0..90 en pasos de 1°
|
||||
rho : float — Densidad del agua [kg/m³]
|
||||
|
||||
Returns
|
||||
-------
|
||||
GZCurve
|
||||
"""
|
||||
from arshipdesign.hydrostatics.upright import compute_upright
|
||||
|
||||
hydro = compute_upright(hull, draft, rho=rho)
|
||||
|
||||
kg_val = float(hull.depth * 0.55 if kg is None else kg)
|
||||
gm = hydro.kmt - kg_val
|
||||
bmt = hydro.bmt
|
||||
|
||||
if angles_deg is None:
|
||||
phi_arr = np.linspace(0.0, 90.0, 91)
|
||||
else:
|
||||
phi_arr = np.asarray(angles_deg, dtype=float)
|
||||
|
||||
points: list[GZPoint] = []
|
||||
for phi_deg in phi_arr:
|
||||
phi_rad = math.radians(phi_deg)
|
||||
sin_phi = math.sin(phi_rad)
|
||||
tan_phi = math.tan(phi_rad) if abs(phi_rad) < math.radians(89.9) else math.tan(math.radians(89.9))
|
||||
gz = sin_phi * (gm + 0.5 * bmt * tan_phi ** 2)
|
||||
points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz)))
|
||||
|
||||
curve = GZCurve(
|
||||
hull_name=hull.name,
|
||||
draft=float(draft),
|
||||
kg=kg_val,
|
||||
gm=gm,
|
||||
bmt=bmt,
|
||||
points=points,
|
||||
)
|
||||
curve._compute_areas()
|
||||
return curve
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
# Método 2: Integración directa (Direct)
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
def _clip_polygon_below_z(
|
||||
poly: list[tuple[float, float]],
|
||||
z_wl: float,
|
||||
) -> list[tuple[float, float]]:
|
||||
"""Sutherland-Hodgman clip: mantiene región z ≤ z_wl.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
poly : list of (y, z)
|
||||
z_wl : float — línea de corte (cota de flotación en el sistema girado)
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of (y, z) — polígono recortado (puede ser vacío)
|
||||
"""
|
||||
if not poly:
|
||||
return []
|
||||
|
||||
output = list(poly)
|
||||
|
||||
def _inside(pt: tuple[float, float]) -> bool:
|
||||
return pt[1] <= z_wl
|
||||
|
||||
def _intersect(
|
||||
a: tuple[float, float], b: tuple[float, float]
|
||||
) -> tuple[float, float]:
|
||||
"""Intersección de segmento (a→b) con z = z_wl."""
|
||||
ya, za = a
|
||||
yb, zb = b
|
||||
dz = zb - za
|
||||
if abs(dz) < 1e-14:
|
||||
return ((ya + yb) * 0.5, z_wl)
|
||||
t = (z_wl - za) / dz
|
||||
return (ya + t * (yb - ya), z_wl)
|
||||
|
||||
# Plano de recorte: z ≤ z_wl (mantener interior)
|
||||
input_list = output
|
||||
output = []
|
||||
|
||||
n = len(input_list)
|
||||
for i in range(n):
|
||||
curr = input_list[i]
|
||||
prev = input_list[i - 1]
|
||||
if _inside(curr):
|
||||
if not _inside(prev):
|
||||
output.append(_intersect(prev, curr))
|
||||
output.append(curr)
|
||||
elif _inside(prev):
|
||||
output.append(_intersect(prev, curr))
|
||||
|
||||
return output
|
||||
|
||||
|
||||
def _polygon_area_centroid(
|
||||
poly: list[tuple[float, float]],
|
||||
) -> tuple[float, float, float]:
|
||||
"""Área y centroide de un polígono cerrado por la fórmula de Shoelace.
|
||||
|
||||
Returns
|
||||
-------
|
||||
area : float — Área absoluta (siempre ≥ 0).
|
||||
yc : float — Coordenada y del centroide.
|
||||
zc : float — Coordenada z del centroide.
|
||||
"""
|
||||
n = len(poly)
|
||||
if n < 3:
|
||||
if n == 0:
|
||||
return 0.0, 0.0, 0.0
|
||||
y0, z0 = poly[0]
|
||||
return 0.0, float(y0), float(z0)
|
||||
|
||||
# Shoelace
|
||||
ys = np.array([p[0] for p in poly], dtype=float)
|
||||
zs = np.array([p[1] for p in poly], dtype=float)
|
||||
|
||||
# Área con signo
|
||||
cross = ys[:-1] * zs[1:] - ys[1:] * zs[:-1]
|
||||
# Cerrar el polígono
|
||||
cross_close = ys[-1] * zs[0] - ys[0] * zs[-1]
|
||||
total_cross = np.sum(cross) + cross_close
|
||||
area_signed = 0.5 * total_cross
|
||||
area = abs(area_signed)
|
||||
|
||||
if area < 1e-15:
|
||||
return 0.0, float(np.mean(ys)), float(np.mean(zs))
|
||||
|
||||
# Centroides
|
||||
cx_terms = (ys[:-1] + ys[1:]) * cross
|
||||
cz_terms = (zs[:-1] + zs[1:]) * cross
|
||||
cx_close = (ys[-1] + ys[0]) * cross_close
|
||||
cz_close = (zs[-1] + zs[0]) * cross_close
|
||||
|
||||
yc = (np.sum(cx_terms) + cx_close) / (6.0 * area_signed)
|
||||
zc = (np.sum(cz_terms) + cz_close) / (6.0 * area_signed)
|
||||
|
||||
return area, float(yc), float(zc)
|
||||
|
||||
|
||||
def _build_section_polygon_extended(
|
||||
half_breadths: np.ndarray,
|
||||
z_positions: np.ndarray,
|
||||
z_extend: float,
|
||||
) -> list[tuple[float, float]]:
|
||||
"""Construye el polígono extendido de una sección transversal.
|
||||
|
||||
Para capturar correctamente el cálculo de GZ a grandes ángulos,
|
||||
el polígono se extiende con paredes verticales (wall-sided) desde
|
||||
la línea de flotación upright hasta z_extend. Esto garantiza que
|
||||
la integral de Sutherland-Hodgman captura el cuño de inmersión /
|
||||
emersión al escorar.
|
||||
|
||||
Geometría:
|
||||
- Estribor: de quilla (j=0) hacia arriba hasta el tope de la sección,
|
||||
luego pared vertical hasta z_extend.
|
||||
- Babor: pared vertical de z_extend hacia abajo, luego espejo de estribor
|
||||
hasta la quilla.
|
||||
- Resultado: 2n+2 vértices (n puntos de sección × 2 bandas + 2 extensiones).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
half_breadths : ndarray, shape (n,)
|
||||
z_positions : ndarray, shape (n,) — creciente desde la quilla
|
||||
z_extend : float — altura máxima del polígono (m)
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of (y, z) — polígono cerrado en sentido antihorario
|
||||
"""
|
||||
n = len(z_positions)
|
||||
y_top = float(half_breadths[-1])
|
||||
z_top = float(z_positions[-1]) # normalmente = draft
|
||||
|
||||
poly: list[tuple[float, float]] = []
|
||||
# Estribor: de quilla (j=0) hacia arriba
|
||||
for j in range(n):
|
||||
poly.append((float(half_breadths[j]), float(z_positions[j])))
|
||||
# Extensión estribor hasta z_extend (wall-sided)
|
||||
if z_extend > z_top + 1e-9:
|
||||
poly.append((y_top, z_extend))
|
||||
# Extensión babor desde z_extend hacia abajo (wall-sided, mirror)
|
||||
if z_extend > z_top + 1e-9:
|
||||
poly.append((-y_top, z_extend))
|
||||
# Babor: de tope hacia abajo (mirror)
|
||||
for j in reversed(range(n)):
|
||||
poly.append((-float(half_breadths[j]), float(z_positions[j])))
|
||||
return poly
|
||||
|
||||
|
||||
def _compute_polygon_volume(
|
||||
upright_polys: list[list[tuple[float, float]]],
|
||||
dx: np.ndarray,
|
||||
z_wl: float,
|
||||
) -> float:
|
||||
"""Suma del volumen de los polígonos recortados en z ≤ z_wl."""
|
||||
vol = 0.0
|
||||
for i, poly in enumerate(upright_polys):
|
||||
clipped = _clip_polygon_below_z(poly, z_wl)
|
||||
if clipped:
|
||||
a, _, _ = _polygon_area_centroid(clipped)
|
||||
vol += a * dx[i]
|
||||
return vol
|
||||
|
||||
|
||||
def compute_gz_direct(
|
||||
hull,
|
||||
draft: float,
|
||||
kg: Optional[float] = None,
|
||||
angles_deg: Optional[np.ndarray] = None,
|
||||
rho: float = 1025.0,
|
||||
) -> GZCurve:
|
||||
"""Calcula la curva GZ por integración directa de secciones.
|
||||
|
||||
Para cada ángulo de escora φ (escora a estribor, φ > 0):
|
||||
1. Construye el polígono extendido (wall-sided por encima del calado)
|
||||
de cada sección para capturar los cuños de inmersión/emersión.
|
||||
2. Aplica la rotación de casco-a-mundo (giro en sentido horario=escora
|
||||
a estribor):
|
||||
y_w = y·cos(φ) + z·sin(φ)
|
||||
z_w = -y·sin(φ) + z·cos(φ)
|
||||
3. Determina z_wl (posición de la flotación en el mundo) por bisección
|
||||
(brentq) tal que el volumen sumergido = V_target.
|
||||
4. Calcula y_B_world = centroide horizontal del volumen recortado.
|
||||
5. GZ = y_B_world − KG·sin(φ)
|
||||
|
||||
La rotación utilizada es la convención de estabilidad naval estándar:
|
||||
+φ = escora a estribor → el lado estribor (y_body > 0) baja.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
hull : Hull
|
||||
draft : float [m]
|
||||
kg : float, optional — Si None, usa hull.depth * 0.55
|
||||
angles_deg : array_like, optional — Default: 0..90 en pasos de 1°
|
||||
rho : float
|
||||
|
||||
Returns
|
||||
-------
|
||||
GZCurve
|
||||
"""
|
||||
from arshipdesign.hydrostatics.upright import compute_upright
|
||||
|
||||
hydro = compute_upright(hull, draft, rho=rho)
|
||||
bmt = hydro.bmt
|
||||
kmt = hydro.kmt
|
||||
kg_val = float(hull.depth * 0.55 if kg is None else kg)
|
||||
gm = kmt - kg_val
|
||||
|
||||
if angles_deg is None:
|
||||
phi_arr = np.linspace(0.0, 90.0, 91)
|
||||
else:
|
||||
phi_arr = np.asarray(angles_deg, dtype=float)
|
||||
|
||||
# ── Preparar secciones y pesos de integración trapezoidal ────────
|
||||
sections = hull.offsets.to_sections()
|
||||
x_sta = np.array([s.x for s in sections], dtype=float)
|
||||
n_sta = len(x_sta)
|
||||
|
||||
dx = np.zeros(n_sta, dtype=float)
|
||||
if n_sta >= 2:
|
||||
dx[0] = (x_sta[1] - x_sta[0]) * 0.5
|
||||
dx[-1] = (x_sta[-1] - x_sta[-2]) * 0.5
|
||||
for i in range(1, n_sta - 1):
|
||||
dx[i] = (x_sta[i + 1] - x_sta[i - 1]) * 0.5
|
||||
|
||||
# ── Polígonos extendidos en posición upright ──────────────────────
|
||||
# Se extienden con paredes verticales hasta z_extend para capturar
|
||||
# el cuño de inmersión a grandes ángulos de escora.
|
||||
# z_extend debe ser suficientemente mayor que draft + beam/2 * sin(90°).
|
||||
z_extend = float(draft) + float(hull.beam) * 1.5 + float(hull.depth) * 0.5
|
||||
upright_polys: list[list[tuple[float, float]]] = []
|
||||
for sec in sections:
|
||||
poly = _build_section_polygon_extended(
|
||||
sec.half_breadths, sec.z_positions, z_extend
|
||||
)
|
||||
upright_polys.append(poly)
|
||||
|
||||
# ── V_target: volumen coherente con los polígonos ─────────────────
|
||||
# Usar el volumen que los polígonos dan al cortarlos en z ≤ draft,
|
||||
# para garantizar coherencia numérica entre phi=0 y phi>0.
|
||||
V_target = _compute_polygon_volume(upright_polys, dx, float(draft))
|
||||
|
||||
# Límites conservadores para la bisección
|
||||
z_bisect_high = z_extend * 0.95 # justo por debajo del tope extendido
|
||||
|
||||
points: list[GZPoint] = []
|
||||
|
||||
for phi_deg in phi_arr:
|
||||
phi_rad = math.radians(float(phi_deg))
|
||||
cos_phi = math.cos(phi_rad)
|
||||
sin_phi = math.sin(phi_rad)
|
||||
|
||||
if abs(phi_rad) < 1e-9:
|
||||
# φ = 0: GZ = 0 exacto por simetría
|
||||
points.append(GZPoint(phi_deg=float(phi_deg), gz=0.0))
|
||||
continue
|
||||
|
||||
# ── Rotación de cuerpo a mundo para escora a estribor ────────
|
||||
# Convención naval: +φ → estribor baja (giro horario en plano y-z)
|
||||
# y_w = y·cos(φ) + z·sin(φ)
|
||||
# z_w = -y·sin(φ) + z·cos(φ)
|
||||
rotated_polys: list[list[tuple[float, float]]] = []
|
||||
for poly in upright_polys:
|
||||
world_poly = [
|
||||
( y * cos_phi + z * sin_phi,
|
||||
-y * sin_phi + z * cos_phi)
|
||||
for y, z in poly
|
||||
]
|
||||
rotated_polys.append(world_poly)
|
||||
|
||||
# ── Bisección para z_wl en el sistema del mundo ───────────────
|
||||
def _vol_minus_target(z_wl: float) -> float:
|
||||
vol = 0.0
|
||||
for i in range(n_sta):
|
||||
clipped = _clip_polygon_below_z(rotated_polys[i], z_wl)
|
||||
if clipped:
|
||||
a, _, _ = _polygon_area_centroid(clipped)
|
||||
vol += a * dx[i]
|
||||
return vol - V_target
|
||||
|
||||
z_b_low = 1e-4
|
||||
z_b_high = z_bisect_high
|
||||
|
||||
f_low = _vol_minus_target(z_b_low)
|
||||
f_high = _vol_minus_target(z_b_high)
|
||||
|
||||
if f_low * f_high > 0:
|
||||
# Fallback: fórmula wall-sided
|
||||
phi_clamp = min(phi_rad, math.radians(89.0))
|
||||
gz_ws = sin_phi * (gm + 0.5 * bmt * math.tan(phi_clamp) ** 2)
|
||||
points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz_ws)))
|
||||
continue
|
||||
|
||||
try:
|
||||
z_wl = brentq(_vol_minus_target, z_b_low, z_b_high,
|
||||
xtol=1e-5, maxiter=150)
|
||||
except ValueError:
|
||||
phi_clamp = min(phi_rad, math.radians(89.0))
|
||||
gz_ws = sin_phi * (gm + 0.5 * bmt * math.tan(phi_clamp) ** 2)
|
||||
points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz_ws)))
|
||||
continue
|
||||
|
||||
# ── Centroide horizontal de carena en el mundo ────────────────
|
||||
vol_total = 0.0
|
||||
moment_y = 0.0
|
||||
|
||||
for i in range(n_sta):
|
||||
clipped = _clip_polygon_below_z(rotated_polys[i], z_wl)
|
||||
if clipped:
|
||||
a, yc, _ = _polygon_area_centroid(clipped)
|
||||
vol_total += a * dx[i]
|
||||
moment_y += a * yc * dx[i]
|
||||
|
||||
if vol_total > 1e-12:
|
||||
y_B_world = moment_y / vol_total
|
||||
else:
|
||||
y_B_world = 0.0
|
||||
|
||||
# GZ = y_B_world - y_G_world
|
||||
# G = (y_body=0, z_body=KG); en mundo (CW): y_G_w = 0*cos + KG*sin = KG*sin
|
||||
gz = y_B_world - kg_val * sin_phi
|
||||
points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz)))
|
||||
|
||||
curve = GZCurve(
|
||||
hull_name=hull.name,
|
||||
draft=float(draft),
|
||||
kg=kg_val,
|
||||
gm=gm,
|
||||
bmt=bmt,
|
||||
points=points,
|
||||
)
|
||||
curve._compute_areas()
|
||||
return curve
|
||||
@@ -0,0 +1,152 @@
|
||||
"""
|
||||
imo_is2008.py — IMO IS Code 2008, Capítulo 2 — criterios de estabilidad intacta.
|
||||
|
||||
Criterios A.2.1 para buques de carga general y pequeñas embarcaciones:
|
||||
2.1.1 Área 0–30° ≥ 0.055 m·rad
|
||||
2.1.2 Área 0–40° ≥ 0.090 m·rad
|
||||
2.1.3 Área 30–40° ≥ 0.030 m·rad
|
||||
2.1.4 GZ a 30° ≥ 0.200 m
|
||||
2.1.5 Ángulo de GZ máximo ≥ 25°
|
||||
2.1.6 GM₀ ≥ 0.150 m
|
||||
|
||||
Referencia: IMO MSC.267(85) — IS Code 2008, Parte A, Cap. 2.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
|
||||
@dataclass
|
||||
class IMOCriterion:
|
||||
"""Un criterio individual del IS Code 2008."""
|
||||
code: str # e.g. "A.2.1.1"
|
||||
description: str # Descripción corta
|
||||
required: float # Valor mínimo requerido
|
||||
achieved: float # Valor obtenido de la curva GZ
|
||||
unit: str # Unidades (m·rad, m, °)
|
||||
passed: bool # True si achieved >= required
|
||||
|
||||
|
||||
@dataclass
|
||||
class IMOResult:
|
||||
"""Resultado completo de la verificación IMO IS Code 2008."""
|
||||
criteria: list[IMOCriterion]
|
||||
overall_passed: bool
|
||||
|
||||
def table_rows(self) -> list[tuple[str, str, str, str, bool]]:
|
||||
"""Devuelve filas para la tabla: (code, description, required_str, achieved_str, passed).
|
||||
|
||||
El formato de los strings varía según las unidades:
|
||||
- m·rad: 4 decimales
|
||||
- m: 3 decimales
|
||||
- °: 1 decimal
|
||||
"""
|
||||
rows = []
|
||||
for c in self.criteria:
|
||||
if c.unit == "m·rad":
|
||||
req_str = f"{c.required:.4f} {c.unit}"
|
||||
ach_str = f"{c.achieved:.4f} {c.unit}"
|
||||
elif c.unit == "m":
|
||||
req_str = f"{c.required:.3f} {c.unit}"
|
||||
ach_str = f"{c.achieved:.3f} {c.unit}"
|
||||
elif c.unit == "°":
|
||||
req_str = f"{c.required:.1f}{c.unit}"
|
||||
ach_str = f"{c.achieved:.1f}{c.unit}"
|
||||
else:
|
||||
req_str = f"{c.required} {c.unit}"
|
||||
ach_str = f"{c.achieved} {c.unit}"
|
||||
rows.append((c.code, c.description, req_str, ach_str, c.passed))
|
||||
return rows
|
||||
|
||||
|
||||
def check_imo_is2008(gz) -> IMOResult:
|
||||
"""Verifica todos los criterios IS Code 2008 Cap.2 para la curva GZ dada.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
gz : GZCurve
|
||||
Curva de estabilidad calculada.
|
||||
|
||||
Returns
|
||||
-------
|
||||
IMOResult
|
||||
Contiene la lista de criterios individuales y el resultado global.
|
||||
"""
|
||||
import numpy as np
|
||||
from arshipdesign.stability.gz_integrator import GZCurve
|
||||
|
||||
def _criterion(
|
||||
code: str,
|
||||
desc: str,
|
||||
req: float,
|
||||
ach: float,
|
||||
unit: str,
|
||||
) -> IMOCriterion:
|
||||
return IMOCriterion(
|
||||
code=code,
|
||||
description=desc,
|
||||
required=req,
|
||||
achieved=ach,
|
||||
unit=unit,
|
||||
passed=(ach >= req),
|
||||
)
|
||||
|
||||
criteria: list[IMOCriterion] = []
|
||||
|
||||
# A.2.1.1 — Área bajo la curva GZ entre 0° y 30°
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.1",
|
||||
"Área 0–30°",
|
||||
0.055,
|
||||
float(gz.area_0_30),
|
||||
"m·rad",
|
||||
))
|
||||
|
||||
# A.2.1.2 — Área bajo la curva GZ entre 0° y 40°
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.2",
|
||||
"Área 0–40°",
|
||||
0.090,
|
||||
float(gz.area_0_40),
|
||||
"m·rad",
|
||||
))
|
||||
|
||||
# A.2.1.3 — Área bajo la curva GZ entre 30° y 40°
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.3",
|
||||
"Área 30–40°",
|
||||
0.030,
|
||||
float(gz.area_30_40),
|
||||
"m·rad",
|
||||
))
|
||||
|
||||
# A.2.1.4 — GZ a 30° de escora
|
||||
gz30 = float(np.interp(30.0, gz.angles_deg, gz.gz_values))
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.4",
|
||||
"GZ a 30°",
|
||||
0.200,
|
||||
gz30,
|
||||
"m",
|
||||
))
|
||||
|
||||
# A.2.1.5 — Ángulo en que se produce el GZ máximo
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.5",
|
||||
"Ángulo GZ máx",
|
||||
25.0,
|
||||
float(gz.phi_gz_max),
|
||||
"°",
|
||||
))
|
||||
|
||||
# A.2.1.6 — Altura metacéntrica inicial GM₀
|
||||
criteria.append(_criterion(
|
||||
"A.2.1.6",
|
||||
"GM₀",
|
||||
0.150,
|
||||
float(gz.gm),
|
||||
"m",
|
||||
))
|
||||
|
||||
overall = all(c.passed for c in criteria)
|
||||
return IMOResult(criteria=criteria, overall_passed=overall)
|
||||
Reference in New Issue
Block a user