PythonModels/components/amesim_mechanical.py

@@ -125,3 +125,44 @@ class AmesimMassFrictionEndstops:
     def viscous_friction_force(self, velocity_m_s: float) -> float:
         return -self.viscous_friction_n_per_m_per_s * velocity_m_s
 
+    def windage_force(self, velocity_m_s: float) -> float:
+        return -self.windage_n_per_m2_per_s2 * velocity_m_s * abs(velocity_m_s)
+
+    def dry_friction_force(self, velocity_m_s: float) -> float:
+        if velocity_m_s > 0.0:
+            return -self.coulomb_friction_n
+        if velocity_m_s < 0.0:
+            return self.coulomb_friction_n
+        return 0.0
+
+    def limit_contact_force(self, displacement_m: float, velocity_m_s: float) -> float:
+        lower_force = self.lower_static_force_magnitude(displacement_m)
+        if lower_force > 0.0:
+            lower_force += max(-self.lower_damping_n_per_m_per_s * velocity_m_s, 0.0)
+
+        upper_force = self.upper_static_force_magnitude(displacement_m)
+        if upper_force > 0.0:
+            upper_force += max(self.upper_damping_n_per_m_per_s * velocity_m_s, 0.0)
+
+        return lower_force - upper_force
+
+    def derivatives(
+        self,
+        *,
+        velocity_m_s: float,
+        displacement_m: float,
+        port_1_force_n: float = 0.0,
+        port_2_force_n: float = 0.0,
+        external_force_n: float = 0.0,
+    ) -> tuple[float, float]:
+        total_force = (
+            port_1_force_n
+            + port_2_force_n
+            + external_force_n
+            + self.viscous_friction_force(velocity_m_s)
+            + self.windage_force(velocity_m_s)
+            + self.dry_friction_force(velocity_m_s)
+            + self.limit_contact_force(displacement_m, velocity_m_s)
+        )
+        return total_force / self.mass_kg, velocity_m_s
+

PythonModels/systems/test_mql.py

@@ -3865,6 +3865,7 @@ class TestMqlSystem:
         self.pnl00r_assembly = self._build_pnl00r_assembly()
         self.node3_assembly = self._build_node3_assembly()
         self.node4_assembly = self._build_node4_assembly()
+        self.mechanical_assembly = self._build_mechanical_assembly()
         self.pneumatic_assembly = self._build_pneumatic_assembly()
         pneumatic_components = self._pneumatic_components_by_alias()
         for spec in COMPONENT_SPECS:
@@ -3887,6 +3888,13 @@ class TestMqlSystem:
 
         return build_test_mql_pneumatic_assembly()
 
