diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index 75553f5..b35b8c2 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -541,6 +541,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): self.port_2 = self.register_declared_port("port_2") self.port_2.p = self.p0 self.port_2.h_outflow = initial_h + self._connected_h: dict[str, float] = {} @staticmethod def _integer_parameter(name: str, value: float) -> int: @@ -860,7 +861,7 @@ class AmesimPnl0002(AmesimPnl0001): """AMESim PNL0002 R-C-R pneumatic pipe with one center compliance.""" MODEL_TYPE = "amesim_pnl0002" - MODEL_VERSION = "0.2.0" + MODEL_VERSION = "0.3.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), @@ -940,13 +941,43 @@ class AmesimPnl0002(AmesimPnl0001): port_pressure: float, center_pressure: float, center_temperature: float, + *, + port_name: str | None = None, ) -> float: - return self.mass_flow(port_pressure, center_pressure, center_temperature) + upstream_temperature = center_temperature + if ( + port_name is not None + and port_name in self._connected_h + and port_pressure > center_pressure + ): + inlet_h = self._connected_h[port_name] + upstream_temperature = self.medium.temperature_from_pressure_enthalpy( + max(port_pressure, 1.0), + inlet_h, + ) + return self.mass_flow( + port_pressure, + center_pressure, + max(upstream_temperature, 1.0), + ) + + def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: + self._connected_h = dict(connected_h) def component_result_values(self) -> Mapping[str, float]: props = self.properties() - flow_1 = self.port_mass_flow(self.port_1.p, props.p, props.T) - flow_2 = self.port_mass_flow(self.port_2.p, props.p, props.T) + flow_1 = self.port_mass_flow( + self.port_1.p, + props.p, + props.T, + port_name="port_1", + ) + flow_2 = self.port_mass_flow( + self.port_2.p, + props.p, + props.T, + port_name="port_2", + ) diagnostic_flow = flow_1 if abs(flow_1) >= abs(flow_2) else flow_2 upstream_pressure = max(self.port_1.p, self.port_2.p, props.p, 1.0) density = max(self.medium.density(upstream_pressure, props.T), 1.0e-12) @@ -980,7 +1011,12 @@ class AmesimPnl0002(AmesimPnl0001): ), role="flow", value=self.port_1.m_flow - - self.port_mass_flow(self.port_1.p, props.p, props.T), + - self.port_mass_flow( + self.port_1.p, + props.p, + props.T, + port_name="port_1", + ), ), EquationResidual( id=f"{self.name}:port_2_pressure_flow_relation", @@ -994,7 +1030,12 @@ class AmesimPnl0002(AmesimPnl0001): ), role="flow", value=self.port_2.m_flow - - self.port_mass_flow(self.port_2.p, props.p, props.T), + - self.port_mass_flow( + self.port_2.p, + props.p, + props.T, + port_name="port_2", + ), ), ) diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py index dfffb0f..f50ba17 100644 --- a/app/simulation/systems/generic.py +++ b/app/simulation/systems/generic.py @@ -40,6 +40,10 @@ class SimulationPreparationError(ValueError): self.issues = issues +class ThermofluidClosureError(RuntimeError): + """Raised when stream enthalpy and pressure-flow do not reach one fixed point.""" + + @dataclass(frozen=True) class GenericSimulationResult: success: bool @@ -260,6 +264,7 @@ class GenericFluidSystem: self.max_algebraic_residual = 0.0 self.max_algebraic_evaluations = 0 self.max_stream_iterations = 0 + self.max_thermofluid_iterations = 0 self.signal_propagation_count = 0 self.pneumatic_volume_propagation_count = 0 @@ -280,20 +285,62 @@ class GenericFluidSystem: for component in self.dynamic_components: component.refresh_thermodynamic_ports() algebraic = self.pressure_flow_solver.solve() + pressure_flow_solve_count = 1 pneumatic_volume = self.pneumatic_volume_resolver.solve() self.pneumatic_volume_propagation_count += pneumatic_volume.propagated if pneumatic_volume.propagated: for component in self.dynamic_components: component.refresh_thermodynamic_ports() algebraic = self.pressure_flow_solver.solve() - stream, connected_h = self.stream_resolver.solve() + pressure_flow_solve_count += 1 + # Some constitutive flow laws recover their upstream temperature from - # the connected stream enthalpy. Stream propagation updates that - # cache after the first pressure-flow pass, so refresh explicit flows - # once more before evaluating state derivatives and result variables. - algebraic = self.pressure_flow_solver.solve() + # connected stream enthalpy, while junction stream mixing itself depends + # on the resulting mass flows. A single stream -> pressure-flow refresh + # leaves that two-way coupling to the next RHS call, making the ODE RHS + # depend on evaluation history and corrupting finite-difference + # Jacobians. Close both layers to one fixed point inside this call. + physical_ports = tuple( + port + for component in self.network.components.values() + for definition in component.port_definitions + if definition.kind == "physical" + for port in (component.get_port(definition.name),) + ) + connected_h: dict[str, dict[str, float]] = {} + max_coupling_iterations = 25 + flow_relative_tolerance = 1.0e-12 + for coupling_iteration in range(1, max_coupling_iterations + 1): + previous_flows = tuple(port.m_flow for port in physical_ports) + stream, connected_h = self.stream_resolver.solve() + for component in self.dynamic_components: + component.update_stream_outflows(connected_h[component.name]) + algebraic = self.pressure_flow_solver.solve() + pressure_flow_solve_count += 1 + current_flows = tuple(port.m_flow for port in physical_ports) + flow_scale = max( + [abs(value) for value in (*previous_flows, *current_flows)] + [1.0] + ) + max_flow_delta = max( + ( + abs(current - previous) + for previous, current in zip(previous_flows, current_flows) + ), + default=0.0, + ) + if max_flow_delta <= flow_relative_tolerance * flow_scale: + break + else: + raise ThermofluidClosureError( + "Stream enthalpy and pressure-flow coupling did not converge " + f"after {max_coupling_iterations} iterations." + ) + self.max_thermofluid_iterations = max( + self.max_thermofluid_iterations, + coupling_iteration, + ) self.mechanical_state_reducer.update_constraint_accelerations() - self.algebraic_solve_count += 2 + int(bool(pneumatic_volume.propagated)) + self.algebraic_solve_count += pressure_flow_solve_count self.max_algebraic_residual = max( self.max_algebraic_residual, algebraic.max_scaled_residual, @@ -467,6 +514,7 @@ class GenericFluidSystem: }, "stream": { "maxIterationsPerSolve": self.max_stream_iterations, + "maxThermofluidIterations": self.max_thermofluid_iterations, "last": ( self.stream_resolver.last_diagnostics.as_dict() if self.stream_resolver.last_diagnostics is not None diff --git a/tests/test_amesim_pnl0002_pnl0003_component.py b/tests/test_amesim_pnl0002_pnl0003_component.py index c7be92d..98f1a97 100644 --- a/tests/test_amesim_pnl0002_pnl0003_component.py +++ b/tests/test_amesim_pnl0002_pnl0003_component.py @@ -45,6 +45,80 @@ class AmesimPnl0002ComponentTests(unittest.TestCase): self.assertGreater(forward, 0.0) self.assertLess(reverse, 0.0) + def test_external_upstream_flow_uses_connected_stream_temperature(self) -> None: + pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=350.0) + center = pipe.properties() + external_pressure = 15.3e6 + external_temperature = 293.15 + external_h = self.medium.specific_enthalpy_at_pressure( + external_pressure, + external_temperature, + ) + pipe.update_stream_outflows( + {"port_1": external_h, "port_2": center.h} + ) + + flow = pipe.port_mass_flow( + external_pressure, + center.p, + center.T, + port_name="port_1", + ) + expected = pipe.mass_flow( + external_pressure, + center.p, + external_temperature, + ) + center_temperature_flow = pipe.mass_flow( + external_pressure, + center.p, + center.T, + ) + + self.assertAlmostEqual(flow, expected) + self.assertNotAlmostEqual(flow, center_temperature_flow) + + def test_missing_stream_cache_preserves_center_temperature_fallback(self) -> None: + pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=350.0) + center = pipe.properties() + external_pressure = 15.3e6 + + flow = pipe.port_mass_flow( + external_pressure, + center.p, + center.T, + port_name="port_1", + ) + + self.assertAlmostEqual( + flow, + pipe.mass_flow(external_pressure, center.p, center.T), + ) + + def test_center_upstream_flow_keeps_center_temperature(self) -> None: + pipe = AmesimPnl0002("pnl_2", self.medium, p0=15.3e6, T0=350.0) + center = pipe.properties() + external_pressure = 1.0e5 + external_h = self.medium.specific_enthalpy_at_pressure( + external_pressure, + 293.15, + ) + pipe.update_stream_outflows( + {"port_1": external_h, "port_2": center.h} + ) + + flow = pipe.port_mass_flow( + external_pressure, + center.p, + center.T, + port_name="port_1", + ) + + self.assertAlmostEqual( + flow, + pipe.mass_flow(external_pressure, center.p, center.T), + ) + def test_pressure_flow_residuals_use_two_port_resistances(self) -> None: pipe = AmesimPnl0002("pnl_2", self.medium, p0=100000.0, T0=293.15) center = pipe.properties() diff --git a/tests/test_component_catalog.py b/tests/test_component_catalog.py index f96dbb0..1bce93c 100644 --- a/tests/test_component_catalog.py +++ b/tests/test_component_catalog.py @@ -464,7 +464,10 @@ class ComponentCatalogTests(unittest.TestCase): with self.subTest(model_type=model_type): component = self.components[model_type] model_parameters = parameters(model_type) - self.assertEqual(component["modelVersion"], "0.2.0") + expected_version = ( + "0.3.0" if model_type == "amesim_pnl0002" else "0.2.0" + ) + self.assertEqual(component["modelVersion"], expected_version) self.assertEqual(model_parameters["mode"]["editor"], "choice") self.assertEqual(model_parameters["mode"]["options"], thermal_options) self.assertEqual( diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py index 74a72d7..4caff61 100644 --- a/tests/test_pressure_flow_solver_initialization.py +++ b/tests/test_pressure_flow_solver_initialization.py @@ -11,6 +11,8 @@ from app.simulation.components.amesim.flow.orifices import ( ) from app.simulation.components.amesim.flow.pipes import AmesimPnl00r from app.simulation.components.amesim.flow.pipes import AmesimPnl0001 +from app.simulation.components.amesim.flow.pipes import AmesimPnl0002 +from app.simulation.components.amesim.junctions.nodes import AmesimPn3Node2 from app.simulation.components.amesim.media.mediums import ( AmesimHeliumPengRobinsonMedium, ) @@ -142,6 +144,49 @@ class PressureFlowSolverInitializationTests(unittest.TestCase): ) self.assertAlmostEqual(valve.port_3.m_flow, -valve.port_2.m_flow, places=12) + def test_stream_dependent_pipe_flow_closes_with_junction_mixing_in_one_rhs(self) -> None: + medium = IdealGasMedium() + hot = Cylinder("hot", medium, V=0.1, p0=500_000.0, T0=400.0) + cold = Cylinder("cold", medium, V=0.1, p0=500_000.0, T0=250.0) + sink = Tank("sink", medium, V=0.1, p0=100_000.0, T0=300.0) + hot_line = AmesimPnl00r("hot_line", medium, diam=0.01, le=0.5) + cold_line = AmesimPnl00r("cold_line", medium, diam=0.01, le=0.5) + junction = AmesimPn3Node2("junction") + pipe = AmesimPnl0002( + "pipe", + medium, + diam=0.01, + le=1.0, + p0=200_000.0, + T0=300.0, + ) + network = SimulationNetwork("stream-flow-fixed-point") + for component in ( + hot, + cold, + sink, + hot_line, + cold_line, + junction, + pipe, + ): + network.add_component(component) + network.connect("hot", "port_b", "hot_line", "port_1") + network.connect("hot_line", "port_2", "junction", "port_1") + network.connect("cold", "port_b", "cold_line", "port_1") + network.connect("cold_line", "port_2", "junction", "port_3") + network.connect("junction", "port_2", "pipe", "port_1") + network.connect("pipe", "port_2", "sink", "port_a") + + system = GenericFluidSystem(network) + state = system.consistent_initial_state_vector() + first = system.rhs(0.0, state) + second = system.rhs(0.0, state) + + self.assertGreater(system.max_thermofluid_iterations, 1) + for first_value, second_value in zip(first, second): + self.assertAlmostEqual(first_value, second_value, places=10) + def test_current_storage_pressure_reseeds_stale_orifice_ports_and_flow(self) -> None: network, medium, high, low, valve = self._near_equal_pressure_network() solver = PressureFlowSolver(network, max_evaluations=10)