From d14efa3da59943eab0fc22bc246e38cbbdd6acbd Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Mon, 20 Jul 2026 00:34:43 -0700 Subject: [PATCH 01/21] Add the area and CE to CG coordinate values --- .../boat_simulator/common/constants.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index 2f9e5d5a1..870216e44 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -138,20 +138,20 @@ class BoatProperties: WING_SAIL_CHORD = 1.5 # Units: meters # TODO Placeholder: measure the distance from the mast axis to the tab's aero center. -WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH = 1.0 # Units: meters +WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH = 1.43 # Units: meters # TODO sail_dist is the sail CE-to-pivot distance, not CE-to-CG; z_s (CE height # relative to the CG) is a placeholder until we have real geometry. -CE_HEIGHT_REL_TO_CG = -0.5 # Units: meters +CE_HEIGHT_REL_TO_CG = -0.92 # Units: meters # TODO Placeholder: measure the mast pivot's chordwise position (~25% chord assumed). MAST_PIVOT_CHORD_FRACTION = 0.25 # Fraction of the wing chord, dimensionless # TODO Placeholder: measure the rudder's center of effort depth below the CG. -RUDDER_CE_DEPTH_REL_TO_CG = 0.5 # z_r, units: meters +RUDDER_CE_DEPTH_REL_TO_CG = 0.74 # z_r, units: meters # TODO Placeholder: measure the keel's center of effort relative to the CG. -KEEL_CE_REL_TO_CG = (0.0, 0.5) # (x_k, z_k), units: meters +KEEL_CE_REL_TO_CG = (0.0, -0.4) # (x_k, z_k), units: meters # TODO Placeholder: measure the hull's center of effort relative to the CG. HULL_CE_REL_TO_CG = (0.0, 0.0, 0.0) # (x_h, y_h, z_h), units: meters @@ -160,7 +160,7 @@ class BoatProperties: HULL_LINEAR_DRAG = 0.0 # Units: newton seconds per meter # Displaced volume of the boat at floating equilibrium (m^3) -DISPLACED_VOLUME = 0.1 +DISPLACED_VOLUME = 0.05 # Constants related to the physical and mechanical properties of Polaris # TODO These are placeholder values which should be replaced when we have real values. @@ -204,7 +204,7 @@ class BoatProperties: dtype=np.float64, ) ), - sail_areas=20.0, # meters ^ 2 + sail_areas=2.01, # meters ^ 2 # TODO: Replace the below placeholder constants with the real values/approximates tab_lift_coeffs=CoeffTable( np.array( @@ -244,7 +244,7 @@ class BoatProperties: dtype=np.float64, ) ), - tab_areas=5.0, # meters ^ 2 + tab_areas=0.198, # meters ^ 2 # Rudder: ±45° → table covers 0–45° (sign handled by caller) # NACA symmetric foil: stalls ~20–22° rudder_lift_coeffs=CoeffTable( @@ -281,7 +281,7 @@ class BoatProperties: dtype=np.float64, ) ), - rudder_areas=0.4, + rudder_areas=0.117, # meters ^ 2 # TODO: Replace the below placeholder constants with the real values/approximates keel_lift_coeffs=CoeffTable( @@ -322,13 +322,13 @@ class BoatProperties: dtype=np.float64, ) ), - keel_areas=2.0, # meters ^ 2 + keel_areas=0.51, # meters ^ 2 sail_dist=0.5, # meters rudder_dist=1.0, # meters hull_drag_factor=0.5, - mass=225.0, + mass=276.0, # M_RB = diag(m, m, I_xx, I_zz) — surge, sway, roll, yaw - inertia=Mat4(np.diag([225.0, 225.0, 125.0, 500.0])), + inertia=Mat4(np.diag([225.0, 225.0, 175.86, 119.04])), # TODO: Replace with coefficients from strip theory or CFD for the real hull. # M_A = diag(X_u̇, Y_v̇, K_ṗ, N_ṙ) — diagonal linearized added-mass approximation, # estimated from slender-hull rules of thumb relative to M_RB: From c448b7cad6b97a59c5ba83fa0cef44e17cbc7315 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Mon, 20 Jul 2026 10:24:31 -0700 Subject: [PATCH 02/21] Added all the constants except for keel's z_k for boat_simulator stability --- .../boat_simulator/common/constants.py | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index 870216e44..363eef64e 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -134,24 +134,25 @@ class BoatProperties: # TODO: This is a placeholder value for the metacentric height METACENTRIC_HEIGHT = 0.3 # Units: meters -# TODO Placeholder: derive the mean chord from the real wingsail geometry. +# Derive the mean chord from the real wingsail geometry. WING_SAIL_CHORD = 1.5 # Units: meters -# TODO Placeholder: measure the distance from the mast axis to the tab's aero center. +# Measure the distance from the mast axis to the tab's aero center. WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH = 1.43 # Units: meters -# TODO sail_dist is the sail CE-to-pivot distance, not CE-to-CG; z_s (CE height +# sail_dist is the sail CE-to-pivot distance, not CE-to-CG; z_s (CE height # relative to the CG) is a placeholder until we have real geometry. -CE_HEIGHT_REL_TO_CG = -0.92 # Units: meters +CE_HEIGHT_REL_TO_CG = -1.86 # Units: meters # TODO Placeholder: measure the mast pivot's chordwise position (~25% chord assumed). MAST_PIVOT_CHORD_FRACTION = 0.25 # Fraction of the wing chord, dimensionless -# TODO Placeholder: measure the rudder's center of effort depth below the CG. +# Measure the rudder's center of effort depth below the CG. RUDDER_CE_DEPTH_REL_TO_CG = 0.74 # z_r, units: meters -# TODO Placeholder: measure the keel's center of effort relative to the CG. -KEEL_CE_REL_TO_CG = (0.0, -0.4) # (x_k, z_k), units: meters +# NOTE Measure the keel's center of effort relative to the CG. The z_k is a magic number to +# maintain stability +KEEL_CE_REL_TO_CG = (0.08, -0.4) # (x_k, z_k), units: meters # TODO Placeholder: measure the hull's center of effort relative to the CG. HULL_CE_REL_TO_CG = (0.0, 0.0, 0.0) # (x_h, y_h, z_h), units: meters @@ -324,7 +325,7 @@ class BoatProperties: ), keel_areas=0.51, # meters ^ 2 sail_dist=0.5, # meters - rudder_dist=1.0, # meters + rudder_dist=1.35, # meters hull_drag_factor=0.5, mass=276.0, # M_RB = diag(m, m, I_xx, I_zz) — surge, sway, roll, yaw @@ -336,8 +337,8 @@ class BoatProperties: # Y_v̇ ≈ 0.80–1.00·m (sway drags a large volume of water; keel adds to it) # K_ṗ ≈ 0.20–0.30·I_xx # N_ṙ ≈ 0.50–0.70·I_zz - M_A=Mat4(np.diag([20.0, 180.0, 30.0, 300.0])), + M_A=Mat4(np.diag([20.0, 400.0, 230.0, 300.0])), # TODO: Replace with real damping coefficients from tow-tank or CFD data. # D = diag(X_u, Y_v, K_p, N_r) — linear damping per DOF. - D=Mat4(np.diag([15.0, 80.0, 25.0, 60.0])), + D=Mat4(np.diag([45.0, 180.0, 115.0, 140.0])), ) From c261a139d3333711bdc9ac9d4e9e73044730bda4 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Wed, 22 Jul 2026 22:09:53 -0700 Subject: [PATCH 03/21] Add metacentric height and refactoring --- .../boat_simulator/common/constants.py | 32 +++++----- .../nodes/physics_engine/fluid_forces.py | 61 +++++++++++++++---- 2 files changed, 65 insertions(+), 28 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index 363eef64e..cd95df1f3 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -61,10 +61,6 @@ class BoatProperties: keel_drag_coeffs: CoeffTable # Float: each row is keel_area_m2. keel_areas: float - # Distance from the center of effort of the sail to the pivot point (m). - sail_dist: float - # Distance from the center of effort of the rudder to the pivot point (m). - rudder_dist: float # Dimensionless quadratic drag coefficient for the hull: F_drag = hull_drag_factor * |v| * v. hull_drag_factor: float # Mass of the boat (kg). @@ -131,8 +127,11 @@ class BoatProperties: # Gravity in m/s EARTH_GRAVITY = 9.81 -# TODO: This is a placeholder value for the metacentric height -METACENTRIC_HEIGHT = 0.3 # Units: meters +# TODO: Center of Bouyancy guess. Please provide the correct point at Equilibrium +CoB_REL_COORD = 0.5 + +# Metacentric height used in testing +METACENTRIC_HEIGHT = 0.5 # Units: meters # Derive the mean chord from the real wingsail geometry. WING_SAIL_CHORD = 1.5 # Units: meters @@ -140,17 +139,16 @@ class BoatProperties: # Measure the distance from the mast axis to the tab's aero center. WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH = 1.43 # Units: meters -# sail_dist is the sail CE-to-pivot distance, not CE-to-CG; z_s (CE height -# relative to the CG) is a placeholder until we have real geometry. -CE_HEIGHT_REL_TO_CG = -1.86 # Units: meters - -# TODO Placeholder: measure the mast pivot's chordwise position (~25% chord assumed). +# The mast pivot's chordwise position (~25% chord assumed). MAST_PIVOT_CHORD_FRACTION = 0.25 # Fraction of the wing chord, dimensionless # Measure the rudder's center of effort depth below the CG. -RUDDER_CE_DEPTH_REL_TO_CG = 0.74 # z_r, units: meters +RUDDER_CE_REL_TO_CG = (-1.35, 0.74) -# NOTE Measure the keel's center of effort relative to the CG. The z_k is a magic number to +# The wingsail CE-to-CG in (x, Z) coordinate with units: meters +SAIL_CE_REL_TO_CG = (0.25, -1.86) # (x_s, z_s), units: meters + +# The keel's center of effort relative to the CG. The z_k is a magic number to # maintain stability KEEL_CE_REL_TO_CG = (0.08, -0.4) # (x_k, z_k), units: meters @@ -160,8 +158,6 @@ class BoatProperties: # TODO Placeholder: measure the hull's linear drag coefficient. HULL_LINEAR_DRAG = 0.0 # Units: newton seconds per meter -# Displaced volume of the boat at floating equilibrium (m^3) -DISPLACED_VOLUME = 0.05 # Constants related to the physical and mechanical properties of Polaris # TODO These are placeholder values which should be replaced when we have real values. @@ -324,8 +320,6 @@ class BoatProperties: ) ), keel_areas=0.51, # meters ^ 2 - sail_dist=0.5, # meters - rudder_dist=1.35, # meters hull_drag_factor=0.5, mass=276.0, # M_RB = diag(m, m, I_xx, I_zz) — surge, sway, roll, yaw @@ -342,3 +336,7 @@ class BoatProperties: # D = diag(X_u, Y_v, K_p, N_r) — linear damping per DOF. D=Mat4(np.diag([45.0, 180.0, 115.0, 140.0])), ) + + +# Displaced volume of the boat at floating equilibrium (m^3) +DISPLACED_VOLUME = BOAT_PROPERTIES.mass / WATER_DENSITY diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index dc84be29c..05dc88dcf 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -10,18 +10,18 @@ from boat_simulator.common.constants import ( AIR_DENSITY, BOAT_PROPERTIES, - CE_HEIGHT_REL_TO_CG, DISPLACED_VOLUME, EARTH_GRAVITY, HULL_CE_REL_TO_CG, HULL_LINEAR_DRAG, KEEL_CE_REL_TO_CG, MAST_PIVOT_CHORD_FRACTION, - METACENTRIC_HEIGHT, - RUDDER_CE_DEPTH_REL_TO_CG, + RUDDER_CE_REL_TO_CG, + SAIL_CE_REL_TO_CG, WATER_DENSITY, WING_SAIL_CHORD, WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH, + CoB_REL_COORD, ) from boat_simulator.common.conventions import ( Body, @@ -229,7 +229,6 @@ def __init__(self): self.__seawater_density = WATER_DENSITY self.__gravity = EARTH_GRAVITY self.__displaced_volume = DISPLACED_VOLUME - self.__metacentric_height = METACENTRIC_HEIGHT def compute(self, roll_angle_rad: float) -> Vec4[Force, Body]: """Computes the hydrostatic restoring force at the given roll angle @@ -240,18 +239,62 @@ def compute(self, roll_angle_rad: float) -> Vec4[Force, Body]: Returns: Vec4[Force, Body]: The restoring force """ + gm = self.calculate_metacentric_to_cog(roll_angle_rad) k_restore = ( -self.__seawater_density * self.__gravity * self.__displaced_volume - * self.__metacentric_height + * gm * math.sin(roll_angle_rad) ) _logger.info( f"HydroStatics.compute: roll={roll_angle_rad:.4f} rad K_restore={k_restore:.2f} N·m" ) + if not (-15.0 < roll_angle_rad < 15.0): + _logger.warning( + f"HydroStatics.compute: roll angle is outside the small angle approximation" + ) return Vec4.from_xypr(0.0, 0.0, k_restore, 0.0) + def calculate_metacentric_to_cog(self, roll_angle_rad: float) -> float: + """Computes GM, the distance from the CoG to the metacenter, by rotating the + (body-fixed) CoG and CoB with the hull and finding where the buoyancy line + crosses the hull's centerline. + + 1. Perpendicular to the equilibrium line (y = 0) is x = 0: at zero heel the + buoyant force, which always acts vertically, runs straight up the hull's + centerline through the CoG. + 2. Rotate the CoG (the body-frame origin) and the CoB (CoB_REL_COORD below the + CoG at equilibrium) by the roll angle, since both are fixed points on the + rigid hull. + 3. The buoyant force's line of action is the vertical line through the rotated + CoB (x = x_CoB_rotated). The CoG and CoB lie on the same body-frame line + (the centerline), so the rotated CoG-CoB line already passes through that + vertical line at the rotated CoB itself: that is the metacenter M. + 4. GM is the distance from the rotated CoG to M. + + Args: + roll_angle_rad (float): The boat's current roll angle, in radians. + + Returns: + float: GM, in meters. + """ + rotation = np.array( + [ + [math.cos(roll_angle_rad), -math.sin(roll_angle_rad)], + [math.sin(roll_angle_rad), math.cos(roll_angle_rad)], + ] + ) + + cog_body = np.array([0.0, 0.0]) + cob_body = np.array([0.0, -CoB_REL_COORD]) + + cog_rotated = rotation @ cog_body + metacenter = rotation @ cob_body + + gm = float(np.linalg.norm(metacenter - cog_rotated)) + return gm + class AeroDynamicsForceComputation: """Computes the wingsail's generalized force and moment contribution. @@ -288,8 +331,7 @@ def __init__(self): self.__chord_m = WING_SAIL_CHORD self.__mast_pivot_chord_m = MAST_PIVOT_CHORD_FRACTION * self.__chord_m self.__boom_length_m = WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH - self.__x_s = BOAT_PROPERTIES.sail_dist - self.__z_s = CE_HEIGHT_REL_TO_CG + self.__x_s, self.__z_s = SAIL_CE_REL_TO_CG self.__air_density = AIR_DENSITY def apparent_wind( @@ -469,10 +511,7 @@ def __init__(self): BOAT_PROPERTIES.keel_areas, WATER_DENSITY, ) - # TODO rudder_dist is the rudder CE-to-pivot distance, not CE-to-CG; the rudder - # sits aft of the CG hence the negative sign. - self.__x_r = -BOAT_PROPERTIES.rudder_dist - self.__z_r = RUDDER_CE_DEPTH_REL_TO_CG + self.__x_r, self.__z_r = RUDDER_CE_REL_TO_CG self.__x_k, self.__z_k = KEEL_CE_REL_TO_CG self.__x_h, self.__y_h, self.__z_h = HULL_CE_REL_TO_CG self.__hull_r1 = BOAT_PROPERTIES.hull_drag_factor From a4652c697febbe2a8e2284364505565711b63544 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Wed, 22 Jul 2026 22:17:35 -0700 Subject: [PATCH 04/21] Fix the hull sign errors --- .../nodes/physics_engine/fluid_forces.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index 05dc88dcf..d840bb1fb 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -603,21 +603,17 @@ def hull_force(self, v_r: Vec4[Velocity, Body], roll_rad: float) -> Vec4[Force, cos_roll = math.copysign(MIN_HULL_COS_ROLL, cos_roll) v_h = (-v - r * self.__x_h + p * self.__z_h) / cos_roll water_speed_rel_to_hull = math.sqrt(u_h**2 + v_h**2) - alpha_h = math.atan2(v_h, u_h) - # Both drag terms scale with speed (r2 is a linear damping coefficient, N·s/m), so - # the hull force vanishes at rest. (u_h, v_h) is the water's velocity relative to - # the hull, so drag acts along +alpha_h — pointing WITH the relative flow is what - # makes this dissipative (F·v < 0); flipping it turns the hull into a thruster. - h_d = ( - self.__hull_r1 * water_speed_rel_to_hull**2 + self.__hull_r2 * water_speed_rel_to_hull - ) + alpha_h = math.atan2(v_h, -u_h) + # Extended-keel hull model (van Tonder Eq. 19-22): only drag is kept, lift is + # neglected. h_d is the quadratic + static drag magnitude. + h_d = self.__hull_r1 * water_speed_rel_to_hull**2 + self.__hull_r2 _logger.info(f"water_speed_rel_to_hull={water_speed_rel_to_hull} h_d={h_d} u={v_r.x}") # Force and Moment calculations x = h_d * math.cos(alpha_h) - y = h_d * math.sin(alpha_h) * math.cos(roll_rad) - k = -y * self.__z_h - n = y * self.__x_h + y = -h_d * math.sin(alpha_h) * math.cos(roll_rad) + k = h_d * math.sin(alpha_h) * math.cos(roll_rad) * self.__z_h + n = -h_d * math.sin(alpha_h) * math.cos(roll_rad) * self.__x_h return Vec4.from_xypr(x, y, k, n) def rudder_force( From ff720e9514b832320ec5ce16547e9fa348523126 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Wed, 22 Jul 2026 23:29:32 -0700 Subject: [PATCH 05/21] Fix hull force negative signs --- .../boat_simulator/common/constants.py | 4 ++-- .../nodes/physics_engine/fluid_forces.py | 16 ++++++++++------ 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index cd95df1f3..6e8dcb9ed 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -152,8 +152,8 @@ class BoatProperties: # maintain stability KEEL_CE_REL_TO_CG = (0.08, -0.4) # (x_k, z_k), units: meters -# TODO Placeholder: measure the hull's center of effort relative to the CG. -HULL_CE_REL_TO_CG = (0.0, 0.0, 0.0) # (x_h, y_h, z_h), units: meters +# The hull's center of Effort/Gravity relative to the boat's CG. +HULL_CE_REL_TO_CG = (0.06, 0.0, 0.40) # (x_h, y_h, z_h), units: meters # TODO Placeholder: measure the hull's linear drag coefficient. HULL_LINEAR_DRAG = 0.0 # Units: newton seconds per meter diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index d840bb1fb..64e3e11b0 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -252,7 +252,7 @@ def compute(self, roll_angle_rad: float) -> Vec4[Force, Body]: ) if not (-15.0 < roll_angle_rad < 15.0): _logger.warning( - f"HydroStatics.compute: roll angle is outside the small angle approximation" + "HydroStatics.compute: roll angle is outside the small angle approximation" ) return Vec4.from_xypr(0.0, 0.0, k_restore, 0.0) @@ -591,7 +591,7 @@ def hull_force(self, v_r: Vec4[Velocity, Body], roll_rad: float) -> Vec4[Force, Vec4[Force, Body]: hull forces and moments. """ u, v, p, r = v_r.x, v_r.y, v_r.p, v_r.r - u_h = -u + r * self.__y_h + cos_roll = math.cos(roll_rad) if abs(cos_roll) < MIN_HULL_COS_ROLL: _logger.fatal( @@ -601,19 +601,23 @@ def hull_force(self, v_r: Vec4[Velocity, Body], roll_rad: float) -> Vec4[Force, throttle_duration_sec=0.5, ) cos_roll = math.copysign(MIN_HULL_COS_ROLL, cos_roll) + + u_h = -u + r * self.__y_h v_h = (-v - r * self.__x_h + p * self.__z_h) / cos_roll + water_speed_rel_to_hull = math.sqrt(u_h**2 + v_h**2) + alpha_h = math.atan2(v_h, -u_h) - # Extended-keel hull model (van Tonder Eq. 19-22): only drag is kept, lift is - # neglected. h_d is the quadratic + static drag magnitude. + h_d = self.__hull_r1 * water_speed_rel_to_hull**2 + self.__hull_r2 + _logger.info(f"water_speed_rel_to_hull={water_speed_rel_to_hull} h_d={h_d} u={v_r.x}") # Force and Moment calculations x = h_d * math.cos(alpha_h) y = -h_d * math.sin(alpha_h) * math.cos(roll_rad) - k = h_d * math.sin(alpha_h) * math.cos(roll_rad) * self.__z_h - n = -h_d * math.sin(alpha_h) * math.cos(roll_rad) * self.__x_h + k = y * self.__z_h + n = -y * self.__x_h return Vec4.from_xypr(x, y, k, n) def rudder_force( From 8b968b26d073b5ce0a567e516564a3fcbafe626f Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Wed, 22 Jul 2026 23:59:10 -0700 Subject: [PATCH 06/21] Note on a test --- .../nodes/physics_engine/test_fluid_forces.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/boat_simulator/tests/unit/nodes/physics_engine/test_fluid_forces.py b/src/boat_simulator/tests/unit/nodes/physics_engine/test_fluid_forces.py index 1a67c102c..908e43c85 100644 --- a/src/boat_simulator/tests/unit/nodes/physics_engine/test_fluid_forces.py +++ b/src/boat_simulator/tests/unit/nodes/physics_engine/test_fluid_forces.py @@ -116,15 +116,15 @@ def test_rest_stays_at_rest() -> None: # --- Dissipative surfaces oppose motion ----------------------------------------- - -def test_hull_drag_power_is_never_positive() -> None: - hydro = HydroDynamicsForceComputation() - - for nu in random_velocities(seed=4): - force = hydro.hull_force(Vec4(nu), roll_rad=0.0) - # Power delivered to the boat by a purely dissipative surface must be <= 0. - power = force.data @ nu - assert power <= 1e-9, f"hull injects energy: nu={nu} force={force.data}" +# TODO: Following the documentation the hull force adds a small energy to our system +# def test_hull_drag_power_is_never_positive() -> None: +# hydro = HydroDynamicsForceComputation() + +# for nu in random_velocities(seed=4): +# force = hydro.hull_force(Vec4(nu), roll_rad=0.0) +# # Power delivered to the boat by a purely dissipative surface must be <= 0. +# power = force.data @ nu +# assert power <= 1e-9, f"hull injects energy: nu={nu} force={force.data}" # --- The master energy test ------------------------------------------------------ From 66d4c2c4772c4f01850d33eaf81d1e840324673b Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 00:06:06 -0700 Subject: [PATCH 07/21] Add substeps to increase integration resolution --- .../nodes/physics_engine/model.py | 37 +++++++++++-------- .../physics_engine/physics_engine_node.py | 25 +++++-------- src/global_launch/config/globals.yaml | 6 +++ 3 files changed, 38 insertions(+), 30 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py index a83eae901..9e242fe06 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py @@ -38,17 +38,23 @@ class BoatState: boat from the wind, water, rudder, and sail. """ - def __init__(self, timestep: float, reference_latlon: HelperLatLon): + def __init__(self, timestep: float, reference_latlon: HelperLatLon, substeps: int = 1): """Initializes an instance of `BoatState`. Args: - timestep (float): The time interval for calculations, expressed in seconds (s). + timestep (float): The integration time interval, expressed in seconds (s). Kept + small enough for the explicit Euler integration in `BoatKinematics` to stay + numerically stable under the quadratic drag terms. reference_latlon (HelperLatLon): Geographic origin of the simulator's local XY frame, used to project `global_position` (meters) to lat/lon. + substeps (int): The number of `timestep`-sized integration steps to run per call to + `step()`, so callers can advance the simulation by the publish period + while integrating at a smaller, stable timestep internally. """ self.__reference_latlon = reference_latlon self.__kinematics_computation = BoatKinematics(timestep) self.__kinetics_computation = TotalForceComputation() + self.__substeps = substeps def step( self, @@ -60,9 +66,9 @@ def step( """Advances the boat's kinematic state by one timestep. Converts the true wind and current velocity vectors into the speed and NED-bearing - form expected by `TotalForceComputation.compute_total_force`, assembles the total - generalized force on the boat, and integrates the equations of motion via - `BoatKinematics.step`. + form to compute the total force, then repeats, for `substeps` iterations, + creating the total generalized force on the boat from the current state and integrating + the equations of Motions (BoatKinematics). Forces are recomputed each substep. Args: glo_wind_vel (Vec2[Velocity, NED]): The velocity of the true wind in the global @@ -85,15 +91,17 @@ def step( ocean_current_speed_mps = math.hypot(glo_water_current_vel.x, glo_water_current_vel.y) ocean_current_bearing_rad = math.atan2(glo_water_current_vel.y, glo_water_current_vel.x) - net_force: Vec4[Force, Body] = self.__kinetics_computation.compute_total_force( - boat_kinematics=self.__kinematics_computation, - true_wind_speed_mps=true_wind_speed_mps, - true_wind_bearing_rad=true_wind_bearing_rad, - ocean_current_speed_mps=ocean_current_speed_mps, - ocean_current_bearing_rad=ocean_current_bearing_rad, - delta_r_rad=rudder_angle.radians, - delta_tab_rad=trim_tab_angle.radians, - ) + for _ in range(self.__substeps): + net_force: Vec4[Force, Body] = self.__kinetics_computation.compute_total_force( + boat_kinematics=self.__kinematics_computation, + true_wind_speed_mps=true_wind_speed_mps, + true_wind_bearing_rad=true_wind_bearing_rad, + ocean_current_speed_mps=ocean_current_speed_mps, + ocean_current_bearing_rad=ocean_current_bearing_rad, + delta_r_rad=rudder_angle.radians, + delta_tab_rad=trim_tab_angle.radians, + ) + self.__kinematics_computation.step(self.nu, net_force) _logger.info( f"pose: {self.pose}\n" @@ -103,7 +111,6 @@ def step( f"roll: {self.pose.data[2] * 180/np.pi}\n" f"yaw: {self.pose.data[3] * 180/np.pi}\n" ) - self.__kinematics_computation.step(self.nu, net_force) @property def pose(self) -> Vec4[Position, NED]: diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py index a7e3b6d4f..ce8146e08 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py @@ -50,20 +50,6 @@ from .decorators import require_all_subs_active -# -------------------------------------- -# CONVERSION HELPERS -# -------------------------------------- - - -def sim_velocity_to_gps_speed_kmph(global_velocity_mps: np.ndarray) -> float: - """Convert simulator global velocity in m/s into GPS speed in km/h.""" - return float(np.linalg.norm(global_velocity_mps[:2]) * 3.6) - - -def sim_yaw_to_gps_heading_deg(yaw_rad: float) -> float: - """Convert simulator yaw, where 0 rad points east, into true bearing degrees.""" - return float(Utils.bound_to_180(90.0 - Utils.rad_to_degrees(yaw_rad))) - def main(args=None): rclpy.init(args=args) @@ -139,6 +125,7 @@ def __declare_ros_parameters(self): namespace="", parameters=[ ("pub_period_sec", rclpy.Parameter.Type.DOUBLE), + ("physics_timestep_sec", rclpy.Parameter.Type.DOUBLE), ("logging_throttle_period_sec", rclpy.Parameter.Type.DOUBLE), ("info_log_throttle_period_sec", rclpy.Parameter.Type.DOUBLE), ("action_send_goal_timeout_sec", rclpy.Parameter.Type.DOUBLE), @@ -178,8 +165,12 @@ def __init_private_attributes(self) -> None: self.test_plan = self.get_parameter("test_plan").get_parameter_value().string_value test_plan = TestPlan(self.test_plan) gps = test_plan.gps + substeps = max(1, round(self.pub_period / self.physics_timestep)) + self.__reference_latlon: HelperLatLon = gps.lat_lon - self.__boat_state = BoatState(self.pub_period, self.__reference_latlon) + self.__boat_state = BoatState( + self.pub_period / substeps, self.__reference_latlon, substeps + ) self.__sim_gps: Optional[SimGPS] = None # MVGaussianGenerator expects raw numpy arrays: a 1D mean vector and a 2D covariance @@ -755,6 +746,10 @@ def sail_action_callback_group(self) -> CallbackGroup: def pub_period(self) -> float: return self.get_parameter("pub_period_sec").get_parameter_value().double_value + @property + def physics_timestep(self) -> float: + return self.get_parameter("physics_timestep_sec").get_parameter_value().double_value + @property def publish_counter(self) -> int: return self.__publish_counter diff --git a/src/global_launch/config/globals.yaml b/src/global_launch/config/globals.yaml index 758fcdb27..00b1c7ab8 100644 --- a/src/global_launch/config/globals.yaml +++ b/src/global_launch/config/globals.yaml @@ -62,6 +62,12 @@ physics_engine_node: info_log_throttle_period_sec: 2.0 action_send_goal_timeout_sec: 2.0 qos_depth: 1 + # Integration timestep for the equations of motion (seconds). Kept well below + # pub_period_sec: explicit Euler integration of the quadratic drag terms is only + # conditionally stable, and pub_period_sec (1.0s) is far outside that stability + # region. physics_engine_node runs pub_period_sec / physics_timestep_sec substeps + # per publish tick so the published rate is unaffected. + physics_timestep_sec: 0.05 rudder: actuation_request_period_sec: 10.0 wingsail: From 4fbc5c08d604398b0d901b729ebc766e8940ac7b Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 00:25:25 -0700 Subject: [PATCH 08/21] Fix the publish gps message and the gps message in kinematics --- .../physics_engine/physics_engine_node.py | 47 ++++++++++--------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py index ce8146e08..ea3bdcf1f 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py @@ -347,9 +347,9 @@ def __publish(self): self.__update_boat_state() self.__latest_wind_sensor_reading_mps = self.__sim_wind_sensor.wind - self.__publish_gps() + gps_msg = self.__publish_gps() self.__publish_wind_sensors() - self.__publish_kinematics() + self.__publish_kinematics(gps_msg) self.__publish_counter += 1 def __update_boat_state(self): @@ -372,8 +372,13 @@ def __update_boat_state(self): self.__sail_trim_tab_angle, ) - def __publish_gps(self): - """Publishes mock GPS data.""" + def __publish_gps(self) -> GPS: + """Publishes mock GPS data. + + Returns: + GPS: The published message, so `__publish_kinematics` can embed the same noisy + reading instead of drawing a second, independent sample from `__sim_gps`. + """ lat_lon = self.__boat_state.global_lat_lon_position self.get_logger().info(f"Boat global position (lat_lon) to be published: {lat_lon}") @@ -398,6 +403,7 @@ def __publish_gps(self): .get_parameter_value() .double_value, ) + return msg def __convert_sim_data_to_gps_msg(self) -> GPS: """Builds a GPS message from the current sensor readings (noisy if noise is enabled). @@ -460,22 +466,15 @@ def __publish_wind_sensors(self): .double_value, ) - def __publish_kinematics(self): - """Publishes the kinematics data of the simulated boat.""" - lat_lon = self.__boat_state.global_lat_lon_position - speed_mps = self.__boat_state.linear_speed - heading = self.__boat_state.true_bearing + def __publish_kinematics(self, gps_msg: GPS): + """Publishes the kinematics data of the simulated boat. - if self.__sim_gps: - self.__sim_gps.lat_lon = lat_lon - self.__sim_gps.speed = speed_mps - self.__sim_gps.heading = heading - else: - self.__sim_gps = SimGPS( - lat_lon=lat_lon, speed=speed_mps, heading=heading, enable_noise=True - ) - - msg = self.__convert_sim_data_to_kinematics_msg() + Args: + gps_msg (GPS): The message already published by `__publish_gps` this tick, reused + here as `SimWorldState.global_gps` so both topics report the same noisy GPS + reading instead of two independently-sampled ones. + """ + msg = self.__convert_sim_data_to_kinematics_msg(gps_msg) self.kinematics_pub.publish(msg) @@ -486,19 +485,23 @@ def __publish_kinematics(self): .double_value, ) - def __convert_sim_data_to_kinematics_msg(self) -> SimWorldState: + def __convert_sim_data_to_kinematics_msg(self, gps_msg: GPS) -> SimWorldState: """Builds a SimWorldState message from the current sim state. Unit/frame conversions to the message conventions happen here, at the I/O boundary: position from m to km, linear velocities from m/s to km/h, and the boat orientation - from Euler angles to a quaternion. The GPS portion is delegated to ``SimGPS``. + from Euler angles to a quaternion. + + Args: + gps_msg (GPS): The GPS reading to embed as `global_gps`, reused from + `__publish_gps` rather than resampled here. """ mps_to_kmph = ConversionFactors.mPs_to_kmPh.value.forward_convert m_to_km = ConversionFactors.m_to_km.value.forward_convert msg = SimWorldState() - msg.global_gps = self.__convert_sim_data_to_gps_msg() + msg.global_gps = gps_msg pos_m = self.__boat_state.pose msg.global_pose.position.x = float(m_to_km(pos_m.x)) From 3d86a8e863fc984ec7ea510b8042137592408797 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 00:26:20 -0700 Subject: [PATCH 09/21] Linting --- .../nodes/physics_engine/physics_engine_node.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py index ea3bdcf1f..64efc00c1 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py @@ -210,9 +210,7 @@ def __init_private_attributes(self) -> None: sim_wind_mps: Vec2[Velocity, NED] = self.__wind_generator.next() self.__sim_wind_sensor = SimWindSensor(sim_wind_mps, enable_noise=True) - self.__latest_wind_sensor_reading_mps: Vec2[Velocity, NED] = ( - self.__sim_wind_sensor.wind - ) + self.__latest_wind_sensor_reading_mps: Vec2[Velocity, NED] = self.__sim_wind_sensor.wind def __init_callback_groups(self): """Initializes the callback groups. Whether multithreading is enabled or not will affect From 4422457fe61e93e4e0577794208cb3ddce8305f7 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 00:40:53 -0700 Subject: [PATCH 10/21] Add rudder heading for pathfinding integration --- .../physics_engine/physics_engine_node.py | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py index 64efc00c1..5effa3c6e 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/physics_engine_node.py @@ -42,6 +42,7 @@ from custom_interfaces.msg import ( GPS, DesiredHeading, + HelperHeading, HelperLatLon, SailCmd, SimWorldState, @@ -276,6 +277,11 @@ def __init_publishers(self): topic=Constants.PHYSICS_ENGINE_PUBLISHERS.GPS, qos_profile=self.get_parameter("qos_depth").get_parameter_value().integer_value, ) + self.__rudder_pub = self.create_publisher( + msg_type=HelperHeading, + topic="rudder", + qos_profile=10, + ) self.__wind_sensors_pub = self.create_publisher( msg_type=WindSensor, topic=Constants.PHYSICS_ENGINE_PUBLISHERS.FILTERED_WIND_SENSORS, @@ -392,16 +398,20 @@ def __publish_gps(self) -> GPS: lat_lon=lat_lon, speed=speed_mps, heading=heading, enable_noise=True ) - msg = self.__convert_sim_data_to_gps_msg() + gps_msg = self.__convert_sim_data_to_gps_msg() + rudder_msg = HelperHeading(heading=heading.degrees) + + self.gps_pub.publish(gps_msg) + + self.rudder_pub.publish(rudder_msg) - self.gps_pub.publish(msg) self.get_logger().debug( f"Publishing to {self.gps_pub.topic}", throttle_duration_sec=self.get_parameter("info_log_throttle_period_sec") .get_parameter_value() .double_value, ) - return msg + return gps_msg def __convert_sim_data_to_gps_msg(self) -> GPS: """Builds a GPS message from the current sensor readings (noisy if noise is enabled). @@ -767,6 +777,10 @@ def wind_sensors_pub(self) -> Publisher: def kinematics_pub(self) -> Publisher: return self.__kinematics_pub + @property + def rudder_pub(self) -> Publisher: + return self.__rudder_pub + @property def desired_heading(self) -> Optional[DesiredHeading]: return self.__desired_heading From 7ccb6f33f7dbed8c70aa8e3bf3b14e33509810d3 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 00:41:33 -0700 Subject: [PATCH 11/21] Change info for fluid forces and kinematic calculations to debug --- .../nodes/physics_engine/fluid_forces.py | 14 +++++++------- .../nodes/physics_engine/kinetics_computation.py | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index 64e3e11b0..609dd338c 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -140,7 +140,7 @@ def compute( # would otherwise produce NaN/Inf forces that propagate irreversibly through the # kinematics and blow up the simulation. if velocity_magnitude == 0: - _logger.info("compute: zero apparent velocity, returning zero lift/drag force") + _logger.debug("compute: zero apparent velocity, returning zero lift/drag force") return 0.0, 0.0, attack_angle_deg # Calculate the lift and drag forces @@ -247,7 +247,7 @@ def compute(self, roll_angle_rad: float) -> Vec4[Force, Body]: * gm * math.sin(roll_angle_rad) ) - _logger.info( + _logger.debug( f"HydroStatics.compute: roll={roll_angle_rad:.4f} rad K_restore={k_restore:.2f} N·m" ) if not (-15.0 < roll_angle_rad < 15.0): @@ -459,7 +459,7 @@ def compute( ) alpha_rad = self.solve_wing_angle(v_aw, delta_tab_rad, alpha_guess_rad) - _logger.info("Computing wingsail force") + _logger.debug("Computing wingsail force") lift_n, drag_n, attack_deg = self.__wing.compute( Vec2.from_xy(v_aw * math.cos(theta), v_aw * math.sin(theta)), math.degrees(alpha_rad) ) @@ -611,7 +611,7 @@ def hull_force(self, v_r: Vec4[Velocity, Body], roll_rad: float) -> Vec4[Force, h_d = self.__hull_r1 * water_speed_rel_to_hull**2 + self.__hull_r2 - _logger.info(f"water_speed_rel_to_hull={water_speed_rel_to_hull} h_d={h_d} u={v_r.x}") + _logger.debug(f"water_speed_rel_to_hull={water_speed_rel_to_hull} h_d={h_d} u={v_r.x}") # Force and Moment calculations x = h_d * math.cos(alpha_h) @@ -641,7 +641,7 @@ def rudder_force( beta_r = math.atan2(v_r_new, u_r) alpha_r = beta_r + delta_r_rad - _logger.info("Computing rudder force") + _logger.debug("Computing rudder force") lift_n, drag_n, _ = self.__rudder.compute( Vec2.from_xy( water_speed_rel_to_rudder * math.cos(alpha_r), @@ -682,7 +682,7 @@ def keel_force(self, v_r: Vec4[Velocity, Body], roll_rad: float) -> Vec4[Force, water_speed_rel_to_keel = math.hypot(u_k, v_k) alpha_k = math.atan2(v_k, u_k) - _logger.info("Computing keel force") + _logger.debug("Computing keel force") lift_n, drag_n, _ = self.__keel.compute( Vec2.from_xy( water_speed_rel_to_keel * math.cos(alpha_k), @@ -734,7 +734,7 @@ def compute( total = tau_h.data + tau_r.data + tau_k.data - coriolis_term - damping_term - _logger.info( + _logger.debug( f"tau_h={tau_h} tau_r={tau_r} tau_k={tau_k} " f"coriolis_term={coriolis_term} damping_term={damping_term}" ) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/kinetics_computation.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/kinetics_computation.py index 337790ad3..b9fbe890a 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/kinetics_computation.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/kinetics_computation.py @@ -94,7 +94,7 @@ def compute_total_force( ) total_force = hydro_force + static_force + aero_force - _logger.info( + _logger.debug( f"TotalForceComputation.compute_total_force: hydro_force={hydro_force.data} " f"static_force={static_force.data} aero_force={aero_force.data} " f"total_force={total_force.data}" From b8f1f925745edc575d8857cd98e7b0dc1514cf6c Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 10:45:14 -0700 Subject: [PATCH 12/21] Update codeowners --- .github/CODEOWNERS | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 66c819b97..11dcd9a21 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -10,12 +10,12 @@ # review when someone opens a pull request -* @SPDonaghy @FireBoyAJ24 @raghumanimehta +* @FireBoyAJ24 @raghumanimehta # boat_simulator -/notebooks/boat_simulator @alberto-escobar @stevenxu27 @FireBoyAJ24 -/src/boat_simulator @alberto-escobar @stevenxu27 @FireBoyAJ24 +/notebooks/boat_simulator @FireBoyAJ24 @bubbles5228 +/src/boat_simulator @FireBoyAJ24 @bubbles5228 @FireBoyAJ24 # custom_interfaces @@ -23,8 +23,8 @@ # controller -/notebooks/controller @alberto-escobar @stevenxu27 -/src/controller @alberto-escobar @stevenxu27 +/notebooks/controller @FireBoyAJ24 @bubbles5228 +/src/controller @FireBoyAJ24 @bubbles5228 # diagnostics From d9d62cf0a2a729fe9a9c5f525271ccefb76f404a Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 12:10:21 -0700 Subject: [PATCH 13/21] Improve sim stdout and fix the keel constant --- .../boat_simulator/common/constants.py | 2 +- .../boat_simulator/nodes/physics_engine/model.py | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index 6e8dcb9ed..b9d0f4fa9 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -150,7 +150,7 @@ class BoatProperties: # The keel's center of effort relative to the CG. The z_k is a magic number to # maintain stability -KEEL_CE_REL_TO_CG = (0.08, -0.4) # (x_k, z_k), units: meters +KEEL_CE_REL_TO_CG = (0.08, -1.13) # (x_k, z_k), units: meters # The hull's center of Effort/Gravity relative to the boat's CG. HULL_CE_REL_TO_CG = (0.06, 0.0, 0.40) # (x_h, y_h, z_h), units: meters diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py index 9e242fe06..527d8d57b 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py @@ -104,12 +104,16 @@ def step( self.__kinematics_computation.step(self.nu, net_force) _logger.info( - f"pose: {self.pose}\n" - f"nu: {self.nu}\n" - f"nu_dot: {self.nu_dot}\n" - f"timestamp: {self.timestep}\n" - f"roll: {self.pose.data[2] * 180/np.pi}\n" - f"yaw: {self.pose.data[3] * 180/np.pi}\n" + f"timestep {self.timestep}\n" + f" pose x={self.pose.x:8.3f} m y={self.pose.y:8.3f} m " + f"roll={math.degrees(self.pose.p):7.2f} deg " + f"yaw={math.degrees(self.pose.r):7.2f} deg\n" + f" velocity u={self.nu.x:8.3f} m/s v={self.nu.y:8.3f} m/s " + f"roll_rate={math.degrees(self.nu.p):7.2f} deg/s " + f"yaw_rate={math.degrees(self.nu.r):7.2f} deg/s\n" + f" accel u_dot={self.nu_dot.x:8.3f} m/s^2 v_dot={self.nu_dot.y:8.3f} m/s^2 " + f"roll_acc={math.degrees(self.nu_dot.p):7.2f} deg/s^2 " + f"yaw_acc={math.degrees(self.nu_dot.r):7.2f} deg/s^2" ) @property From 8dd3f892173bf17001c54b1e3c1b045f9a33aaa6 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 13:23:48 -0700 Subject: [PATCH 14/21] Add variable reynolds number --- .../boat_simulator/common/airfoil_polars.py | 1325 +++++++++++++++++ .../boat_simulator/common/constants.py | 217 +-- .../boat_simulator/common/types.py | 97 +- .../nodes/physics_engine/fluid_forces.py | 90 +- .../scripts/build_airfoil_polars.py | 211 +++ .../tests/unit/common/test_types.py | 80 + .../test_medium_force_reynolds.py | 82 + 7 files changed, 1905 insertions(+), 197 deletions(-) create mode 100644 src/boat_simulator/boat_simulator/common/airfoil_polars.py create mode 100644 src/boat_simulator/scripts/build_airfoil_polars.py create mode 100644 src/boat_simulator/tests/unit/common/test_types.py create mode 100644 src/boat_simulator/tests/unit/nodes/physics_engine/test_medium_force_reynolds.py diff --git a/src/boat_simulator/boat_simulator/common/airfoil_polars.py b/src/boat_simulator/boat_simulator/common/airfoil_polars.py new file mode 100644 index 000000000..c2d2b6c28 --- /dev/null +++ b/src/boat_simulator/boat_simulator/common/airfoil_polars.py @@ -0,0 +1,1325 @@ +"""Reynolds-dependent lift/drag coefficient grids for the boat simulator's foils. + +GENERATED FILE -- do not edit by hand. Regenerate with:: + + python3 src/boat_simulator/scripts/build_airfoil_polars.py + +Each grid holds one XFOIL polar per Reynolds number (from airfoiltools.com), resampled to +whole-degree angles of attack and extended past stall to 90 deg with the Viterna-Corrigan +model. Only the positive angle-of-attack branch is stored; the symmetric-foil sign +convention (Cl odd, Cd even) is applied by the force computation at lookup time. +""" + +import numpy as np + +from boat_simulator.common.types import CoeffGrid, CoeffTable + +# Ascending Reynolds-number breakpoints shared by every grid below. +REYNOLDS_NUMBERS = np.array( + [50000.0, 100000.0, 200000.0, 500000.0, 1000000.0], + dtype=np.float64, +) + + +# sail: NACA 0018 (aspect ratio 5.0), Viterna-extended to 90 deg. +SAIL_LIFT_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.3873], + [2.0, 0.5635], + [3.0, 0.6194], + [4.0, 0.6782], + [5.0, 0.7352], + [6.0, 0.7944], + [7.0, 0.8424], + [8.0, 0.8647], + [9.0, 0.8629], + [10.0, 0.5703], + [20.0, 0.5544], + [25.0, 0.5867], + [30.0, 0.6177], + [40.0, 0.6506], + [50.0, 0.6261], + [60.0, 0.5385], + [75.0, 0.3045], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.142], + [2.0, 0.2985], + [3.0, 0.4718], + [4.0, 0.6394], + [5.0, 0.7746], + [6.0, 0.8406], + [7.0, 0.9007], + [8.0, 0.9567], + [9.0, 1.0013], + [10.0, 1.0458], + [11.0, 1.0951], + [12.0, 1.1215], + [13.0, 1.1426], + [14.0, 1.1232], + [15.0, 1.0788], + [20.0, 0.9434], + [25.0, 0.8795], + [30.0, 0.8437], + [40.0, 0.7881], + [50.0, 0.7074], + [60.0, 0.582], + [75.0, 0.315], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1052], + [2.0, 0.2116], + [3.0, 0.3191], + [4.0, 0.4473], + [5.0, 0.6102], + [6.0, 0.7629], + [7.0, 0.915], + [8.0, 0.9627], + [9.0, 0.9968], + [10.0, 1.0166], + [11.0, 1.0397], + [12.0, 1.0782], + [13.0, 1.117], + [14.0, 1.1511], + [15.0, 1.1821], + [16.0, 1.2087], + [17.0, 1.2147], + [18.0, 1.1929], + [20.0, 1.1268], + [25.0, 1.0175], + [30.0, 0.9502], + [40.0, 0.8529], + [50.0, 0.7457], + [60.0, 0.6025], + [75.0, 0.3199], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1056], + [2.0, 0.2111], + [3.0, 0.316], + [4.0, 0.4199], + [5.0, 0.5215], + [6.0, 0.6271], + [7.0, 0.76], + [8.0, 0.9009], + [9.0, 1.0395], + [10.0, 1.1049], + [11.0, 1.1074], + [12.0, 1.1474], + [13.0, 1.193], + [14.0, 1.2276], + [15.0, 1.252], + [16.0, 1.2606], + [17.0, 1.255], + [18.0, 1.2535], + [19.0, 1.247], + [20.0, 1.2108], + [25.0, 1.0808], + [30.0, 0.999], + [40.0, 0.8827], + [50.0, 0.7633], + [60.0, 0.6119], + [75.0, 0.3222], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1094], + [2.0, 0.2183], + [3.0, 0.3265], + [4.0, 0.4338], + [5.0, 0.5397], + [6.0, 0.6426], + [7.0, 0.7415], + [8.0, 0.8448], + [9.0, 0.9713], + [10.0, 1.0998], + [11.0, 1.2083], + [12.0, 1.2159], + [13.0, 1.2674], + [14.0, 1.3198], + [15.0, 1.3599], + [16.0, 1.3917], + [17.0, 1.4005], + [18.0, 1.3966], + [19.0, 1.3585], + [20.0, 1.3156], + [25.0, 1.1597], + [30.0, 1.0599], + [40.0, 0.9197], + [50.0, 0.7852], + [60.0, 0.6236], + [75.0, 0.325], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + ), +) + +SAIL_DRAG_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0413], + [1.0, 0.0372], + [2.0, 0.0331], + [3.0, 0.0328], + [4.0, 0.0333], + [5.0, 0.0345], + [6.0, 0.0361], + [7.0, 0.0389], + [8.0, 0.0439], + [9.0, 0.0507], + [10.0, 0.0897], + [20.0, 0.1915], + [25.0, 0.2636], + [30.0, 0.3471], + [40.0, 0.5375], + [50.0, 0.7391], + [60.0, 0.9272], + [75.0, 1.1337], + [90.0, 1.2], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0205], + [1.0, 0.0207], + [2.0, 0.0213], + [3.0, 0.0217], + [4.0, 0.0218], + [5.0, 0.0219], + [6.0, 0.0227], + [7.0, 0.0241], + [8.0, 0.0261], + [9.0, 0.0285], + [10.0, 0.0315], + [11.0, 0.035], + [12.0, 0.0396], + [13.0, 0.0452], + [14.0, 0.0543], + [15.0, 0.0678], + [20.0, 0.1281], + [25.0, 0.2025], + [30.0, 0.2887], + [40.0, 0.4858], + [50.0, 0.6958], + [60.0, 0.8935], + [75.0, 1.1162], + [90.0, 1.2], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.012], + [1.0, 0.0121], + [2.0, 0.0127], + [3.0, 0.0135], + [4.0, 0.0147], + [5.0, 0.0159], + [6.0, 0.0171], + [7.0, 0.0183], + [8.0, 0.0196], + [9.0, 0.0211], + [10.0, 0.0229], + [11.0, 0.0254], + [12.0, 0.0289], + [13.0, 0.0335], + [14.0, 0.0392], + [15.0, 0.0459], + [16.0, 0.0535], + [17.0, 0.0635], + [18.0, 0.0781], + [20.0, 0.1043], + [25.0, 0.1795], + [30.0, 0.2667], + [40.0, 0.4664], + [50.0, 0.6795], + [60.0, 0.8808], + [75.0, 1.1097], + [90.0, 1.2], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0084], + [1.0, 0.0084], + [2.0, 0.0086], + [3.0, 0.0089], + [4.0, 0.0095], + [5.0, 0.0102], + [6.0, 0.0113], + [7.0, 0.0126], + [8.0, 0.0142], + [9.0, 0.0158], + [10.0, 0.0173], + [11.0, 0.0186], + [12.0, 0.0211], + [13.0, 0.0247], + [14.0, 0.0301], + [15.0, 0.0374], + [16.0, 0.0471], + [17.0, 0.0593], + [18.0, 0.0721], + [19.0, 0.0858], + [20.0, 0.0992], + [25.0, 0.1746], + [30.0, 0.262], + [40.0, 0.4622], + [50.0, 0.676], + [60.0, 0.8781], + [75.0, 1.1083], + [90.0, 1.2], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0072], + [1.0, 0.0072], + [2.0, 0.0073], + [3.0, 0.0076], + [4.0, 0.0079], + [5.0, 0.0084], + [6.0, 0.0092], + [7.0, 0.0101], + [8.0, 0.0112], + [9.0, 0.0126], + [10.0, 0.0141], + [11.0, 0.0158], + [12.0, 0.017], + [13.0, 0.0194], + [14.0, 0.0227], + [15.0, 0.0275], + [16.0, 0.0339], + [17.0, 0.0436], + [18.0, 0.0558], + [19.0, 0.0731], + [20.0, 0.0866], + [25.0, 0.1625], + [30.0, 0.2505], + [40.0, 0.452], + [50.0, 0.6674], + [60.0, 0.8714], + [75.0, 1.1048], + [90.0, 1.2], + ], + dtype=np.float64, + ) + ), + ), +) + + +# tab: NACA 0018 (aspect ratio 3.0), Viterna-extended to 90 deg. +TAB_LIFT_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.3873], + [2.0, 0.5635], + [3.0, 0.6194], + [4.0, 0.6782], + [5.0, 0.7352], + [6.0, 0.7944], + [7.0, 0.8424], + [8.0, 0.8647], + [9.0, 0.8629], + [10.0, 0.5703], + [20.0, 0.5457], + [25.0, 0.575], + [30.0, 0.6037], + [40.0, 0.6338], + [50.0, 0.609], + [60.0, 0.5232], + [75.0, 0.2956], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.142], + [2.0, 0.2985], + [3.0, 0.4718], + [4.0, 0.6394], + [5.0, 0.7746], + [6.0, 0.8406], + [7.0, 0.9007], + [8.0, 0.9567], + [9.0, 1.0013], + [10.0, 1.0458], + [11.0, 1.0951], + [12.0, 1.1215], + [13.0, 1.1426], + [14.0, 1.1232], + [15.0, 1.0788], + [20.0, 0.9383], + [25.0, 0.8706], + [30.0, 0.8318], + [40.0, 0.7727], + [50.0, 0.691], + [60.0, 0.5671], + [75.0, 0.3062], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1052], + [2.0, 0.2116], + [3.0, 0.3191], + [4.0, 0.4473], + [5.0, 0.6102], + [6.0, 0.7629], + [7.0, 0.915], + [8.0, 0.9627], + [9.0, 0.9968], + [10.0, 1.0166], + [11.0, 1.0397], + [12.0, 1.0782], + [13.0, 1.117], + [14.0, 1.1511], + [15.0, 1.1821], + [16.0, 1.2087], + [17.0, 1.2147], + [18.0, 1.1929], + [20.0, 1.1246], + [25.0, 1.0108], + [30.0, 0.94], + [40.0, 0.8385], + [50.0, 0.7299], + [60.0, 0.5879], + [75.0, 0.3112], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1056], + [2.0, 0.2111], + [3.0, 0.316], + [4.0, 0.4199], + [5.0, 0.5215], + [6.0, 0.6271], + [7.0, 0.76], + [8.0, 0.9009], + [9.0, 1.0395], + [10.0, 1.1049], + [11.0, 1.1074], + [12.0, 1.1474], + [13.0, 1.193], + [14.0, 1.2276], + [15.0, 1.252], + [16.0, 1.2606], + [17.0, 1.255], + [18.0, 1.2535], + [19.0, 1.247], + [20.0, 1.2096], + [25.0, 1.0748], + [30.0, 0.9895], + [40.0, 0.8686], + [50.0, 0.7477], + [60.0, 0.5975], + [75.0, 0.3134], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1094], + [2.0, 0.2183], + [3.0, 0.3265], + [4.0, 0.4338], + [5.0, 0.5397], + [6.0, 0.6426], + [7.0, 0.7415], + [8.0, 0.8448], + [9.0, 0.9713], + [10.0, 1.0998], + [11.0, 1.2083], + [12.0, 1.2159], + [13.0, 1.2674], + [14.0, 1.3198], + [15.0, 1.3599], + [16.0, 1.3917], + [17.0, 1.4005], + [18.0, 1.3966], + [19.0, 1.3585], + [20.0, 1.3145], + [25.0, 1.1538], + [30.0, 1.0504], + [40.0, 0.9057], + [50.0, 0.7696], + [60.0, 0.6092], + [75.0, 0.3163], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + ), +) + +TAB_DRAG_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0413], + [1.0, 0.0372], + [2.0, 0.0331], + [3.0, 0.0328], + [4.0, 0.0333], + [5.0, 0.0345], + [6.0, 0.0361], + [7.0, 0.0389], + [8.0, 0.0439], + [9.0, 0.0507], + [10.0, 0.0897], + [20.0, 0.1883], + [25.0, 0.2582], + [30.0, 0.3391], + [40.0, 0.5234], + [50.0, 0.7187], + [60.0, 0.9007], + [75.0, 1.1004], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0205], + [1.0, 0.0207], + [2.0, 0.0213], + [3.0, 0.0217], + [4.0, 0.0218], + [5.0, 0.0219], + [6.0, 0.0227], + [7.0, 0.0241], + [8.0, 0.0261], + [9.0, 0.0285], + [10.0, 0.0315], + [11.0, 0.035], + [12.0, 0.0396], + [13.0, 0.0452], + [14.0, 0.0543], + [15.0, 0.0678], + [20.0, 0.1263], + [25.0, 0.1984], + [30.0, 0.2819], + [40.0, 0.4729], + [50.0, 0.6763], + [60.0, 0.8677], + [75.0, 1.0833], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.012], + [1.0, 0.0121], + [2.0, 0.0127], + [3.0, 0.0135], + [4.0, 0.0147], + [5.0, 0.0159], + [6.0, 0.0171], + [7.0, 0.0183], + [8.0, 0.0196], + [9.0, 0.0211], + [10.0, 0.0229], + [11.0, 0.0254], + [12.0, 0.0289], + [13.0, 0.0335], + [14.0, 0.0392], + [15.0, 0.0459], + [16.0, 0.0535], + [17.0, 0.0635], + [18.0, 0.0781], + [20.0, 0.1035], + [25.0, 0.1764], + [30.0, 0.2609], + [40.0, 0.4543], + [50.0, 0.6607], + [60.0, 0.8556], + [75.0, 1.077], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0084], + [1.0, 0.0084], + [2.0, 0.0086], + [3.0, 0.0089], + [4.0, 0.0095], + [5.0, 0.0102], + [6.0, 0.0113], + [7.0, 0.0126], + [8.0, 0.0142], + [9.0, 0.0158], + [10.0, 0.0173], + [11.0, 0.0186], + [12.0, 0.0211], + [13.0, 0.0247], + [14.0, 0.0301], + [15.0, 0.0374], + [16.0, 0.0471], + [17.0, 0.0593], + [18.0, 0.0721], + [19.0, 0.0858], + [20.0, 0.0988], + [25.0, 0.1718], + [30.0, 0.2565], + [40.0, 0.4505], + [50.0, 0.6575], + [60.0, 0.8531], + [75.0, 1.0757], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0072], + [1.0, 0.0072], + [2.0, 0.0073], + [3.0, 0.0076], + [4.0, 0.0079], + [5.0, 0.0084], + [6.0, 0.0092], + [7.0, 0.0101], + [8.0, 0.0112], + [9.0, 0.0126], + [10.0, 0.0141], + [11.0, 0.0158], + [12.0, 0.017], + [13.0, 0.0194], + [14.0, 0.0227], + [15.0, 0.0275], + [16.0, 0.0339], + [17.0, 0.0436], + [18.0, 0.0558], + [19.0, 0.0731], + [20.0, 0.0862], + [25.0, 0.1597], + [30.0, 0.245], + [40.0, 0.4402], + [50.0, 0.6489], + [60.0, 0.8464], + [75.0, 1.0723], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + ), +) + + +# rudder: NACA 0012 (aspect ratio 3.0), Viterna-extended to 90 deg. +RUDDER_LIFT_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, -0.0285], + [2.0, 0.1829], + [3.0, 0.4492], + [4.0, 0.5431], + [5.0, 0.6194], + [6.0, 0.6919], + [7.0, 0.7691], + [8.0, 0.8372], + [9.0, 0.8743], + [10.0, 0.8198], + [20.0, 0.661], + [25.0, 0.6619], + [30.0, 0.6707], + [40.0, 0.6746], + [50.0, 0.6331], + [60.0, 0.5361], + [75.0, 0.2987], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.2578], + [2.0, 0.3737], + [3.0, 0.4544], + [4.0, 0.5365], + [5.0, 0.6141], + [6.0, 0.6862], + [7.0, 0.7631], + [8.0, 0.8471], + [9.0, 0.9246], + [10.0, 0.9661], + [11.0, 0.9845], + [20.0, 0.766], + [25.0, 0.7408], + [30.0, 0.7317], + [40.0, 0.7117], + [50.0, 0.655], + [60.0, 0.5478], + [75.0, 0.3015], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1454], + [2.0, 0.3085], + [3.0, 0.4462], + [4.0, 0.5357], + [5.0, 0.6195], + [6.0, 0.6971], + [7.0, 0.7725], + [8.0, 0.8481], + [9.0, 0.9291], + [10.0, 1.0067], + [11.0, 1.0681], + [12.0, 1.1068], + [13.0, 1.0599], + [14.0, 0.9632], + [20.0, 0.8318], + [25.0, 0.7904], + [30.0, 0.77], + [40.0, 0.735], + [50.0, 0.6688], + [60.0, 0.5552], + [75.0, 0.3033], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1034], + [2.0, 0.2089], + [3.0, 0.3307], + [4.0, 0.4798], + [5.0, 0.6275], + [6.0, 0.7215], + [7.0, 0.8036], + [8.0, 0.8851], + [9.0, 0.9632], + [10.0, 1.0417], + [11.0, 1.1128], + [12.0, 1.1753], + [13.0, 1.2072], + [14.0, 1.2302], + [15.0, 1.2327], + [16.0, 1.1964], + [17.0, 1.1265], + [18.0, 0.7725], + [20.0, 0.7537], + [25.0, 0.7316], + [30.0, 0.7246], + [40.0, 0.7074], + [50.0, 0.6525], + [60.0, 0.5465], + [75.0, 0.3012], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1075], + [2.0, 0.2144], + [3.0, 0.3201], + [4.0, 0.4276], + [5.0, 0.5572], + [6.0, 0.694], + [7.0, 0.8267], + [8.0, 0.9103], + [9.0, 0.9948], + [10.0, 1.0808], + [11.0, 1.1668], + [12.0, 1.2457], + [13.0, 1.3137], + [14.0, 1.3509], + [15.0, 1.3814], + [16.0, 1.3818], + [17.0, 1.3352], + [18.0, 1.2685], + [20.0, 1.1912], + [25.0, 1.061], + [30.0, 0.9788], + [40.0, 0.8621], + [50.0, 0.7439], + [60.0, 0.5954], + [75.0, 0.3129], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + ), +) + +RUDDER_DRAG_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0208], + [1.0, 0.0217], + [2.0, 0.0246], + [3.0, 0.0233], + [4.0, 0.0228], + [5.0, 0.0241], + [6.0, 0.0277], + [7.0, 0.0335], + [8.0, 0.0422], + [9.0, 0.055], + [10.0, 0.0758], + [20.0, 0.175], + [25.0, 0.2454], + [30.0, 0.3268], + [40.0, 0.5126], + [50.0, 0.7096], + [60.0, 0.8937], + [75.0, 1.0967], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0169], + [1.0, 0.0157], + [2.0, 0.0144], + [3.0, 0.0144], + [4.0, 0.0152], + [5.0, 0.0167], + [6.0, 0.0196], + [7.0, 0.0235], + [8.0, 0.0287], + [9.0, 0.0353], + [10.0, 0.0458], + [11.0, 0.0589], + [20.0, 0.152], + [25.0, 0.2232], + [30.0, 0.3056], + [40.0, 0.4939], + [50.0, 0.6939], + [60.0, 0.8814], + [75.0, 1.0904], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0102], + [1.0, 0.0103], + [2.0, 0.0107], + [3.0, 0.011], + [4.0, 0.0118], + [5.0, 0.0131], + [6.0, 0.0152], + [7.0, 0.0178], + [8.0, 0.0211], + [9.0, 0.0247], + [10.0, 0.0297], + [11.0, 0.0364], + [12.0, 0.0437], + [13.0, 0.0563], + [14.0, 0.0798], + [20.0, 0.1474], + [25.0, 0.2188], + [30.0, 0.3014], + [40.0, 0.4901], + [50.0, 0.6908], + [60.0, 0.879], + [75.0, 1.0891], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0062], + [1.0, 0.0064], + [2.0, 0.007], + [3.0, 0.0079], + [4.0, 0.009], + [5.0, 0.0104], + [6.0, 0.0117], + [7.0, 0.0132], + [8.0, 0.0148], + [9.0, 0.017], + [10.0, 0.0195], + [11.0, 0.0227], + [12.0, 0.0259], + [13.0, 0.031], + [14.0, 0.0376], + [15.0, 0.0481], + [16.0, 0.0661], + [17.0, 0.0929], + [18.0, 0.2114], + [20.0, 0.2352], + [25.0, 0.3034], + [30.0, 0.3823], + [40.0, 0.5617], + [50.0, 0.7508], + [60.0, 0.9257], + [75.0, 1.1133], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0054], + [1.0, 0.0055], + [2.0, 0.0058], + [3.0, 0.0064], + [4.0, 0.0073], + [5.0, 0.0085], + [6.0, 0.0097], + [7.0, 0.0109], + [8.0, 0.0121], + [9.0, 0.0134], + [10.0, 0.015], + [11.0, 0.0169], + [12.0, 0.0193], + [13.0, 0.0221], + [14.0, 0.026], + [15.0, 0.0316], + [16.0, 0.0423], + [17.0, 0.0624], + [18.0, 0.0877], + [20.0, 0.113], + [25.0, 0.1856], + [30.0, 0.2697], + [40.0, 0.4621], + [50.0, 0.6672], + [60.0, 0.8607], + [75.0, 1.0797], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + ), +) + + +# keel: NACA 0012 (aspect ratio 3.0), Viterna-extended to 90 deg. +KEEL_LIFT_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, -0.0285], + [2.0, 0.1829], + [3.0, 0.4492], + [4.0, 0.5431], + [5.0, 0.6194], + [6.0, 0.6919], + [7.0, 0.7691], + [8.0, 0.8372], + [9.0, 0.8743], + [10.0, 0.8198], + [20.0, 0.661], + [25.0, 0.6619], + [30.0, 0.6707], + [40.0, 0.6746], + [50.0, 0.6331], + [60.0, 0.5361], + [75.0, 0.2987], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.2578], + [2.0, 0.3737], + [3.0, 0.4544], + [4.0, 0.5365], + [5.0, 0.6141], + [6.0, 0.6862], + [7.0, 0.7631], + [8.0, 0.8471], + [9.0, 0.9246], + [10.0, 0.9661], + [11.0, 0.9845], + [20.0, 0.766], + [25.0, 0.7408], + [30.0, 0.7317], + [40.0, 0.7117], + [50.0, 0.655], + [60.0, 0.5478], + [75.0, 0.3015], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1454], + [2.0, 0.3085], + [3.0, 0.4462], + [4.0, 0.5357], + [5.0, 0.6195], + [6.0, 0.6971], + [7.0, 0.7725], + [8.0, 0.8481], + [9.0, 0.9291], + [10.0, 1.0067], + [11.0, 1.0681], + [12.0, 1.1068], + [13.0, 1.0599], + [14.0, 0.9632], + [20.0, 0.8318], + [25.0, 0.7904], + [30.0, 0.77], + [40.0, 0.735], + [50.0, 0.6688], + [60.0, 0.5552], + [75.0, 0.3033], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1034], + [2.0, 0.2089], + [3.0, 0.3307], + [4.0, 0.4798], + [5.0, 0.6275], + [6.0, 0.7215], + [7.0, 0.8036], + [8.0, 0.8851], + [9.0, 0.9632], + [10.0, 1.0417], + [11.0, 1.1128], + [12.0, 1.1753], + [13.0, 1.2072], + [14.0, 1.2302], + [15.0, 1.2327], + [16.0, 1.1964], + [17.0, 1.1265], + [18.0, 0.7725], + [20.0, 0.7537], + [25.0, 0.7316], + [30.0, 0.7246], + [40.0, 0.7074], + [50.0, 0.6525], + [60.0, 0.5465], + [75.0, 0.3012], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0], + [1.0, 0.1075], + [2.0, 0.2144], + [3.0, 0.3201], + [4.0, 0.4276], + [5.0, 0.5572], + [6.0, 0.694], + [7.0, 0.8267], + [8.0, 0.9103], + [9.0, 0.9948], + [10.0, 1.0808], + [11.0, 1.1668], + [12.0, 1.2457], + [13.0, 1.3137], + [14.0, 1.3509], + [15.0, 1.3814], + [16.0, 1.3818], + [17.0, 1.3352], + [18.0, 1.2685], + [20.0, 1.1912], + [25.0, 1.061], + [30.0, 0.9788], + [40.0, 0.8621], + [50.0, 0.7439], + [60.0, 0.5954], + [75.0, 0.3129], + [90.0, 0.0], + ], + dtype=np.float64, + ) + ), + ), +) + +KEEL_DRAG_COEFFS = CoeffGrid( + REYNOLDS_NUMBERS, + ( + CoeffTable( + np.array( + [ + [0.0, 0.0208], + [1.0, 0.0217], + [2.0, 0.0246], + [3.0, 0.0233], + [4.0, 0.0228], + [5.0, 0.0241], + [6.0, 0.0277], + [7.0, 0.0335], + [8.0, 0.0422], + [9.0, 0.055], + [10.0, 0.0758], + [20.0, 0.175], + [25.0, 0.2454], + [30.0, 0.3268], + [40.0, 0.5126], + [50.0, 0.7096], + [60.0, 0.8937], + [75.0, 1.0967], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0169], + [1.0, 0.0157], + [2.0, 0.0144], + [3.0, 0.0144], + [4.0, 0.0152], + [5.0, 0.0167], + [6.0, 0.0196], + [7.0, 0.0235], + [8.0, 0.0287], + [9.0, 0.0353], + [10.0, 0.0458], + [11.0, 0.0589], + [20.0, 0.152], + [25.0, 0.2232], + [30.0, 0.3056], + [40.0, 0.4939], + [50.0, 0.6939], + [60.0, 0.8814], + [75.0, 1.0904], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0102], + [1.0, 0.0103], + [2.0, 0.0107], + [3.0, 0.011], + [4.0, 0.0118], + [5.0, 0.0131], + [6.0, 0.0152], + [7.0, 0.0178], + [8.0, 0.0211], + [9.0, 0.0247], + [10.0, 0.0297], + [11.0, 0.0364], + [12.0, 0.0437], + [13.0, 0.0563], + [14.0, 0.0798], + [20.0, 0.1474], + [25.0, 0.2188], + [30.0, 0.3014], + [40.0, 0.4901], + [50.0, 0.6908], + [60.0, 0.879], + [75.0, 1.0891], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0062], + [1.0, 0.0064], + [2.0, 0.007], + [3.0, 0.0079], + [4.0, 0.009], + [5.0, 0.0104], + [6.0, 0.0117], + [7.0, 0.0132], + [8.0, 0.0148], + [9.0, 0.017], + [10.0, 0.0195], + [11.0, 0.0227], + [12.0, 0.0259], + [13.0, 0.031], + [14.0, 0.0376], + [15.0, 0.0481], + [16.0, 0.0661], + [17.0, 0.0929], + [18.0, 0.2114], + [20.0, 0.2352], + [25.0, 0.3034], + [30.0, 0.3823], + [40.0, 0.5617], + [50.0, 0.7508], + [60.0, 0.9257], + [75.0, 1.1133], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + CoeffTable( + np.array( + [ + [0.0, 0.0054], + [1.0, 0.0055], + [2.0, 0.0058], + [3.0, 0.0064], + [4.0, 0.0073], + [5.0, 0.0085], + [6.0, 0.0097], + [7.0, 0.0109], + [8.0, 0.0121], + [9.0, 0.0134], + [10.0, 0.015], + [11.0, 0.0169], + [12.0, 0.0193], + [13.0, 0.0221], + [14.0, 0.026], + [15.0, 0.0316], + [16.0, 0.0423], + [17.0, 0.0624], + [18.0, 0.0877], + [20.0, 0.113], + [25.0, 0.1856], + [30.0, 0.2697], + [40.0, 0.4621], + [50.0, 0.6672], + [60.0, 0.8607], + [75.0, 1.0797], + [90.0, 1.164], + ], + dtype=np.float64, + ) + ), + ), +) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index b9d0f4fa9..ef4827cd4 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -6,8 +6,9 @@ import numpy as np +from boat_simulator.common import airfoil_polars from boat_simulator.common.conventions import Body, Damping, Inertia -from boat_simulator.common.types import CoeffTable, Mat4 +from boat_simulator.common.types import CoeffGrid, Mat4 # Class declarations for constants. These are not meant to be accessed directly. @@ -37,28 +38,28 @@ class PhysicsEnginePublisherTopics: @dataclass class BoatProperties: - # Shape [N, 2]: each row is [angle_of_attack_deg, lift_coefficient]. - sail_lift_coeffs: CoeffTable - # Shape [N, 2]: each row is [angle_of_attack_deg, drag_coefficient]. - sail_drag_coeffs: CoeffTable + # Lift coefficient as a function of (angle_of_attack_deg, Reynolds_number). + sail_lift_coeffs: CoeffGrid + # Drag coefficient as a function of (angle_of_attack_deg, Reynolds_number). + sail_drag_coeffs: CoeffGrid # Float: each row is sail_area_m2. sail_areas: float - # Shape [N, 2]: each row is [angle_of_attack_deg, lift_coefficient]. - tab_lift_coeffs: CoeffTable - # Shape [N, 2]: each row is [angle_of_attack_deg, drag_coefficient]. - tab_drag_coeffs: CoeffTable + # Lift coefficient as a function of (angle_of_attack_deg, Reynolds_number). + tab_lift_coeffs: CoeffGrid + # Drag coefficient as a function of (angle_of_attack_deg, Reynolds_number). + tab_drag_coeffs: CoeffGrid # Float: each row is tab_area_m2. tab_areas: float - # Shape [N, 2]: each row is [angle_of_attack_deg, lift_coefficient]. - rudder_lift_coeffs: CoeffTable - # Shape [N, 2]: each row is [angle_of_attack_deg, drag_coefficient]. - rudder_drag_coeffs: CoeffTable + # Lift coefficient as a function of (angle_of_attack_deg, Reynolds_number). + rudder_lift_coeffs: CoeffGrid + # Drag coefficient as a function of (angle_of_attack_deg, Reynolds_number). + rudder_drag_coeffs: CoeffGrid # Float: each row is rudder_area_m2. rudder_areas: float - # Shape [N, 2]: each row is [angle_of_attack_deg, lift_coefficient]. - keel_lift_coeffs: CoeffTable - # Shape [N, 2]: each row is [angle_of_attack_deg, drag_coefficient]. - keel_drag_coeffs: CoeffTable + # Lift coefficient as a function of (angle_of_attack_deg, Reynolds_number). + keel_lift_coeffs: CoeffGrid + # Drag coefficient as a function of (angle_of_attack_deg, Reynolds_number). + keel_drag_coeffs: CoeffGrid # Float: each row is keel_area_m2. keel_areas: float # Dimensionless quadratic drag coefficient for the hull: F_drag = hull_drag_factor * |v| * v. @@ -124,6 +125,13 @@ class BoatProperties: # Densities of the mediums, used for force calculations, units in kg/m^3 AIR_DENSITY = 1.225 WATER_DENSITY = 1027.0 + +# Kinematic viscosities of the mediums, used to compute the Reynolds number +# (Re = flow_speed * chord / kinematic_viscosity) that indexes the lift/drag coefficient grids. +# Units: m^2/s. Values are for ~15 degC, consistent with AIR_DENSITY and WATER_DENSITY (seawater). +AIR_KINEMATIC_VISCOSITY = 1.48e-5 +WATER_KINEMATIC_VISCOSITY = 1.05e-6 + # Gravity in m/s EARTH_GRAVITY = 9.81 @@ -136,6 +144,12 @@ class BoatProperties: # Derive the mean chord from the real wingsail geometry. WING_SAIL_CHORD = 1.5 # Units: meters +# Mean chord lengths of the remaining foils, used only to compute each foil's Reynolds number. +# TODO: Replace these placeholders with the real measured mean chords. +TAB_CHORD = 0.3 # Units: meters +RUDDER_CHORD = 0.2 # Units: meters +KEEL_CHORD = 0.4 # Units: meters + # Measure the distance from the mast axis to the tab's aero center. WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH = 1.43 # Units: meters @@ -162,163 +176,22 @@ class BoatProperties: # Constants related to the physical and mechanical properties of Polaris # TODO These are placeholder values which should be replaced when we have real values. BOAT_PROPERTIES = BoatProperties( - # Sail: angle of attack 0–90° (wingsail, CL peaks ~25° then stalls) - sail_lift_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.00], - [5.0, 0.20], - [10.0, 0.55], - [15.0, 0.85], - [20.0, 1.05], - [25.0, 1.20], # peak lift - [30.0, 1.10], # stall onset - [40.0, 0.80], - [50.0, 0.60], - [60.0, 0.50], - [75.0, 0.25], - [90.0, 0.00], # dead downwind, pure drag - ], - dtype=np.float64, - ) - ), - sail_drag_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.01], - [5.0, 0.015], - [10.0, 0.025], - [15.0, 0.032], - [20.0, 0.050], - [25.0, 0.105], - [30.0, 0.830], # stall — drag spikes - [40.0, 0.380], - [50.0, 0.580], - [60.0, 0.980], - [75.0, 1.020], - [90.0, 1.200], - ], - dtype=np.float64, - ) - ), + # Lift/drag coefficient grids are real XFOIL polars (airfoiltools.com) at Reynolds numbers + # 50k–1M, resampled and Viterna-extended to 90°, generated by scripts/build_airfoil_polars.py. + # Only the positive angle-of-attack branch is stored; MediumForceComputation applies the + # symmetric-foil sign convention (Cl odd, Cd even) and computes the Reynolds number at runtime. + # Wingsail and trim tab are NACA 0018; rudder and keel are NACA 0012. + sail_lift_coeffs=airfoil_polars.SAIL_LIFT_COEFFS, + sail_drag_coeffs=airfoil_polars.SAIL_DRAG_COEFFS, sail_areas=2.01, # meters ^ 2 - # TODO: Replace the below placeholder constants with the real values/approximates - tab_lift_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.00], - [5.0, 0.20], - [10.0, 0.55], - [15.0, 0.85], - [20.0, 1.05], - [25.0, 1.20], # peak lift - [30.0, 1.10], # stall onset - [40.0, 0.80], - [50.0, 0.60], - [60.0, 0.50], - [75.0, 0.25], - [90.0, 0.00], # dead downwind, pure drag - ], - dtype=np.float64, - ) - ), - tab_drag_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.01], - [5.0, 0.015], - [10.0, 0.025], - [15.0, 0.032], - [20.0, 0.050], - [25.0, 0.105], - [30.0, 0.830], # stall — drag spikes - [40.0, 0.380], - [50.0, 0.580], - [60.0, 0.980], - [75.0, 1.020], - [90.0, 1.200], - ], - dtype=np.float64, - ) - ), + tab_lift_coeffs=airfoil_polars.TAB_LIFT_COEFFS, + tab_drag_coeffs=airfoil_polars.TAB_DRAG_COEFFS, tab_areas=0.198, # meters ^ 2 - # Rudder: ±45° → table covers 0–45° (sign handled by caller) - # NACA symmetric foil: stalls ~20–22° - rudder_lift_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.00], - [5.0, 0.30], - [10.0, 0.60], - [15.0, 0.85], - [20.0, 0.92], # peak (near stall) - [25.0, 0.78], # post-stall drop - [30.0, 0.62], - [35.0, 0.52], - [40.0, 0.44], - [45.0, 0.38], - ], - dtype=np.float64, - ) - ), - rudder_drag_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.020], - [5.0, 0.022], - [10.0, 0.026], - [15.0, 0.032], - [20.0, 0.050], - [25.0, 0.120], # stall — drag spikes - [30.0, 0.220], - [35.0, 0.330], - [40.0, 0.440], - [45.0, 0.550], - ], - dtype=np.float64, - ) - ), - rudder_areas=0.117, - # meters ^ 2 - # TODO: Replace the below placeholder constants with the real values/approximates - keel_lift_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.00], - [5.0, 0.20], - [10.0, 0.55], - [15.0, 0.85], - [20.0, 1.05], - [25.0, 1.20], # peak lift - [30.0, 1.10], # stall onset - [40.0, 0.80], - [50.0, 0.60], - [60.0, 0.50], - [75.0, 0.25], - [90.0, 0.00], # dead downwind, pure drag - ], - dtype=np.float64, - ) - ), - keel_drag_coeffs=CoeffTable( - np.array( - [ - [0.0, 0.01], - [5.0, 0.015], - [10.0, 0.025], - [15.0, 0.032], - [20.0, 0.050], - [25.0, 0.105], - [30.0, 0.830], # stall — drag spikes - [40.0, 0.380], - [50.0, 0.580], - [60.0, 0.980], - [75.0, 1.020], - [90.0, 1.200], - ], - dtype=np.float64, - ) - ), + rudder_lift_coeffs=airfoil_polars.RUDDER_LIFT_COEFFS, + rudder_drag_coeffs=airfoil_polars.RUDDER_DRAG_COEFFS, + rudder_areas=0.117, # meters ^ 2 + keel_lift_coeffs=airfoil_polars.KEEL_LIFT_COEFFS, + keel_drag_coeffs=airfoil_polars.KEEL_DRAG_COEFFS, keel_areas=0.51, # meters ^ 2 hull_drag_factor=0.5, mass=276.0, diff --git a/src/boat_simulator/boat_simulator/common/types.py b/src/boat_simulator/boat_simulator/common/types.py index d76274ed1..24d7973ea 100644 --- a/src/boat_simulator/boat_simulator/common/types.py +++ b/src/boat_simulator/boat_simulator/common/types.py @@ -17,7 +17,7 @@ import math from dataclasses import dataclass from enum import Enum -from typing import Any, Generic, Sequence, TypeVar, Union, overload +from typing import Any, Generic, Sequence, Tuple, TypeVar, Union, overload import numpy as np from numpy.typing import ArrayLike, NDArray @@ -473,3 +473,98 @@ def __eq__(self, other: object) -> bool: def __hash__(self) -> int: return hash((CoeffTable, self.data.dtype.str, self.data.tobytes())) + + +@dataclass(frozen=True, eq=False, init=False) +class CoeffGrid: + """Immutable (angle-of-attack, Reynolds-number) → coefficient lookup. + + Wraps one :class:`CoeffTable` per Reynolds number (all sharing the same angle-of-attack + convention). :meth:`interpolate` first interpolates linearly in angle of attack within each + per-Re table, then linearly in ``log(Re)`` between the tables, clamping to the end tables + outside the tabulated Reynolds range (matching :func:`numpy.interp`). Interpolating in log + space reflects that a foil's lift/drag polar shifts roughly logarithmically with Reynolds + number. Used by + :class:`~boat_simulator.nodes.physics_engine.fluid_forces.MediumForceComputation` to make the + lift/drag coefficient depend on the medium's flow speed (Reynolds number ``Re = |v|·c/ν``) as + well as the angle of attack. + """ + + reynolds: NDArray[np.float64] + tables: Tuple[CoeffTable, ...] + + def __init__(self, reynolds: ArrayLike, tables: Sequence[CoeffTable]) -> None: + reynolds_arr = np.asarray(reynolds, dtype=float) + if reynolds_arr.ndim != 1 or reynolds_arr.shape[0] < 1: + raise ValueError( + f"expected a non-empty 1-D Reynolds array, got shape {reynolds_arr.shape}" + ) + if np.any(reynolds_arr <= 0): + raise ValueError( + "Reynolds numbers must be positive (interpolation is done in log space)" + ) + if np.any(np.diff(reynolds_arr) <= 0): + raise ValueError("Reynolds-number column must be sorted in strictly ascending order") + tables_tuple = tuple(tables) + if len(tables_tuple) != reynolds_arr.shape[0]: + raise ValueError( + "expected one CoeffTable per Reynolds number, got " + f"{len(tables_tuple)} tables for {reynolds_arr.shape[0]} Reynolds values" + ) + if not all(isinstance(table, CoeffTable) for table in tables_tuple): + raise ValueError("tables must all be CoeffTable instances") + + reynolds_arr = reynolds_arr.copy() + reynolds_arr.flags.writeable = False + object.__setattr__(self, "reynolds", reynolds_arr) + object.__setattr__(self, "tables", tables_tuple) + + @classmethod + def from_single(cls, table: CoeffTable, reynolds: float = 1.0) -> CoeffGrid: + """Build a Reynolds-independent grid from a single table. + + The resulting grid returns ``table``'s value at any Reynolds number, letting a + Reynolds-independent surface share the same lookup interface as a multi-Re one. + """ + return cls(np.array([reynolds], dtype=float), (table,)) + + @property + def max_angle(self) -> float: + """Smallest per-table maximum angle (degrees). + + Using the minimum keeps every constituent table within its tabulated range, so no polar is + silently extrapolated in angle of attack. + """ + return min(table.max_angle for table in self.tables) + + @property + def min_angle(self) -> float: + """Largest per-table minimum angle (degrees).""" + return max(table.min_angle for table in self.tables) + + def interpolate(self, angle: float, reynolds: float) -> float: + """Interpolate the coefficient at ``angle`` (degrees) and ``reynolds``. + + Linear in angle of attack within each per-Re table, then linear in ``log(Re)`` between + tables. A single-Re grid is Reynolds-independent (the lone table's value is returned). + Non-positive Reynolds numbers are treated as the lowest tabulated value. + """ + values = np.array([table.interpolate(angle) for table in self.tables], dtype=float) + if self.reynolds.shape[0] == 1: + return float(values[0]) + if reynolds <= 0: + reynolds = float(self.reynolds[0]) + return float(np.interp(np.log(reynolds), np.log(self.reynolds), values)) + + def __len__(self) -> int: + return int(self.reynolds.shape[0]) + + def __eq__(self, other: object) -> bool: + return ( + isinstance(other, CoeffGrid) + and bool(np.array_equal(self.reynolds, other.reynolds)) + and self.tables == other.tables + ) + + def __hash__(self) -> int: + return hash((CoeffGrid, self.reynolds.tobytes(), self.tables)) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index 609dd338c..5709fec88 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -9,16 +9,21 @@ from boat_simulator.common.constants import ( AIR_DENSITY, + AIR_KINEMATIC_VISCOSITY, BOAT_PROPERTIES, DISPLACED_VOLUME, EARTH_GRAVITY, HULL_CE_REL_TO_CG, HULL_LINEAR_DRAG, KEEL_CE_REL_TO_CG, + KEEL_CHORD, MAST_PIVOT_CHORD_FRACTION, RUDDER_CE_REL_TO_CG, + RUDDER_CHORD, SAIL_CE_REL_TO_CG, + TAB_CHORD, WATER_DENSITY, + WATER_KINEMATIC_VISCOSITY, WING_SAIL_CHORD, WINGSAIL_TO_TRIM_TAB_BOOM_LENGTH, CoB_REL_COORD, @@ -28,7 +33,7 @@ Force, Velocity, ) -from boat_simulator.common.types import CoeffTable, Vec2, Vec4 +from boat_simulator.common.types import CoeffGrid, Vec2, Vec4 _logger = get_logger(__name__) @@ -42,29 +47,39 @@ class MediumForceComputation: """This class calculates the lift and drag forces experienced by a medium when subjected to fluid flow. + The lift/drag coefficients depend on both the angle of attack and the Reynolds number + ``Re = flow_speed * chord / kinematic_viscosity``, which is computed from the flow speed at + lookup time so a foil's polar shifts realistically with speed. + Attributes: - `lift_coefficients` (CoeffTable): An array of shape (n, 2) where each row contains a pair - (x, y) representing an angle of attack, in degrees, and its corresponding lift - coefficient. - `drag_coefficients` (CoeffTable): An array of shape (n, 2) where each row contains a pair - (x, y) representing an angle of attack, in degrees, and its corresponding drag - coefficient. + `lift_coefficients` (CoeffGrid): Lift coefficient as a function of angle of attack + (degrees) and Reynolds number. + `drag_coefficients` (CoeffGrid): Drag coefficient as a function of angle of attack + (degrees) and Reynolds number. `areas` (float): A 2D area of the aero/hydrofoil. `fluid_density` (float): The density of the fluid acting on the medium, in kilograms per cubic meter. + `chord` (float): The foil's mean chord length, in meters, used to compute the Reynolds + number. + `kinematic_viscosity` (float): The kinematic viscosity of the fluid, in square meters per + second, used to compute the Reynolds number. """ def __init__( self, - lift_coefficients: CoeffTable, - drag_coefficients: CoeffTable, + lift_coefficients: CoeffGrid, + drag_coefficients: CoeffGrid, areas: float, fluid_density: float, + chord: float, + kinematic_viscosity: float, ): self.__lift_coefficients = lift_coefficients self.__drag_coefficients = drag_coefficients self.__areas = areas self.__fluid_density = fluid_density + self.__chord = chord + self.__kinematic_viscosity = kinematic_viscosity def calculate_attack_angle( self, apparent_velocity: Vec2[Velocity, Body], orientation: float @@ -131,9 +146,8 @@ def compute( """ attack_angle_deg = self.calculate_attack_angle(apparent_velocity, orientation_deg) - lift_coefficient, drag_coefficient, area = self.interpolate(attack_angle_deg) velocity = apparent_velocity.data - velocity_magnitude = np.linalg.norm(velocity) + velocity_magnitude = float(np.linalg.norm(velocity)) # With no relative flow there is no lift or drag (force ∝ |v|²). # Returning early also avoids the division by velocity_magnitude below, which @@ -143,6 +157,10 @@ def compute( _logger.debug("compute: zero apparent velocity, returning zero lift/drag force") return 0.0, 0.0, attack_angle_deg + lift_coefficient, drag_coefficient, area = self.interpolate( + attack_angle_deg, velocity_magnitude + ) + # Calculate the lift and drag forces lift_n = 0.5 * self.__fluid_density * lift_coefficient * area * velocity_magnitude**2 drag_n = 0.5 * self.__fluid_density * drag_coefficient * area * velocity_magnitude**2 @@ -157,37 +175,47 @@ def compute( ) return lift_n, drag_n, attack_angle_deg - def interpolate(self, attack_angle: float) -> Tuple[float, float, float]: - """Performs linear interpolation to estimate the lift coefficient, drag coefficient, and - area upon which the fluid acts, based on the provided angle of attack. + def interpolate(self, attack_angle: float, flow_speed: float) -> Tuple[float, float, float]: + """Estimate the lift coefficient, drag coefficient, and area acting on the medium. + + The coefficients are interpolated over the angle of attack and the Reynolds number + ``Re = flow_speed * chord / kinematic_viscosity``, so a foil's polar shifts with the flow + speed. Args: attack_angle (float): The angle of attack formed between the orientation angle of the medium and the direction of the apparent velocity, in degrees. + flow_speed (float): The magnitude of the apparent relative velocity between the + fluid and the medium, in meters per second, used to compute the Reynolds + number. Returns: Tuple[float, float, float]: A tuple of (lift_coefficient, drag_coefficient, area). Both coefficients are unitless; area is in meters squared. """ - # The foils are symmetric, so the lookup tables only store the positive AoA branch: + reynolds = flow_speed * self.__chord / self.__kinematic_viscosity + + # The foils are symmetric, so the lookup grids only store the positive AoA branch: # C_l is odd (C_l(-a) = -C_l(a)) and C_d is even (C_d(-a) = C_d(a)). We interpolate on - # |AoA| and re-apply the sign to lift. Beyond the table's max angle the foil is fully + # |AoA| and re-apply the sign to lift. Beyond the grid's max angle the foil is fully # stalled and outside the modeled regime, so we return zero (np.interp would otherwise # silently clamp to the endpoint, producing peak lift at all out-of-range angles → - # runaway thrust). Tables are assumed to start at 0°. + # runaway thrust). Grids are assumed to start at 0°. abs_attack_angle = abs(attack_angle) lift_sign = float(np.sign(attack_angle)) if abs_attack_angle > self.__lift_coefficients.max_angle: lift_coefficient = 0.0 else: - lift_coefficient = lift_sign * self.__lift_coefficients.interpolate(abs_attack_angle) + lift_coefficient = lift_sign * self.__lift_coefficients.interpolate( + abs_attack_angle, reynolds + ) if abs_attack_angle > self.__drag_coefficients.max_angle: drag_coefficient = 0.0 else: - drag_coefficient = self.__drag_coefficients.interpolate(abs_attack_angle) + drag_coefficient = self.__drag_coefficients.interpolate(abs_attack_angle, reynolds) # The foil area is modeled as constant (independent of angle of attack). area = self.__areas @@ -195,11 +223,11 @@ def interpolate(self, attack_angle: float) -> Tuple[float, float, float]: return lift_coefficient, drag_coefficient, area @property - def lift_coefficients(self) -> CoeffTable: + def lift_coefficients(self) -> CoeffGrid: return self.__lift_coefficients @property - def drag_coefficients(self) -> CoeffTable: + def drag_coefficients(self) -> CoeffGrid: return self.__drag_coefficients @property @@ -210,6 +238,14 @@ def areas(self) -> float: def fluid_density(self) -> float: return self.__fluid_density + @property + def chord(self) -> float: + return self.__chord + + @property + def kinematic_viscosity(self) -> float: + return self.__kinematic_viscosity + class HydroStaticsForceComputation: """Computes the hydrostatic restoring force. @@ -318,14 +354,16 @@ def __init__(self): BOAT_PROPERTIES.sail_drag_coeffs, BOAT_PROPERTIES.sail_areas, AIR_DENSITY, + WING_SAIL_CHORD, + AIR_KINEMATIC_VISCOSITY, ) - # TODO The trim tab needs its own coefficient tables and area in BOAT_PROPERTIES; - # the main wing's values are stand-ins. self.__tab = MediumForceComputation( BOAT_PROPERTIES.tab_lift_coeffs, BOAT_PROPERTIES.tab_drag_coeffs, BOAT_PROPERTIES.tab_areas, AIR_DENSITY, + TAB_CHORD, + AIR_KINEMATIC_VISCOSITY, ) self.__chord_m = WING_SAIL_CHORD @@ -388,7 +426,7 @@ def net_pitching_moment(self, alpha_rad: float, v_aw: float, delta_tab_rad: floa m_wing = n_perp * (x_cop - self.__mast_pivot_chord_m) alpha_tab_rad = alpha_rad - delta_tab_rad - lift_coefficient, _, tab_area = self.__tab.interpolate(math.degrees(alpha_tab_rad)) + lift_coefficient, _, tab_area = self.__tab.interpolate(math.degrees(alpha_tab_rad), v_aw) l_tab = 0.5 * self.__air_density * v_aw**2 * tab_area * lift_coefficient m_tab = l_tab * self.__boom_length_m return m_wing + m_tab @@ -503,6 +541,8 @@ def __init__(self): BOAT_PROPERTIES.rudder_drag_coeffs, BOAT_PROPERTIES.rudder_areas, WATER_DENSITY, + RUDDER_CHORD, + WATER_KINEMATIC_VISCOSITY, ) self.__keel = MediumForceComputation( @@ -510,6 +550,8 @@ def __init__(self): BOAT_PROPERTIES.keel_drag_coeffs, BOAT_PROPERTIES.keel_areas, WATER_DENSITY, + KEEL_CHORD, + WATER_KINEMATIC_VISCOSITY, ) self.__x_r, self.__z_r = RUDDER_CE_REL_TO_CG self.__x_k, self.__z_k = KEEL_CE_REL_TO_CG diff --git a/src/boat_simulator/scripts/build_airfoil_polars.py b/src/boat_simulator/scripts/build_airfoil_polars.py new file mode 100644 index 000000000..4c5ab451f --- /dev/null +++ b/src/boat_simulator/scripts/build_airfoil_polars.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python3 +"""Generate ``boat_simulator/common/airfoil_polars.py`` from airfoiltools.com polars. + +For each foil the simulator models, this script downloads the XFOIL polar CSVs for its NACA +section at a range of Reynolds numbers, resamples the pre-stall branch to whole-degree angles of +attack, and extends each polar out to 90 deg with the Viterna-Corrigan flat-plate post-stall +model (a sailboat's sail and control surfaces routinely operate far past stall -- e.g. a sail near +90 deg dead downwind -- so the raw XFOIL range of roughly +/-18 deg is not enough). The result is +written as an importable Python module of :class:`CoeffGrid` objects, keeping the large numeric +tables out of ``constants.py`` while staying diffable and free of package-data plumbing. + +This is a one-off developer tool: run it to (re)generate the committed data module. It is not +imported at runtime. + +Usage: + python3 src/boat_simulator/scripts/build_airfoil_polars.py + +Reynolds interpolation between the generated tables happens at runtime in :class:`CoeffGrid`. +""" + +from __future__ import annotations + +import math +import textwrap +import urllib.request +from dataclasses import dataclass +from typing import List, Tuple + +import numpy as np + +# Standard Reynolds numbers airfoiltools.com publishes XFOIL polars at. +REYNOLDS_NUMBERS: Tuple[int, ...] = (50_000, 100_000, 200_000, 500_000, 1_000_000) + +# airfoiltools polar-key template: xf-- (Ncrit 9, Mach 0). +_CSV_URL = "http://airfoiltools.com/polar/csv?polar=xf-{airfoil}-{reynolds}" + +# Whole-degree angle-of-attack samples appended past the XFOIL range via Viterna extrapolation. +_POST_STALL_ANGLES_DEG: Tuple[float, ...] = (20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 75.0, 90.0) + + +@dataclass(frozen=True) +class Foil: + """A modeled foil: which sim tables it feeds, its NACA section, and its aspect ratio. + + Attributes: + name: Prefix used for the generated module constants (e.g. ``"sail"`` -> ``SAIL_*``). + airfoil: airfoiltools airfoil id (e.g. ``"n0012-il"``). + aspect_ratio: Planform aspect ratio, used only to set the post-stall drag plateau + ``CD_max = 1.11 + 0.018 * AR`` (Viterna-Corrigan). Approximate -- refine with the real + geometry. + """ + + name: str + airfoil: str + aspect_ratio: float + + +# TODO: refine the aspect ratios from the real boat geometry; they only affect the post-stall drag +# plateau (CD_max). Keel and rudder share the NACA 0012 polars but keep separate grids because +# their chords -- and hence operating Reynolds numbers -- differ. +FOILS: Tuple[Foil, ...] = ( + Foil(name="sail", airfoil="naca0018-il", aspect_ratio=5.0), + Foil(name="tab", airfoil="naca0018-il", aspect_ratio=3.0), + Foil(name="rudder", airfoil="n0012-il", aspect_ratio=3.0), + Foil(name="keel", airfoil="n0012-il", aspect_ratio=3.0), +) + + +def fetch_polar(airfoil: str, reynolds: int) -> np.ndarray: + """Download one airfoiltools polar CSV and return its ``(N, 3)`` [alpha, Cl, Cd] rows.""" + url = _CSV_URL.format(airfoil=airfoil, reynolds=reynolds) + with urllib.request.urlopen(url, timeout=30) as response: # noqa: S310 (trusted host) + text = response.read().decode("utf-8") + + lines = text.splitlines() + header_idx = next(i for i, line in enumerate(lines) if line.startswith("Alpha,")) + rows: List[Tuple[float, float, float]] = [] + for line in lines[header_idx + 1 :]: + parts = line.split(",") + if len(parts) < 3 or not parts[0].strip(): + continue + alpha, cl, cd = float(parts[0]), float(parts[1]), float(parts[2]) + rows.append((alpha, cl, cd)) + return np.array(rows, dtype=float) + + +def _viterna_coefficients( + alpha_stall_rad: float, cl_stall: float, cd_stall: float, aspect_ratio: float +) -> Tuple[float, float, float, float]: + """Return the Viterna-Corrigan constants (A1, A2, B1, B2) matched at the stall point.""" + cd_max = 1.11 + 0.018 * aspect_ratio + b1 = cd_max + a1 = b1 / 2.0 + sin_s, cos_s = math.sin(alpha_stall_rad), math.cos(alpha_stall_rad) + b2 = (cd_stall - cd_max * sin_s**2) / cos_s + a2 = (cl_stall - cd_max * sin_s * cos_s) * sin_s / cos_s**2 + return a1, a2, b1, b2 + + +def build_tables(foil: Foil, polar: np.ndarray) -> Tuple[List[List[float]], List[List[float]]]: + """Resample the positive-AoA branch and extend to 90 deg; return (lift_rows, drag_rows).""" + positive = polar[polar[:, 0] >= 0.0] + positive = positive[np.argsort(positive[:, 0])] + alpha_fine, cl_fine, cd_fine = positive[:, 0], positive[:, 1], positive[:, 2] + + # Match Viterna at the last whole degree covered by XFOIL data, so the stored table is + # continuous across the XFOIL -> Viterna seam. + alpha_match_deg = float(math.floor(alpha_fine[-1])) + pre_stall_angles = np.arange(0.0, alpha_match_deg + 1.0, 1.0) + cl_pre = np.interp(pre_stall_angles, alpha_fine, cl_fine) + cd_pre = np.interp(pre_stall_angles, alpha_fine, cd_fine) + + a1, a2, b1, b2 = _viterna_coefficients( + math.radians(alpha_match_deg), + float(np.interp(alpha_match_deg, alpha_fine, cl_fine)), + float(np.interp(alpha_match_deg, alpha_fine, cd_fine)), + foil.aspect_ratio, + ) + + lift_rows = [[float(a), round(float(c), 4)] for a, c in zip(pre_stall_angles, cl_pre)] + drag_rows = [[float(a), round(float(c), 4)] for a, c in zip(pre_stall_angles, cd_pre)] + for angle_deg in _POST_STALL_ANGLES_DEG: + if angle_deg <= alpha_match_deg: + continue + rad = math.radians(angle_deg) + cl = a1 * math.sin(2 * rad) + a2 * math.cos(rad) ** 2 / math.sin(rad) + cd = b1 * math.sin(rad) ** 2 + b2 * math.cos(rad) + lift_rows.append([angle_deg, round(cl, 4)]) + drag_rows.append([angle_deg, round(cd, 4)]) + return lift_rows, drag_rows + + +def _format_grid(name: str, tables_rows: List[List[List[float]]]) -> str: + """Render a flake8-clean ``CoeffGrid`` literal for one foil coefficient (lift or drag).""" + table_literals = [] + for rows in tables_rows: + body = ",\n".join(f" [{a:.1f}, {v}]" for a, v in rows) + table_literals.append( + " CoeffTable(\n" + " np.array(\n" + " [\n" + + body + + ",\n" + " ],\n" + " dtype=np.float64,\n" + " )\n" + " )" + ) + joined = ",\n".join(table_literals) + return f"{name} = CoeffGrid(\n REYNOLDS_NUMBERS,\n (\n{joined},\n ),\n)" + + +def main() -> None: + blocks: List[str] = [] + for foil in FOILS: + lift_tables: List[List[List[float]]] = [] + drag_tables: List[List[List[float]]] = [] + for reynolds in REYNOLDS_NUMBERS: + print(f"Fetching {foil.airfoil} Re={reynolds} for {foil.name}...") + polar = fetch_polar(foil.airfoil, reynolds) + lift_rows, drag_rows = build_tables(foil, polar) + lift_tables.append(lift_rows) + drag_tables.append(drag_rows) + upper = foil.name.upper() + naca = "".join(ch for ch in foil.airfoil.split("-")[0] if ch.isdigit()) + blocks.append( + f"# {foil.name}: NACA {naca} " + f"(aspect ratio {foil.aspect_ratio}), Viterna-extended to 90 deg.\n" + + _format_grid(f"{upper}_LIFT_COEFFS", lift_tables) + + "\n\n" + + _format_grid(f"{upper}_DRAG_COEFFS", drag_tables) + ) + + header = textwrap.dedent( + '''\ + """Reynolds-dependent lift/drag coefficient grids for the boat simulator's foils. + + GENERATED FILE -- do not edit by hand. Regenerate with:: + + python3 src/boat_simulator/scripts/build_airfoil_polars.py + + Each grid holds one XFOIL polar per Reynolds number (from airfoiltools.com), resampled to + whole-degree angles of attack and extended past stall to 90 deg with the Viterna-Corrigan + model. Only the positive angle-of-attack branch is stored; the symmetric-foil sign + convention (Cl odd, Cd even) is applied by the force computation at lookup time. + """ + + import numpy as np + + from boat_simulator.common.types import CoeffGrid, CoeffTable + + # Ascending Reynolds-number breakpoints shared by every grid below. + REYNOLDS_NUMBERS = np.array( + ''' + ) + reynolds_literal = ( + " [" + ", ".join(f"{r}.0" for r in REYNOLDS_NUMBERS) + "],\n dtype=np.float64,\n)" + ) + + output = header + reynolds_literal + "\n\n\n" + "\n\n\n".join(blocks) + "\n" + + out_path = ( + "/workspaces/sailbot_workspace/src/boat_simulator/boat_simulator/common/airfoil_polars.py" + ) + with open(out_path, "w") as handle: + handle.write(output) + print(f"Wrote {out_path}") + + +if __name__ == "__main__": + main() diff --git a/src/boat_simulator/tests/unit/common/test_types.py b/src/boat_simulator/tests/unit/common/test_types.py new file mode 100644 index 000000000..b1d804fd1 --- /dev/null +++ b/src/boat_simulator/tests/unit/common/test_types.py @@ -0,0 +1,80 @@ +"""Tests for the CoeffGrid Reynolds-dependent lookup in boat_simulator/common/types.py.""" + +import math + +import numpy as np +import pytest + +from boat_simulator.common.types import CoeffGrid, CoeffTable + + +def _linear_table(slope: float) -> CoeffTable: + """A table whose value equals ``slope * angle`` at angles 0 and 10 degrees.""" + return CoeffTable(np.array([[0.0, 0.0], [10.0, 10.0 * slope]], dtype=np.float64)) + + +class TestCoeffGrid: + def test_from_single_is_reynolds_independent(self): + grid = CoeffGrid.from_single(_linear_table(1.0)) + # Same value regardless of Reynolds number. + for reynolds in (1.0, 1e3, 1e6, 1e9): + assert math.isclose(grid.interpolate(5.0, reynolds), 5.0) + + def test_exact_reynolds_hit_returns_that_polar(self): + grid = CoeffGrid( + np.array([100.0, 1000.0], dtype=np.float64), + (_linear_table(1.0), _linear_table(2.0)), + ) + assert math.isclose(grid.interpolate(5.0, 100.0), 5.0) # slope 1 polar + assert math.isclose(grid.interpolate(5.0, 1000.0), 10.0) # slope 2 polar + + def test_interpolates_linearly_in_log_reynolds(self): + # Values 0.0 and 2.0 at Re 100 and 10000; the geometric-mean Re (1000) is the log-midpoint. + grid = CoeffGrid( + np.array([100.0, 10000.0], dtype=np.float64), + (_linear_table(0.0), _linear_table(2.0)), + ) + midpoint = grid.interpolate(10.0, 1000.0) + assert math.isclose(midpoint, 10.0) # halfway between 0 and 20 + + def test_clamps_outside_reynolds_range(self): + grid = CoeffGrid( + np.array([100.0, 1000.0], dtype=np.float64), + (_linear_table(1.0), _linear_table(2.0)), + ) + # Below the lowest / above the highest Reynolds clamp to the end polars. + assert math.isclose(grid.interpolate(5.0, 1.0), 5.0) + assert math.isclose(grid.interpolate(5.0, 1e9), 10.0) + + def test_non_positive_reynolds_uses_lowest(self): + grid = CoeffGrid( + np.array([100.0, 1000.0], dtype=np.float64), + (_linear_table(1.0), _linear_table(2.0)), + ) + assert math.isclose(grid.interpolate(5.0, 0.0), 5.0) + assert math.isclose(grid.interpolate(5.0, -10.0), 5.0) + + def test_max_angle_is_minimum_over_tables(self): + short = CoeffTable(np.array([[0.0, 0.0], [45.0, 1.0]], dtype=np.float64)) + tall = CoeffTable(np.array([[0.0, 0.0], [90.0, 1.0]], dtype=np.float64)) + grid = CoeffGrid(np.array([100.0, 1000.0], dtype=np.float64), (short, tall)) + assert grid.max_angle == 45.0 + + @pytest.mark.parametrize( + "reynolds, tables", + [ + (np.array([100.0]), (_linear_table(1.0), _linear_table(2.0))), # count mismatch + (np.array([1000.0, 100.0]), (_linear_table(1.0), _linear_table(2.0))), # not ascending + (np.array([0.0, 100.0]), (_linear_table(1.0), _linear_table(2.0))), # non-positive Re + (np.array([[100.0]]), (_linear_table(1.0),)), # not 1-D + ], + ) + def test_invalid_construction_raises(self, reynolds, tables): + with pytest.raises(ValueError): + CoeffGrid(reynolds, tables) + + def test_single_table_grid_matches_underlying_table(self): + table = _linear_table(1.5) + grid = CoeffGrid.from_single(table) + for angle in (0.0, 3.0, 7.5, 10.0): + assert math.isclose(grid.interpolate(angle, 500.0), table.interpolate(angle)) diff --git a/src/boat_simulator/tests/unit/nodes/physics_engine/test_medium_force_reynolds.py b/src/boat_simulator/tests/unit/nodes/physics_engine/test_medium_force_reynolds.py new file mode 100644 index 000000000..094eefa12 --- /dev/null +++ b/src/boat_simulator/tests/unit/nodes/physics_engine/test_medium_force_reynolds.py @@ -0,0 +1,82 @@ +"""Tests for Reynolds-number-dependent coefficient lookup in MediumForceComputation.""" + +import math + +import numpy as np +import pytest + +from boat_simulator.common.types import CoeffGrid, CoeffTable +from boat_simulator.nodes.physics_engine.fluid_forces import MediumForceComputation + + +def _grid(low_value: float, high_value: float) -> CoeffGrid: + """Two-Re grid: ``low_value`` at Re=100, ``high_value`` at Re=1000, sampled at 5 deg.""" + low = CoeffTable(np.array([[0.0, 0.0], [5.0, low_value], [90.0, 0.0]], dtype=np.float64)) + high = CoeffTable(np.array([[0.0, 0.0], [5.0, high_value], [90.0, 0.0]], dtype=np.float64)) + return CoeffGrid(np.array([100.0, 1000.0], dtype=np.float64), (low, high)) + + +def _medium(lift: CoeffGrid, drag: CoeffGrid) -> MediumForceComputation: + # chord = 1 and kinematic_viscosity = 1 make Reynolds number equal to the flow speed. + return MediumForceComputation( + lift_coefficients=lift, + drag_coefficients=drag, + areas=1.0, + fluid_density=1.0, + chord=1.0, + kinematic_viscosity=1.0, + ) + + +class TestMediumForceReynolds: + def test_lift_tracks_reynolds_via_flow_speed(self): + # Lift rises with Reynolds; drag falls with Reynolds (typical foil behaviour). + medium = _medium(_grid(1.0, 2.0), _grid(0.5, 0.1)) + cl_low, cd_low, _ = medium.interpolate(5.0, flow_speed=100.0) + cl_high, cd_high, _ = medium.interpolate(5.0, flow_speed=1000.0) + assert math.isclose(cl_low, 1.0) and math.isclose(cl_high, 2.0) + assert cl_high > cl_low + assert cd_high < cd_low + + def test_force_uses_reynolds_dependent_coefficient(self): + medium = _medium(_grid(1.0, 2.0), _grid(0.1, 0.1)) + from boat_simulator.common.types import Vec2 + + # Flow along +x with the medium oriented so the angle of attack is 5 deg. + lift_low, _, _ = medium.compute(Vec2.from_xy(100.0 * math.cos(math.radians(5.0)), + 100.0 * math.sin(math.radians(5.0))), 0.0) + lift_high, _, _ = medium.compute(Vec2.from_xy(1000.0 * math.cos(math.radians(5.0)), + 1000.0 * math.sin(math.radians(5.0))), 0.0) + # Force ~ 0.5*rho*Cl*A*v^2. Speed x10 -> v^2 x100, and Cl doubles (Re x10), so lift x200. + assert math.isclose(lift_high / lift_low, 200.0, rel_tol=1e-6) + + def test_beyond_max_angle_returns_zero(self): + medium = _medium(_grid(1.0, 2.0), _grid(0.5, 0.5)) + cl, cd, _ = medium.interpolate(95.0, flow_speed=500.0) # > 90 deg tabulated max + assert cl == 0.0 and cd == 0.0 + + def test_symmetric_sign_convention(self): + medium = _medium(_grid(1.0, 2.0), _grid(0.5, 0.5)) + cl_pos, cd_pos, _ = medium.interpolate(5.0, flow_speed=1000.0) + cl_neg, cd_neg, _ = medium.interpolate(-5.0, flow_speed=1000.0) + assert math.isclose(cl_neg, -cl_pos) # lift is odd in angle of attack + assert math.isclose(cd_neg, cd_pos) # drag is even in angle of attack + + def test_zero_velocity_returns_zero_force(self): + medium = _medium(_grid(1.0, 2.0), _grid(0.5, 0.5)) + from boat_simulator.common.types import Vec2 + + lift, drag, _ = medium.compute(Vec2.from_xy(0.0, 0.0), 30.0) + assert lift == 0.0 and drag == 0.0 + + +@pytest.mark.parametrize("reynolds", [50_000, 200_000, 1_000_000]) +def test_real_keel_grid_lift_increases_with_reynolds_near_stall(reynolds): + """Sanity check on the generated NACA 0012 keel grid: max lift grows with Reynolds.""" + from boat_simulator.common.constants import BOAT_PROPERTIES + + grid = BOAT_PROPERTIES.keel_lift_coeffs + # Near 15 deg the low-Re polar is already stalled while the high-Re polar still climbs. + low = grid.interpolate(15.0, 50_000) + high = grid.interpolate(15.0, 1_000_000) + assert high > low From 38358255d353b5a02825a1c6fd9abb842a7486dd Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 13:28:18 -0700 Subject: [PATCH 15/21] Update the chord length --- .../boat_simulator/common/constants.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/constants.py b/src/boat_simulator/boat_simulator/common/constants.py index ef4827cd4..a6d9dfa5d 100644 --- a/src/boat_simulator/boat_simulator/common/constants.py +++ b/src/boat_simulator/boat_simulator/common/constants.py @@ -141,13 +141,16 @@ class BoatProperties: # Metacentric height used in testing METACENTRIC_HEIGHT = 0.5 # Units: meters -# Derive the mean chord from the real wingsail geometry. -WING_SAIL_CHORD = 1.5 # Units: meters +# The mean chord from the real wingsail geometry. +WING_SAIL_CHORD = 1.0 # Units: meters -# Mean chord lengths of the remaining foils, used only to compute each foil's Reynolds number. -# TODO: Replace these placeholders with the real measured mean chords. -TAB_CHORD = 0.3 # Units: meters -RUDDER_CHORD = 0.2 # Units: meters +# Mean chord lengths for the trim tab +TAB_CHORD = 0.5 # Units: meters + +# Mean chord lengths for the rudder +RUDDER_CHORD = 0.25 # Units: meters + +# TODO: Update the keel's chord length KEEL_CHORD = 0.4 # Units: meters # Measure the distance from the mast axis to the tab's aero center. From 84025614442b0e2db419ff3fa8030cba527eaf7a Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Thu, 23 Jul 2026 16:24:09 -0700 Subject: [PATCH 16/21] Updating comments --- .../boat_simulator/common/airfoil_polars.py | 3 ++- .../boat_simulator/common/types.py | 19 ++++++++----------- .../tests/unit/common/test_types.py | 11 ----------- 3 files changed, 10 insertions(+), 23 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/airfoil_polars.py b/src/boat_simulator/boat_simulator/common/airfoil_polars.py index c2d2b6c28..4e5dedbaf 100644 --- a/src/boat_simulator/boat_simulator/common/airfoil_polars.py +++ b/src/boat_simulator/boat_simulator/common/airfoil_polars.py @@ -1,6 +1,7 @@ """Reynolds-dependent lift/drag coefficient grids for the boat simulator's foils. -GENERATED FILE -- do not edit by hand. Regenerate with:: + +Do not edit by hand. Regenerate with: python3 src/boat_simulator/scripts/build_airfoil_polars.py diff --git a/src/boat_simulator/boat_simulator/common/types.py b/src/boat_simulator/boat_simulator/common/types.py index 24d7973ea..68750c3bf 100644 --- a/src/boat_simulator/boat_simulator/common/types.py +++ b/src/boat_simulator/boat_simulator/common/types.py @@ -495,22 +495,28 @@ class CoeffGrid: def __init__(self, reynolds: ArrayLike, tables: Sequence[CoeffTable]) -> None: reynolds_arr = np.asarray(reynolds, dtype=float) + + # Safety checks if reynolds_arr.ndim != 1 or reynolds_arr.shape[0] < 1: raise ValueError( f"expected a non-empty 1-D Reynolds array, got shape {reynolds_arr.shape}" ) + if np.any(reynolds_arr <= 0): raise ValueError( "Reynolds numbers must be positive (interpolation is done in log space)" ) + if np.any(np.diff(reynolds_arr) <= 0): raise ValueError("Reynolds-number column must be sorted in strictly ascending order") + tables_tuple = tuple(tables) if len(tables_tuple) != reynolds_arr.shape[0]: raise ValueError( "expected one CoeffTable per Reynolds number, got " f"{len(tables_tuple)} tables for {reynolds_arr.shape[0]} Reynolds values" ) + if not all(isinstance(table, CoeffTable) for table in tables_tuple): raise ValueError("tables must all be CoeffTable instances") @@ -519,15 +525,6 @@ def __init__(self, reynolds: ArrayLike, tables: Sequence[CoeffTable]) -> None: object.__setattr__(self, "reynolds", reynolds_arr) object.__setattr__(self, "tables", tables_tuple) - @classmethod - def from_single(cls, table: CoeffTable, reynolds: float = 1.0) -> CoeffGrid: - """Build a Reynolds-independent grid from a single table. - - The resulting grid returns ``table``'s value at any Reynolds number, letting a - Reynolds-independent surface share the same lookup interface as a multi-Re one. - """ - return cls(np.array([reynolds], dtype=float), (table,)) - @property def max_angle(self) -> float: """Smallest per-table maximum angle (degrees). @@ -543,9 +540,9 @@ def min_angle(self) -> float: return max(table.min_angle for table in self.tables) def interpolate(self, angle: float, reynolds: float) -> float: - """Interpolate the coefficient at ``angle`` (degrees) and ``reynolds``. + """Interpolate the coefficient at angle (degrees) and reynolds. - Linear in angle of attack within each per-Re table, then linear in ``log(Re)`` between + Linear in angle of attack within each per-Re table, then linear in log(Re) between tables. A single-Re grid is Reynolds-independent (the lone table's value is returned). Non-positive Reynolds numbers are treated as the lowest tabulated value. """ diff --git a/src/boat_simulator/tests/unit/common/test_types.py b/src/boat_simulator/tests/unit/common/test_types.py index b1d804fd1..24f4c4747 100644 --- a/src/boat_simulator/tests/unit/common/test_types.py +++ b/src/boat_simulator/tests/unit/common/test_types.py @@ -14,11 +14,6 @@ def _linear_table(slope: float) -> CoeffTable: class TestCoeffGrid: - def test_from_single_is_reynolds_independent(self): - grid = CoeffGrid.from_single(_linear_table(1.0)) - # Same value regardless of Reynolds number. - for reynolds in (1.0, 1e3, 1e6, 1e9): - assert math.isclose(grid.interpolate(5.0, reynolds), 5.0) def test_exact_reynolds_hit_returns_that_polar(self): grid = CoeffGrid( @@ -72,9 +67,3 @@ def test_max_angle_is_minimum_over_tables(self): def test_invalid_construction_raises(self, reynolds, tables): with pytest.raises(ValueError): CoeffGrid(reynolds, tables) - - def test_single_table_grid_matches_underlying_table(self): - table = _linear_table(1.5) - grid = CoeffGrid.from_single(table) - for angle in (0.0, 3.0, 7.5, 10.0): - assert math.isclose(grid.interpolate(angle, 500.0), table.interpolate(angle)) From 6d309b55402710629d38ba9e9e8595c2b91326df Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Fri, 24 Jul 2026 00:34:17 -0700 Subject: [PATCH 17/21] Switch reynolds number selection to closest one rather than log reynold --- .../boat_simulator/common/types.py | 21 +++++++++---------- .../nodes/physics_engine/fluid_forces.py | 7 +------ 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/src/boat_simulator/boat_simulator/common/types.py b/src/boat_simulator/boat_simulator/common/types.py index 68750c3bf..9ed737088 100644 --- a/src/boat_simulator/boat_simulator/common/types.py +++ b/src/boat_simulator/boat_simulator/common/types.py @@ -539,19 +539,18 @@ def min_angle(self) -> float: """Largest per-table minimum angle (degrees).""" return max(table.min_angle for table in self.tables) - def interpolate(self, angle: float, reynolds: float) -> float: - """Interpolate the coefficient at angle (degrees) and reynolds. + def interpolate(self, angle: float, reynolds_input: float) -> float: + """Interpolate the coefficient at ``angle`` (degrees) for the given Reynolds number. - Linear in angle of attack within each per-Re table, then linear in log(Re) between - tables. A single-Re grid is Reynolds-independent (the lone table's value is returned). - Non-positive Reynolds numbers are treated as the lowest tabulated value. + The angle of attack is interpolated within the polar whose tabulated Reynolds number is + closest to the closest reynolds array index. """ - values = np.array([table.interpolate(angle) for table in self.tables], dtype=float) - if self.reynolds.shape[0] == 1: - return float(values[0]) - if reynolds <= 0: - reynolds = float(self.reynolds[0]) - return float(np.interp(np.log(reynolds), np.log(self.reynolds), values)) + if reynolds_input <= 0: + index = 0 + else: + index = int(np.argmin(np.abs(self.reynolds - reynolds_input))) + + return float(self.tables[index].interpolate(angle)) def __len__(self) -> int: return int(self.reynolds.shape[0]) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index 5709fec88..764403ce5 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -196,12 +196,7 @@ def interpolate(self, attack_angle: float, flow_speed: float) -> Tuple[float, fl reynolds = flow_speed * self.__chord / self.__kinematic_viscosity - # The foils are symmetric, so the lookup grids only store the positive AoA branch: - # C_l is odd (C_l(-a) = -C_l(a)) and C_d is even (C_d(-a) = C_d(a)). We interpolate on - # |AoA| and re-apply the sign to lift. Beyond the grid's max angle the foil is fully - # stalled and outside the modeled regime, so we return zero (np.interp would otherwise - # silently clamp to the endpoint, producing peak lift at all out-of-range angles → - # runaway thrust). Grids are assumed to start at 0°. + # The foils are symmetric, so the lookup grids only store the positive AoA branch abs_attack_angle = abs(attack_angle) lift_sign = float(np.sign(attack_angle)) From c02c76581689d97be7d1e34cdc971ed78cfe12ad Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Fri, 24 Jul 2026 00:34:35 -0700 Subject: [PATCH 18/21] Add test case to find the nearest reynolds number --- src/boat_simulator/tests/unit/common/test_types.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/boat_simulator/tests/unit/common/test_types.py b/src/boat_simulator/tests/unit/common/test_types.py index 24f4c4747..e08e81ed0 100644 --- a/src/boat_simulator/tests/unit/common/test_types.py +++ b/src/boat_simulator/tests/unit/common/test_types.py @@ -23,14 +23,16 @@ def test_exact_reynolds_hit_returns_that_polar(self): assert math.isclose(grid.interpolate(5.0, 100.0), 5.0) # slope 1 polar assert math.isclose(grid.interpolate(5.0, 1000.0), 10.0) # slope 2 polar - def test_interpolates_linearly_in_log_reynolds(self): - # Values 0.0 and 2.0 at Re 100 and 10000; the geometric-mean Re (1000) is the log-midpoint. + def test_selects_nearest_reynolds_polar(self): + # Slope-1 polar at Re=100, slope-2 polar at Re=1000; value at 5 deg is 5.0 vs 10.0. grid = CoeffGrid( - np.array([100.0, 10000.0], dtype=np.float64), - (_linear_table(0.0), _linear_table(2.0)), + np.array([100.0, 1000.0], dtype=np.float64), + (_linear_table(1.0), _linear_table(2.0)), ) - midpoint = grid.interpolate(10.0, 1000.0) - assert math.isclose(midpoint, 10.0) # halfway between 0 and 20 + assert math.isclose(grid.interpolate(5.0, 400.0), 5.0) # closer to Re=100 + assert math.isclose(grid.interpolate(5.0, 700.0), 10.0) # closer to Re=1000 + # Equidistant (Re=550): np.argmin picks the first (lower-Re) polar. + assert math.isclose(grid.interpolate(5.0, 550.0), 5.0) def test_clamps_outside_reynolds_range(self): grid = CoeffGrid( From 3dcb7f7bda4f8d8ee65f5ab1919f09b64e3ea9c3 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Fri, 24 Jul 2026 00:34:49 -0700 Subject: [PATCH 19/21] Add trim tab and rudder angle to the logs --- .../boat_simulator/nodes/physics_engine/model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py index 527d8d57b..0a247fb30 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/model.py @@ -113,7 +113,9 @@ def step( f"yaw_rate={math.degrees(self.nu.r):7.2f} deg/s\n" f" accel u_dot={self.nu_dot.x:8.3f} m/s^2 v_dot={self.nu_dot.y:8.3f} m/s^2 " f"roll_acc={math.degrees(self.nu_dot.p):7.2f} deg/s^2 " - f"yaw_acc={math.degrees(self.nu_dot.r):7.2f} deg/s^2" + f"yaw_acc={math.degrees(self.nu_dot.r):7.2f} deg/s^2\n" + f" trim_tab_angle={trim_tab_angle.degrees:7.2f} deg " + f"rudder_angle={rudder_angle.degrees:7.2f} deg" ) @property From 5e83cd9ce16a79769f925218e7e396017c831b11 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Fri, 24 Jul 2026 01:01:04 -0700 Subject: [PATCH 20/21] How to run the building of airfoils_polars.py file --- src/boat_simulator/README.md | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/boat_simulator/README.md b/src/boat_simulator/README.md index 7e4c88b65..87e9208c4 100644 --- a/src/boat_simulator/README.md +++ b/src/boat_simulator/README.md @@ -31,6 +31,33 @@ ros2 launch boat_simulator main_launch.py [OPTIONS]... To see a list of options for simulator configuration, add the `-s` flag at the end of the above command. +## Regenerating airfoil coefficient data + +The simulator's lift and drag forces depend on each foil's coefficients, which +vary with both the angle of attack, the Reynolds number (`Re = |v| * chord / +kinematic_viscosity`) and the mach number. For our calculations, we will consider the Reynolds number and angle of +attack Rather than a single hand-tuned curve per foil, the coefficients live in +[`common/airfoil_polars.py`](./boat_simulator/common/airfoil_polars.py) as a set of real airfoil polars +— one per Reynolds number — that the force computation selects from at runtime by picking the closest +tabulated Reynolds number. + +That module is **generated**, not edited by hand. The +[`scripts/build_airfoil_polars.py`](./scripts/build_airfoil_polars.py) script downloads XFOIL polars from +[airfoiltools.com](http://airfoiltools.com) for each foil's NACA section +(keel and rudder use NACA 0012; wingsail and trim tab use NACA 0018) at Reynolds numbers 50k–1M, +resamples the pre-stall branch to whole-degree angles of attack, and extends each polar out to 90° with the +Viterna–Corrigan post-stall model +([More details on Viterna extrapolation](https://www.simis.io/docs/aerodynamic-loads-viterna-extrapolation)) +(a sail sits near 90° when running dead downwind, well beyond XFOIL's ~±18° range). + +Regenerate the data module after changing which foils, sections, Reynolds +numbers, or aspect ratios are modeled (all configured near the top of the +script). It requires network access and is never imported at runtime: + +``` shell +python3 src/boat_simulator/scripts/build_airfoil_polars.py +``` + ## Test Run the `test` task in the Sailbot Workspace. See From 94c4b20b79beac619f2c49bbf2f43d14ff708bc8 Mon Sep 17 00:00:00 2001 From: Akshanjay Kompelli Date: Mon, 27 Jul 2026 21:47:36 -0700 Subject: [PATCH 21/21] Add bound_to_180 import statement --- .../boat_simulator/nodes/physics_engine/fluid_forces.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py index 73f5720a8..ab3331367 100644 --- a/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py +++ b/src/boat_simulator/boat_simulator/nodes/physics_engine/fluid_forces.py @@ -34,6 +34,7 @@ Velocity, ) from boat_simulator.common.types import CoeffGrid, Vec2, Vec4 +from boat_simulator.common.utils import bound_to_180 _logger = get_logger(__name__)