PythonModels/components/amesim_pneumatic.py

@@ -115,6 +115,9 @@ class AmesimPneumaticVolume(DynamicComponent):
     def volume_cm3(self) -> float:
         return m3_to_cm3(self.volume)
 
+    def volume_rate_m3_s(self) -> float:
+        return 0.0
+
     def gas_mass_g(self) -> float:
         return kg_to_g(self.state.m)
 
@@ -146,6 +149,7 @@ class AmesimPneumaticVolume(DynamicComponent):
         port_b_m_flow: float,
         connected_h_b: float,
         internal_h: float,
+        volume_rate_m3_s: float | None = None,
     ) -> VolumeState:
         inlet_h_a = self.connection_inlet_enthalpy(
             port_m_flow=port_a_m_flow,
@@ -159,7 +163,15 @@ class AmesimPneumaticVolume(DynamicComponent):
         )
         return VolumeState(
             m=port_a_m_flow + port_b_m_flow,
-            U=port_a_m_flow * inlet_h_a + port_b_m_flow * inlet_h_b,
+            U=(
+                port_a_m_flow * inlet_h_a
+                + port_b_m_flow * inlet_h_b
+                - self.properties().p * (
+                    self.volume_rate_m3_s()
+                    if volume_rate_m3_s is None
+                    else volume_rate_m3_s
+                )
+            ),
         )
 
 
@@ -177,10 +189,11 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
     ) -> None:
         if dead_volume <= 0.0:
             raise ValueError("dead_volume must be positive.")
-        if external_volume < 0.0:
+        if dead_volume + external_volume <= 0.0:
-            raise ValueError("external_volume must be non-negative.")
+            raise ValueError("total volume must be positive.")
         self.dead_volume = dead_volume
         self.external_volume = external_volume
+        self.external_volume_rate = 0.0
         super().__init__(
             name=name,
             volume=dead_volume + external_volume,
@@ -208,10 +221,18 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
             external_volume=liters_to_m3(external_volume_liters),
         )
 
-    def set_external_volume_m3(self, external_volume: float) -> None:
+    def volume_rate_m3_s(self) -> float:
-        if external_volume < 0.0:
+        return self.external_volume_rate
-            raise ValueError("external_volume must be non-negative.")
+
+    def set_external_volume_m3(
+        self,
+        external_volume: float,
+        external_volume_rate_m3_s: float = 0.0,
+    ) -> None:
+        if self.dead_volume + external_volume <= 0.0:
+            raise ValueError("total volume must be positive.")
         self.external_volume = external_volume
+        self.external_volume_rate = external_volume_rate_m3_s
         self.volume = self.dead_volume + self.external_volume
 
 

PythonModels/systems/test_mql.py

@@ -3907,8 +3907,9 @@ class TestMqlFullStateClosure:
     def snapshot(self, state_vector: list[float] | None = None) -> TestMqlFullStateSnapshot:
         values = self.initial_state_vector() if state_vector is None else list(state_vector)
         pneumatic_state, mechanical_state = self._split_state(values)
-        pneumatic_snapshot = self.pneumatic_closure.snapshot(pneumatic_state)
         mechanical_snapshot = self.mechanical_closure.snapshot(mechanical_state)
+        self._sync_variable_chamber_kinematics(mechanical_snapshot)
+        pneumatic_snapshot = self.pneumatic_closure.snapshot(pneumatic_state)
         return TestMqlFullStateSnapshot(
             pneumatic=pneumatic_snapshot,
             mechanical=mechanical_snapshot,
@@ -3925,8 +3926,9 @@ class TestMqlFullStateClosure:
 
     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)
+        self._sync_variable_chamber_kinematics(mechanical_snapshot)
+        pneumatic_snapshot = self.pneumatic_closure.snapshot(pneumatic_state)
         return [
             *self.pneumatic_closure.rhs(pneumatic_state),
             *self.mechanical_closure.rhs(
@@ -3943,6 +3945,22 @@ class TestMqlFullStateClosure:
             ),
         ]
 
