PythonModels/systems/test_mql.py

@@ -3921,6 +3921,9 @@ class TestMqlFullStateClosure:
         )
 
     def rhs(self, state_vector: list[float]) -> list[float]:
+        return self.rhs_at(0.0, state_vector)
+
+    def rhs_at(self, time_s: float, state_vector: list[float]) -> list[float]:
         pneumatic_state, mechanical_state = self._split_state(list(state_vector))
         pneumatic_snapshot = self.pneumatic_closure.snapshot(pneumatic_state)
         mechanical_snapshot = self.mechanical_closure.snapshot(mechanical_state)
@@ -3929,14 +3932,20 @@ class TestMqlFullStateClosure:
             *self.mechanical_closure.rhs(
                 mechanical_state,
                 force_by_mass_alias=self._force_by_mass_alias(
+                    time_s,
                     pneumatic_snapshot,
                     mechanical_snapshot,
                 ),
+                constrained_mass_aliases={
+                    "mass_friction_endstops_18",
+                    "mass_friction_endstops_19",
+                },
             ),
         ]
 
     def _force_by_mass_alias(
         self,
+        time_s: float,
         pneumatic_snapshot: object,
         mechanical_snapshot: object,
     ) -> dict[str, float]:
@@ -3948,6 +3957,13 @@ class TestMqlFullStateClosure:
             force_by_mass_alias[mass_alias] = (
                 force_by_mass_alias.get(mass_alias, 0.0) - contact_force
             )
+        force_by_mass_alias.update(
+            self._large_mass_force_by_alias(
+                time_s,
+                pneumatic_snapshot,
+                mechanical_snapshot,
+            )
+        )
         return force_by_mass_alias
 
     def _mechanical_node_total_force_by_alias(
@@ -3964,6 +3980,30 @@ class TestMqlFullStateClosure:
             ),
         }
 
+    def _large_mass_force_by_alias(
+        self,
+        time_s: float,
+        pneumatic_snapshot: object,
+        mechanical_snapshot: object,
+    ) -> dict[str, float]:
+        node_force = self._mechanical_node_total_force_by_alias(
+            pneumatic_snapshot,
+            mechanical_snapshot,
+        )
+        forcecon_1 = self.mechanical_closure.assembly.force_connectors["forcecon_1"]
+        forcecon_2 = self.mechanical_closure.assembly.force_connectors["forcecon_2"]
+        signals = self.mechanical_closure.assembly.piecewise_signals
+        return {
+            "mass_friction_endstops_18": (
+                node_force["dynamic_mechanical_node_alternative_2"]
+                - forcecon_2.force_at(time_s, signals)
+            ),
+            "mass_friction_endstops_19": (
+                forcecon_1.force_at(time_s, signals)
+                - node_force["dynamic_mechanical_node_alternative_3"]
+            ),
+        }
+
     def _piston_force_by_mass_alias(self, pneumatic_snapshot: object) -> dict[str, float]:
         from PythonModels.systems.test_mql_pneumatic import pressure_to_amesim_gauge_pa
 
@@ -4342,7 +4382,7 @@ class TestMqlSystem:
             resistance_boundary_temperature_k=resistance_boundary_temperature_k,
         )
         return integrate_ode(
-            rhs=lambda t, state: closure.rhs(state),
+            rhs=lambda t, state: closure.rhs_at(t, state),
             initial_state=closure.initial_state_vector(),
             config=config or SolveIVPConfig(t_stop=1.0e-4, max_step=1.0e-5),
             t_eval=t_eval,

PythonModels/systems/test_mql_mechanical.py

@@ -274,17 +274,28 @@ class TestMqlMechanicalMassClosure:
         state_vector: list[float],
         *,
         force_by_mass_alias: dict[str, float] | None = None,
+        constrained_mass_aliases: set[str] | None = None,
     ) -> list[float]:
         snapshot = self.snapshot(state_vector)
         force_by_mass_alias = force_by_mass_alias or {}
+        constrained_mass_aliases = constrained_mass_aliases or set()
         derivatives: list[float] = []
         for state in snapshot.states:
-            mass = self.assembly.masses[state.alias].endstop()
+            spec = self.assembly.masses[state.alias]
+            mass = spec.endstop()
+            applied_force = force_by_mass_alias.get(state.alias, 0.0)
             acceleration, velocity = mass.derivatives(
                 velocity_m_s=state.velocity_m_s,
                 displacement_m=state.displacement_m,
-                port_1_force_n=force_by_mass_alias.get(state.alias, 0.0),
+                port_1_force_n=applied_force,
             )
+            if state.alias in constrained_mass_aliases and _limit_constraint_holds(
+                spec,
+                state,
+                applied_force,
+            ):
+                acceleration = 0.0
+                velocity = 0.0
             derivatives.extend([acceleration, velocity])
         return derivatives
 
@@ -414,6 +425,20 @@ def _build_mechanical_node(
     )
 
 
+def _limit_constraint_holds(
+    spec: TestMqlMassEndstopSpec,
+    state: TestMqlMechanicalMassState,
+    applied_force_n: float,
+) -> bool:
+    if abs(state.velocity_m_s) > spec.stick_velocity_threshold_m_s:
+        return False
+    at_lower_limit = state.displacement_m <= spec.xmin_m + spec.min_penetration_m
+    at_upper_limit = state.displacement_m >= spec.xmax_m - spec.max_penetration_m
+    return (at_lower_limit and applied_force_n <= 0.0) or (
+        at_upper_limit and applied_force_n >= 0.0
+    )
+
+
 def _build_piecewise_signal(
     component: TestMqlResolvedComponent,
     variable_catalog: TestMqlVariableCatalog | None,

tests/test_test_mql_pnl0001_segment.py

@@ -597,6 +597,41 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
             expected_contact_force,
         )
 
+    def test_full_state_closure_constrains_large_masses_at_active_limits(self) -> None:
+        closure = self.system.full_state_closure_from_spec(
+            self.spec,
+            inlet_node_pressure_pa=15.31e6,
+            resistance_boundary_pressure_pa=15.29e6,
+        )
+
+        rhs = closure.rhs(closure.initial_state_vector())
+
+        self.assertAlmostEqual(rhs[128], 0.0)
+        self.assertAlmostEqual(rhs[129], 0.0)
+        self.assertAlmostEqual(rhs[130], 0.0)
+        self.assertAlmostEqual(rhs[131], 0.0)
+
+    def test_full_state_closure_releases_large_mass_when_node_force_points_inward(self) -> None:
+        closure = self.system.full_state_closure_from_spec(
+            self.spec,
+            inlet_node_pressure_pa=15.31e6,
+            resistance_boundary_pressure_pa=15.29e6,
+        )
+        state = closure.initial_state_vector()
+        state[113] = 1.0e-3
+
+        rhs = closure.rhs_at(0.8000000000000009, state)
+        endstop = closure.mechanical_closure.assembly.elastic_endstops[
+            "elasticendstop_8"
+        ].endstop()
+        expected_contact_force = endstop.contact_force(-1.0, 0.0)
+        support_mass = closure.mechanical_closure.assembly.masses[
+            "mass_friction_endstops_18"
+        ].mass_kg
+
+        self.assertAlmostEqual(rhs[128], expected_contact_force / support_mass)
+        self.assertAlmostEqual(rhs[129], 0.0)
+
     def test_simulates_full_132_state_closure(self) -> None:
         solution = self.system.simulate_full_state_from_spec(
             self.spec,