+    @staticmethod
+    def _build_mechanical_assembly():
+        # Local import avoids a module cycle: test_mql_config reads constants from this module.
+        from PythonModels.systems.test_mql_mechanical import build_test_mql_mechanical_assembly
+
+        return build_test_mql_mechanical_assembly()
+
     def _build_pnl0001_assembly(self):
         from PythonModels.systems.test_mql_pneumatic_lines import (
             build_test_mql_pnl0001_assembly,
@@ -4112,6 +4120,27 @@ class TestMqlSystem:
     def pneumatic_state_vector(self) -> list[float]:
         return self.network.initial_state_vector()
 
+    def mechanical_mass_closure(self):
+        from PythonModels.systems.test_mql_mechanical import TestMqlMechanicalMassClosure
+
+        return TestMqlMechanicalMassClosure(self.mechanical_assembly)
+
+    def mechanical_state_vector(self) -> list[float]:
+        return self.mechanical_mass_closure().initial_state_vector()
+
+    def simulate_mechanical_masses(
+        self,
+        config: SolveIVPConfig | None = None,
+        t_eval: list[float] | None = None,
+    ):
+        closure = self.mechanical_mass_closure()
+        return integrate_ode(
+            rhs=lambda t, state: closure.rhs(state),
+            initial_state=closure.initial_state_vector(),
+            config=config or SolveIVPConfig(),
+            t_eval=t_eval,
+        )
+
     @staticmethod
     def pneumatic_branch_spec(
         *,

PythonModels/systems/test_mql_mechanical.py

@@ -63,6 +63,8 @@ class TestMqlMassEndstopSpec:
     stribeck_constant_m_s: float
     use_friction: bool
     stop_type: int
+    initial_velocity_m_s: float
+    initial_displacement_m: float
     data_paths: tuple[str, ...]
 
     def endstop(self) -> AmesimMassFrictionEndstops:
@@ -141,6 +143,64 @@ class TestMqlMechanicalAssembly:
         )
 
 
+@dataclass(frozen=True)
+class TestMqlMechanicalMassState:
+    alias: str
+    velocity_m_s: float
+    displacement_m: float
+
+    def as_vector(self) -> list[float]:
+        return [self.velocity_m_s, self.displacement_m]
+
+
+@dataclass(frozen=True)
+class TestMqlMechanicalMassSnapshot:
+    states: tuple[TestMqlMechanicalMassState, ...]
+
+    @property
+    def state_count(self) -> int:
+        return 2 * len(self.states)
+
+
+class TestMqlMechanicalMassClosure:
+    def __init__(self, assembly: TestMqlMechanicalAssembly) -> None:
+        self.assembly = assembly
+        self.mass_aliases = tuple(assembly.masses)
+
+    def initial_state_vector(self) -> list[float]:
+        state: list[float] = []
+        for alias in self.mass_aliases:
+            spec = self.assembly.masses[alias]
+            state.extend([spec.initial_velocity_m_s, spec.initial_displacement_m])
+        return state
+
+    def snapshot(self, state_vector: list[float] | None = None) -> TestMqlMechanicalMassSnapshot:
+        values = self.initial_state_vector() if state_vector is None else list(state_vector)
+        if len(values) != 2 * len(self.mass_aliases):
+            raise ValueError("mechanical mass state vector requires two values per mass")
+        states = tuple(
+            TestMqlMechanicalMassState(
+                alias=alias,
+                velocity_m_s=values[2 * index],
+                displacement_m=values[2 * index + 1],
+            )
+            for index, alias in enumerate(self.mass_aliases)
+        )
+        return TestMqlMechanicalMassSnapshot(states=states)
+
+    def rhs(self, state_vector: list[float]) -> list[float]:
+        snapshot = self.snapshot(state_vector)
+        derivatives: list[float] = []
+        for state in snapshot.states:
+            mass = self.assembly.masses[state.alias].endstop()
+            acceleration, velocity = mass.derivatives(
+                velocity_m_s=state.velocity_m_s,
+                displacement_m=state.displacement_m,
+            )
+            derivatives.extend([acceleration, velocity])
+        return derivatives
+
+
 def build_test_mql_mechanical_assembly(
     config: TestMqlConfig | None = None,
     amesim_results: AmesimResults | None = None,
@@ -155,7 +215,7 @@ def build_test_mql_mechanical_assembly(
         for component in config.components_by_submodel("PNRP17")
     }
     masses = {
-        component.alias: _build_mass(component, variable_catalog)
+        component.alias: _build_mass(component, variable_catalog, amesim_results)
         for component in config.components_by_submodel("MECMAS21")
     }
     elastic_endstops = {
@@ -202,6 +262,7 @@ def _build_piston(
 def _build_mass(
     component: TestMqlResolvedComponent,
     variable_catalog: TestMqlVariableCatalog | None,
+    amesim_results: AmesimResults | None,
 ) -> TestMqlMassEndstopSpec:
     return TestMqlMassEndstopSpec(
         alias=component.alias,
@@ -224,6 +285,8 @@ def _build_mass(
         stribeck_constant_m_s=component.parameter_value("astrib"),
         use_friction=bool(int(component.parameter_value("useFriction"))),
         stop_type=int(component.parameter_value("stoptype")),
+        initial_velocity_m_s=_initial_value(amesim_results, f"v1@{component.alias}"),
+        initial_displacement_m=_initial_value(amesim_results, f"x1@{component.alias}"),
         data_paths=_data_paths(variable_catalog, component.alias),
     )
 
@@ -263,6 +326,12 @@ def n_per_mm_per_s_to_n_per_m_per_s(value: float) -> float:
     return value * N_PER_MM_PER_S_TO_N_PER_M_PER_S
 
 
+def _initial_value(amesim_results: AmesimResults | None, data_path: str) -> float:
+    if amesim_results is None:
+        return 0.0
+    return float(amesim_results.series(data_path)[0])
+
+
 def _data_paths(
     variable_catalog: TestMqlVariableCatalog | None,
     alias: str,

tests/test_test_mql_mechanical.py

@@ -4,8 +4,11 @@ import unittest
 from pathlib import Path
 
 from PythonModels.reporting.amesim_results import load_test_mql_amesim_results
+from PythonModels.core.solver import SolveIVPConfig
 from PythonModels.reporting.test_mql_variables import build_test_mql_variable_catalog
+from PythonModels.systems.test_mql import TestMqlSystem
 from PythonModels.systems.test_mql_mechanical import (
+    TestMqlMechanicalMassClosure,
     build_test_mql_mechanical_assembly,
     circular_area,
     mm_to_m,
@@ -37,6 +40,43 @@ class TestMqlMechanicalAssemblyTests(unittest.TestCase):
         self.assertEqual(len(self.assembly.force_connectors), 2)
         self.assertEqual(self.assembly.component_count, 46)
 
+    def test_mechanical_mass_closure_exposes_all_mecmas21_states(self) -> None:
+        closure = TestMqlMechanicalMassClosure(self.assembly)
+        state = closure.initial_state_vector()
+        snapshot = closure.snapshot(state)
+        rhs = closure.rhs(state)
+
+        self.assertEqual(closure.mass_aliases, tuple(self.assembly.masses))
+        self.assertEqual(len(closure.mass_aliases), 10)
+        self.assertEqual(len(state), 20)
+        self.assertEqual(snapshot.state_count, 20)
+        self.assertEqual(len(rhs), 20)
+        self.assertEqual(snapshot.states[0].alias, "mass_friction_endstops_10")
+        self.assertAlmostEqual(
+            snapshot.states[0].velocity_m_s,
+            self.amesim_results.series("v1@mass_friction_endstops_10")[0],
+        )
+        self.assertAlmostEqual(
+            snapshot.states[0].displacement_m,
+            self.amesim_results.series("x1@mass_friction_endstops_10")[0],
+        )
+        for index, mass_state in enumerate(snapshot.states):
+            self.assertAlmostEqual(rhs[2 * index + 1], mass_state.velocity_m_s)
+
+    def test_system_simulates_mechanical_mass_state_closure(self) -> None:
+        system = TestMqlSystem()
+
+        solution = system.simulate_mechanical_masses(
+            config=SolveIVPConfig(t_start=0.0, t_stop=1.0e-5, max_step=1.0e-6),
+            t_eval=[0.0, 1.0e-5],
+        )
+
+        self.assertTrue(solution.success)
+        self.assertEqual(len(system.mechanical_state_vector()), 20)
+        self.assertEqual(len(solution.t), 2)
+        self.assertEqual(len(solution.y), 20)
+        self.assertEqual([len(row) for row in solution.y], [2] * 20)
+
     def test_piston_geometry_is_converted_to_si_units(self) -> None:
         piston = self.assembly.pistons["pn_brp2_8"]