diff --git a/PythonModels/components/amesim_pneumatic.py b/PythonModels/components/amesim_pneumatic.py index 68ad0d0..98cec55 100644 --- a/PythonModels/components/amesim_pneumatic.py +++ b/PythonModels/components/amesim_pneumatic.py @@ -82,12 +82,24 @@ class AmesimPneumaticVolume(DynamicComponent): gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, p0: float = 101_325.0, T0: float = 293.15, + heat_transfer_coefficient: float = 0.0, + heat_transfer_area: float = 0.0, + external_temperature_k: float = 293.15, ) -> None: if volume <= 0.0: raise ValueError("volume must be positive.") + if heat_transfer_coefficient < 0.0: + raise ValueError("heat_transfer_coefficient must be non-negative.") + if heat_transfer_area < 0.0: + raise ValueError("heat_transfer_area must be non-negative.") + if external_temperature_k <= 0.0: + raise ValueError("external_temperature_k must be positive.") super().__init__(name=name) self.volume = volume self.gas = gas + self.heat_transfer_coefficient = heat_transfer_coefficient + self.heat_transfer_area = heat_transfer_area + self.external_temperature = external_temperature_k rho0 = gas.density(p0, T0) m0 = rho0 * volume U0 = m0 * gas.specific_internal_energy(T0) @@ -103,8 +115,20 @@ class AmesimPneumaticVolume(DynamicComponent): gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, p0: float = 101_325.0, T0: float = 293.15, + heat_transfer_coefficient: float = 0.0, + heat_transfer_area: float = 0.0, + external_temperature_k: float = 293.15, ) -> "AmesimPneumaticVolume": - return cls(name=name, volume=liters_to_m3(volume_liters), gas=gas, p0=p0, T0=T0) + return cls( + name=name, + volume=liters_to_m3(volume_liters), + gas=gas, + p0=p0, + T0=T0, + heat_transfer_coefficient=heat_transfer_coefficient, + heat_transfer_area=heat_transfer_area, + external_temperature_k=external_temperature_k, + ) def get_state_vector(self) -> list[float]: return self.state.as_vector() @@ -118,6 +142,14 @@ class AmesimPneumaticVolume(DynamicComponent): def volume_rate_m3_s(self) -> float: return 0.0 + def thermal_energy_flow_w(self, temperature_k: float | None = None) -> float: + temperature = self.properties().T if temperature_k is None else temperature_k + return ( + self.heat_transfer_coefficient + * self.heat_transfer_area + * (self.external_temperature - temperature) + ) + def gas_mass_g(self) -> float: return kg_to_g(self.state.m) @@ -139,7 +171,10 @@ class AmesimPneumaticVolume(DynamicComponent): return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h) def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: - return VolumeState(m=m_flow, U=m_flow * inlet_h) + return VolumeState( + m=m_flow, + U=m_flow * inlet_h + self.thermal_energy_flow_w(), + ) def derivatives_from_two_connections( self, @@ -151,6 +186,7 @@ class AmesimPneumaticVolume(DynamicComponent): internal_h: float, volume_rate_m3_s: float | None = None, ) -> VolumeState: + properties = self.properties() inlet_h_a = self.connection_inlet_enthalpy( port_m_flow=port_a_m_flow, connected_h=connected_h_a, @@ -166,7 +202,8 @@ class AmesimPneumaticVolume(DynamicComponent): U=( port_a_m_flow * inlet_h_a + port_b_m_flow * inlet_h_b - - self.properties().p * ( + + self.thermal_energy_flow_w(properties.T) + - properties.p * ( self.volume_rate_m3_s() if volume_rate_m3_s is None else volume_rate_m3_s @@ -186,6 +223,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume): p0: float = 101_325.0, T0: float = 293.15, external_volume: float = 0.0, + heat_transfer_coefficient: float = 0.0, + heat_transfer_area: float = 0.0, + external_temperature_k: float = 293.15, ) -> None: if dead_volume <= 0.0: raise ValueError("dead_volume must be positive.") @@ -200,6 +240,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume): gas=gas, p0=p0, T0=T0, + heat_transfer_coefficient=heat_transfer_coefficient, + heat_transfer_area=heat_transfer_area, + external_temperature_k=external_temperature_k, ) @classmethod @@ -211,6 +254,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume): p0: float = 101_325.0, T0: float = 293.15, external_volume_liters: float = 0.0, + heat_transfer_coefficient: float = 0.0, + heat_transfer_area: float = 0.0, + external_temperature_k: float = 293.15, ) -> "AmesimVariablePneumaticVolume": return cls( name=name, @@ -219,6 +265,9 @@ class AmesimVariablePneumaticVolume(AmesimPneumaticVolume): p0=p0, T0=T0, external_volume=liters_to_m3(external_volume_liters), + heat_transfer_coefficient=heat_transfer_coefficient, + heat_transfer_area=heat_transfer_area, + external_temperature_k=external_temperature_k, ) def volume_rate_m3_s(self) -> float: diff --git a/PythonModels/components/amesim_pneumatic_line.py b/PythonModels/components/amesim_pneumatic_line.py index 7ce54f8..b4b1d37 100644 --- a/PythonModels/components/amesim_pneumatic_line.py +++ b/PythonModels/components/amesim_pneumatic_line.py @@ -48,7 +48,7 @@ class _DarcyPipeResistanceMixin: if upper > 1.0e3: raise ValueError("unable to bracket pneumatic pipe resistance flow") lower = 0.0 - for _ in range(80): + for _ in range(48): middle = 0.5 * (lower + upper) if self._darcy_pressure_drop( middle, @@ -255,15 +255,19 @@ class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent): connected_h_2: float, ) -> VolumeState: internal = self.properties() - inlet_h_1 = self.connection_inlet_enthalpy( - port_m_flow=port_1_m_flow, - connected_h=connected_h_1, - internal_h=internal.h, + # PNL0001 is a fixed-volume distributed line store. Its transported + # energy variable therefore follows specific internal energy, not the + # chamber-style stagnation enthalpy contract. For this ideal gas, + # h = gamma * u. Outflow always carries the local u. + inlet_u_1 = ( + connected_h_1 / self.gas.gamma + if port_1_m_flow > 0.0 + else internal.u ) - inlet_h_2 = self.connection_inlet_enthalpy( - port_m_flow=port_2_m_flow, - connected_h=connected_h_2, - internal_h=internal.h, + inlet_u_2 = ( + connected_h_2 / self.gas.gamma + if port_2_m_flow > 0.0 + else internal.u ) heat_flow = ( self.heat_transfer_coefficient @@ -272,7 +276,7 @@ class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent): ) return VolumeState( m=port_1_m_flow + port_2_m_flow, - U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow, + U=port_1_m_flow * inlet_u_1 + port_2_m_flow * inlet_u_2 + heat_flow, ) diff --git a/PythonModels/core/solver.py b/PythonModels/core/solver.py index 3e101e1..ce755d9 100644 --- a/PythonModels/core/solver.py +++ b/PythonModels/core/solver.py @@ -10,8 +10,9 @@ class SolveIVPConfig: t_stop: float = 20.0 method: str = "BDF" rtol: float = 1e-6 - atol: float = 1e-8 + atol: float = 1e-10 max_step: float = 1e-3 + first_step: float | None = None @dataclass(frozen=True) @@ -91,12 +92,16 @@ def integrate_ode( except ImportError: return _runge_kutta_4(rhs, initial_state, config, t_eval) - return solve_ivp( - fun=rhs, - t_span=(config.t_start, config.t_stop), - y0=initial_state, - method=config.method, - rtol=config.rtol, - atol=config.atol, - t_eval=t_eval, - ) + solve_options = { + "fun": rhs, + "t_span": (config.t_start, config.t_stop), + "y0": initial_state, + "method": config.method, + "rtol": config.rtol, + "atol": config.atol, + "max_step": config.max_step, + "t_eval": t_eval, + } + if config.first_step is not None: + solve_options["first_step"] = config.first_step + return solve_ivp(**solve_options) diff --git a/PythonModels/scripts/run_test_mql_full_state_comparison.py b/PythonModels/scripts/run_test_mql_full_state_comparison.py index 131a412..68c17a7 100644 --- a/PythonModels/scripts/run_test_mql_full_state_comparison.py +++ b/PythonModels/scripts/run_test_mql_full_state_comparison.py @@ -559,6 +559,8 @@ def format_test_mql_full_state_comparison_summary( f" - mass_derivative_kg_s={chamber_diagnostic.mass_derivative_kg_s}", f" - port_a_energy_flow_w={chamber_diagnostic.port_a_energy_flow_w}", f" - boundary_work_w={chamber_diagnostic.boundary_work_w}", + f" - thermal_energy_flow_w=" + f"{chamber_diagnostic.thermal_energy_flow_w}", f" - energy_derivative_w={chamber_diagnostic.energy_derivative_w}", ] ) diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 55aabcd..9a19ba9 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3870,6 +3870,7 @@ class TestMqlVariableChamberRhsDiagnostic: port_a_energy_flow_w: float port_b_energy_flow_w: float boundary_work_w: float + thermal_energy_flow_w: float energy_derivative_w: float @@ -4120,6 +4121,9 @@ class TestMqlFullStateClosure: port_a_energy_flow = port_a_mass_flow * port_a_inlet_h port_b_energy_flow = port_b_mass_flow * port_b_inlet_h boundary_work = -chamber_properties.p * chamber.volume_rate_m3_s() + thermal_energy_flow = chamber.thermal_energy_flow_w( + chamber_properties.T + ) return TestMqlVariableChamberRhsDiagnostic( chamber_alias=chamber_alias, piston_alias=piston_alias, @@ -4133,8 +4137,12 @@ class TestMqlFullStateClosure: port_a_energy_flow_w=port_a_energy_flow, port_b_energy_flow_w=port_b_energy_flow, boundary_work_w=boundary_work, + thermal_energy_flow_w=thermal_energy_flow, energy_derivative_w=( - port_a_energy_flow + port_b_energy_flow + boundary_work + port_a_energy_flow + + port_b_energy_flow + + boundary_work + + thermal_energy_flow ), ) diff --git a/PythonModels/systems/test_mql_pneumatic.py b/PythonModels/systems/test_mql_pneumatic.py index a203989..da29f5b 100644 --- a/PythonModels/systems/test_mql_pneumatic.py +++ b/PythonModels/systems/test_mql_pneumatic.py @@ -229,6 +229,9 @@ def _build_chamber( gas=gas, p0=initial_pressure_pa, T0=_component_temperature(component), + heat_transfer_coefficient=component.parameter_value("kth"), + heat_transfer_area=component.parameter_value("sth"), + external_temperature_k=_component_temperature(component), ) return AmesimPneumaticVolume.from_liters( name=component.alias, @@ -236,6 +239,9 @@ def _build_chamber( gas=gas, p0=initial_pressure_pa, T0=_component_temperature(component), + heat_transfer_coefficient=component.parameter_value("kth"), + heat_transfer_area=component.parameter_value("sth"), + external_temperature_k=_component_temperature(component), ) diff --git a/tests/test_amesim_pneumatic_components.py b/tests/test_amesim_pneumatic_components.py index 97f62cb..1238be7 100644 --- a/tests/test_amesim_pneumatic_components.py +++ b/tests/test_amesim_pneumatic_components.py @@ -79,6 +79,31 @@ class AmesimPneumaticComponentsTest(unittest.TestCase): ) self.assertAlmostEqual(volume.port_a.p, volume.port_b.p) + def test_volume_applies_amesim_heat_exchange_to_energy_derivative(self) -> None: + volume = AmesimPneumaticVolume.from_liters( + name="heated_chamber", + volume_liters=15.0, + p0=100000.0, + T0=300.0, + heat_transfer_coefficient=1500.0, + heat_transfer_area=0.7, + external_temperature_k=293.15, + ) + props = volume.properties() + + 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.U, + 1500.0 * 0.7 * (293.15 - 300.0), + ) + def test_variable_volume_tracks_external_volume(self) -> None: volume = AmesimVariablePneumaticVolume.from_liters( name="pn_c1_8", diff --git a/tests/test_amesim_pnl0001_pipe.py b/tests/test_amesim_pnl0001_pipe.py index 2c6e78f..fc326a5 100644 --- a/tests/test_amesim_pnl0001_pipe.py +++ b/tests/test_amesim_pnl0001_pipe.py @@ -63,7 +63,8 @@ class AmesimPnl0001PipeTests(unittest.TestCase): self.assertAlmostEqual(derivative.m, 0.1) self.assertAlmostEqual( derivative.U, - 0.2 * (internal.h + 1000.0) - 0.1 * internal.h, + 0.2 * (internal.h + 1000.0) / self.pipe.gas.gamma + - 0.1 * internal.u, ) diff --git a/tests/test_core_solver.py b/tests/test_core_solver.py new file mode 100644 index 0000000..ede7faa --- /dev/null +++ b/tests/test_core_solver.py @@ -0,0 +1,68 @@ +import sys +import types +import unittest +from unittest.mock import patch + +from PythonModels.core.solver import SolveIVPConfig, integrate_ode + + +class IntegrateOdeTests(unittest.TestCase): + def test_scipy_solver_receives_step_size_controls(self) -> None: + calls: list[dict[str, object]] = [] + + def fake_solve_ivp(**kwargs): + calls.append(kwargs) + return object() + + scipy_module = types.ModuleType("scipy") + integrate_module = types.ModuleType("scipy.integrate") + integrate_module.solve_ivp = fake_solve_ivp + scipy_module.integrate = integrate_module + + with patch.dict( + sys.modules, + {"scipy": scipy_module, "scipy.integrate": integrate_module}, + ): + integrate_ode( + rhs=lambda _time, state: state, + initial_state=[1.0], + config=SolveIVPConfig( + t_start=0.0, + t_stop=1.0, + max_step=1.0e-4, + first_step=1.0e-8, + ), + t_eval=[0.0, 1.0], + ) + + self.assertEqual(calls[0]["max_step"], 1.0e-4) + self.assertEqual(calls[0]["first_step"], 1.0e-8) + + def test_scipy_solver_omits_unset_first_step(self) -> None: + calls: list[dict[str, object]] = [] + + def fake_solve_ivp(**kwargs): + calls.append(kwargs) + return object() + + scipy_module = types.ModuleType("scipy") + integrate_module = types.ModuleType("scipy.integrate") + integrate_module.solve_ivp = fake_solve_ivp + scipy_module.integrate = integrate_module + + with patch.dict( + sys.modules, + {"scipy": scipy_module, "scipy.integrate": integrate_module}, + ): + integrate_ode( + rhs=lambda _time, state: state, + initial_state=[1.0], + config=SolveIVPConfig(t_start=0.0, t_stop=1.0), + ) + + self.assertEqual(calls[0]["max_step"], 1.0e-3) + self.assertNotIn("first_step", calls[0]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_run_test_mql_full_state_comparison.py b/tests/test_run_test_mql_full_state_comparison.py index e04aa2f..6fecd5f 100644 --- a/tests/test_run_test_mql_full_state_comparison.py +++ b/tests/test_run_test_mql_full_state_comparison.py @@ -89,7 +89,8 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): chamber_diagnostic.energy_derivative_w, chamber_diagnostic.port_a_energy_flow_w + chamber_diagnostic.port_b_energy_flow_w - + chamber_diagnostic.boundary_work_w, + + chamber_diagnostic.boundary_work_w + + chamber_diagnostic.thermal_energy_flow_w, ) self.assertTrue(csv_summary_exists) self.assertIn("Mode: Python 132 full-state closure comparison", summary) @@ -112,6 +113,7 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): self.assertIn("piston_alias=pn_brp2_8", summary) self.assertIn("mass_derivative_kg_s=", summary) self.assertIn("boundary_work_w=", summary) + self.assertIn("thermal_energy_flow_w=", summary) self.assertIn("final_python=", summary) self.assertIn("final_amesim=", summary) self.assertIn("python.press@pn_c1_8", csv_header) diff --git a/tests/test_test_mql_pneumatic.py b/tests/test_test_mql_pneumatic.py index 226ff1d..22fc84a 100644 --- a/tests/test_test_mql_pneumatic.py +++ b/tests/test_test_mql_pneumatic.py @@ -55,6 +55,9 @@ class TestMqlPneumaticAssemblyTests(unittest.TestCase): self.assertAlmostEqual(fixed_chamber.volume, 0.057) self.assertAlmostEqual(variable_chamber.dead_volume, 0.015) self.assertAlmostEqual(variable_chamber.volume, 0.015) + self.assertAlmostEqual(variable_chamber.heat_transfer_coefficient, 1500.0) + self.assertAlmostEqual(variable_chamber.heat_transfer_area, 0.7) + self.assertAlmostEqual(variable_chamber.external_temperature, 293.15) self.assertIs(fixed_chamber.gas, HELIUM_PNEUMATIC_GAS) self.assertIs(variable_chamber.gas, HELIUM_PNEUMATIC_GAS) self.assertAlmostEqual(fixed_chamber.properties().T, 293.15) diff --git a/tests/test_test_mql_pnl0001_segment.py b/tests/test_test_mql_pnl0001_segment.py index bf06d55..680d619 100644 --- a/tests/test_test_mql_pnl0001_segment.py +++ b/tests/test_test_mql_pnl0001_segment.py @@ -709,7 +709,8 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase): diagnostic.energy_derivative_w, diagnostic.port_a_energy_flow_w + diagnostic.port_b_energy_flow_w - + diagnostic.boundary_work_w, + + diagnostic.boundary_work_w + + diagnostic.thermal_energy_flow_w, ) def test_full_state_closure_applies_step_controls_to_variable_orifices(self) -> None: @@ -731,6 +732,27 @@ class TestMqlPn3P4NodeChamberSegmentTests(unittest.TestCase): open_snapshot.p4_port3_remote_node_to_primary_line_flow, 0.0, ) + node_inlet_h = ( + open_snapshot.p4_port3_remote_balance.port_2_enthalpy_flow_w + / open_snapshot.p4_port3_remote_node_to_primary_line_flow + ) + self.assertAlmostEqual( + node_inlet_h, + open_snapshot.p4_port3_remote_orifice_line_port_2.h, + ) + rhs = closure.rhs_at(0.04, state) + line = closure.pneumatic_closure.components.p4_port3_remote_primary_line + expected_line_derivative = line.derivatives_from_connections( + port_1_m_flow=( + open_snapshot.p4_port3_remote_chamber_to_line_flow + ), + connected_h_1=open_snapshot.p4_port3_remote_primary_chamber.h, + port_2_m_flow=( + open_snapshot.p4_port3_remote_node_to_primary_line_flow + ), + connected_h_2=node_inlet_h, + ) + self.assertAlmostEqual(rhs[45], expected_line_derivative.U) self.assertAlmostEqual( self.system.pneumatic_assembly.variable_orifices[ "pn_morifice_1"