+    def _sync_variable_chamber_kinematics(self, mechanical_snapshot: object) -> None:
+        for _mass_alias, piston_alias, chamber_field in _PISTON_FORCE_BINDINGS:
+            piston = self.mechanical_closure.assembly.pistons[piston_alias]
+            kinematics = mechanical_snapshot.piston_kinematics_by_alias[piston_alias]
+            chamber = getattr(self.pneumatic_closure.components, chamber_field)
+            chamber.set_external_volume_m3(
+                piston.geometry().chamber_volume_m3(
+                    kinematics.port_3_displacement_m,
+                    kinematics.port_2_displacement_m,
+                ),
+                piston.geometry().chamber_volume_rate_m3_s(
+                    kinematics.port_3_velocity_m_s,
+                    kinematics.port_2_velocity_m_s,
+                ),
+            )
+
     def _force_by_mass_alias(
         self,
         time_s: float,

tests/test_amesim_pneumatic_components.py

@@ -91,6 +91,27 @@ class AmesimPneumaticComponentsTest(unittest.TestCase):
         volume.set_external_volume_m3(cm3_to_m3(34242.54636512914))
         self.assertAlmostEqual(volume.volume_cm3(), 49242.54636512914)
 
+    def test_variable_volume_boundary_work_uses_external_volume_rate(self) -> None:
+        volume = AmesimVariablePneumaticVolume.from_liters(
+            name="pn_c1_8",
+            dead_volume_liters=15.0,
+            p0=100000.0,
+            T0=293.15,
+        )
+        props = volume.properties()
+        volume.set_external_volume_m3(0.0, external_volume_rate_m3_s=2.0e-6)
+
+        derivative = volume.derivatives_from_two_connections(
+            port_a_m_flow=0.0,
+            connected_h_a=props.h,
+            port_b_m_flow=0.0,
+            connected_h_b=props.h,
+            internal_h=props.h,
+        )
+
+        self.assertAlmostEqual(derivative.m, 0.0)
+        self.assertAlmostEqual(derivative.U, -props.p * 2.0e-6)
+
     def test_orifice_returns_signed_mass_flow(self) -> None:
         orifice = AmesimPneumaticOrifice.from_mm2(
             name="pn_orifice_18",

tests/test_test_mql_pnl0001_segment.py

@@ -566,6 +566,43 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
         )
         self.assertAlmostEqual(rhs[112], -0.8167936695810917)
 
+    def test_full_state_closure_feeds_piston_kinematics_to_variable_chambers(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] = 0.05
+
+        snapshot = closure.snapshot(state)
+        chamber = closure.pneumatic_closure.components.p4_port3_remote_primary_chamber
+        piston = closure.mechanical_closure.assembly.pistons["pn_brp2_8"]
+        expected_external_volume = piston.geometry().chamber_volume_m3(0.0, 0.05)
+
+        self.assertAlmostEqual(
+            chamber.volume,
+            chamber.dead_volume + expected_external_volume,
+        )
+        self.assertLess(
+            snapshot.pneumatic.p4_port3_remote_primary_chamber.p,
+            self.system.pneumatic_assembly.variable_initial_absolute_pressure_pa,
+        )
+
+        base_rhs = closure.rhs(closure.initial_state_vector())
+        moving_state = closure.initial_state_vector()
+        moving_state[112] = 0.1
+        moving_snapshot = closure.snapshot(moving_state)
+        rhs = closure.rhs(moving_state)
+        expected_volume_rate = piston.geometry().chamber_volume_rate_m3_s(0.0, 0.1)
+
+        self.assertAlmostEqual(chamber.volume_rate_m3_s(), expected_volume_rate)
+        self.assertAlmostEqual(
+            rhs[47] - base_rhs[47],
+            -moving_snapshot.pneumatic.p4_port3_remote_primary_chamber.p
+            * expected_volume_rate,
+        )
+
     def test_full_state_closure_applies_elastic_endstop_force_to_small_masses(self) -> None:
         closure = self.system.full_state_closure_from_spec(
             self.spec,
@@ -577,7 +614,6 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
 
         snapshot = closure.snapshot(state)
         rhs = closure.rhs(state)
-        initial_rhs = closure.rhs(closure.initial_state_vector())
         endstop = closure.mechanical_closure.assembly.elastic_endstops[
             "elasticendstop_8"
         ].endstop()
@@ -585,10 +621,13 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase):
         mass = closure.mechanical_closure.assembly.masses[
             "mass_friction_endstops_10"
         ].mass_kg
+        expected_piston_force = closure._piston_force_by_mass_alias(
+            snapshot.pneumatic
+        )["mass_friction_endstops_10"]
 
         self.assertAlmostEqual(
             rhs[112],
-            initial_rhs[112] - expected_contact_force / mass,
+            (expected_piston_force - expected_contact_force) / mass,
         )
         self.assertAlmostEqual(
             snapshot.mechanical_node_total_force_by_alias[