From 0fa166c8e5c0a6000d7f0d299e0bc8d97eef2b3c Mon Sep 17 00:00:00 2001 From: huojiarong Date: Mon, 17 Aug 2026 09:23:40 +0000 Subject: [PATCH] Fix pneumatic node zero-flow reversal --- .../components/amesim/junctions/nodes.py | 43 ++++++++- .../test_amesim_pneumatic_node_components.py | 31 +++++++ tests/test_mechanical_solver_causalization.py | 88 ++++++++++++++++++- 3 files changed, 159 insertions(+), 3 deletions(-) diff --git a/app/simulation/components/amesim/junctions/nodes.py b/app/simulation/components/amesim/junctions/nodes.py index be21525..88a5d34 100644 --- a/app/simulation/components/amesim/junctions/nodes.py +++ b/app/simulation/components/amesim/junctions/nodes.py @@ -9,6 +9,24 @@ from app.simulation.core.medium import IdealGasMedium from app.simulation.core.ports import PortDefinition, PortState +_REFERENCE_OUTFLOW_REGULARIZATION_RATIO = 0.05 + + +def _regularized_inverse_outflow(flow: float, transition_flow: float) -> float: + """Return a C1 inverse that tends to zero as a negative flow vanishes.""" + + if flow >= 0.0: + return 0.0 + transition_flow = max(float(transition_flow), 1.0e-12) + if -flow >= transition_flow: + return 1.0 / flow + return ( + flow + * (2.0 * transition_flow * transition_flow - flow * flow) + / transition_flow**4 + ) + + class _AmesimPneumaticNode(AlgebraicComponent): """Shared implementation for AMESim pneumatic junction submodels. @@ -102,7 +120,7 @@ class _AmesimPneumaticNode(AlgebraicComponent): else self.temperature_reference_h ) - if reference_port.m_flow < -1e-12: + if reference_port.m_flow < 0.0: energy_without_reference = sum( port.m_flow * ( @@ -113,8 +131,29 @@ class _AmesimPneumaticNode(AlgebraicComponent): for name, port in self.ports.items() if name != self.REFERENCE_PORT ) + non_reference_flow_scale = sum( + abs(port.m_flow) + for name, port in self.ports.items() + if name != self.REFERENCE_PORT + ) + transition_flow = ( + _REFERENCE_OUTFLOW_REGULARIZATION_RATIO + * non_reference_flow_scale + ) + # Port 2 carries AMESim's residual-energy causality. Exact + # division is singular when its outflow reverses through zero, so + # use a C1 band that matches the exact balance at its boundary and + # tends to the mixed enthalpy at zero flow. + inverse_flow = _regularized_inverse_outflow( + reference_port.m_flow, + transition_flow, + ) + energy_residual_at_mixed_h = ( + energy_without_reference + + reference_port.m_flow * mixed_h + ) reference_port.h_outflow = ( - -energy_without_reference / reference_port.m_flow + mixed_h - energy_residual_at_mixed_h * inverse_flow ) diff --git a/tests/test_amesim_pneumatic_node_components.py b/tests/test_amesim_pneumatic_node_components.py index b2d43d7..4a9b04e 100644 --- a/tests/test_amesim_pneumatic_node_components.py +++ b/tests/test_amesim_pneumatic_node_components.py @@ -123,6 +123,37 @@ class AmesimPneumaticNodeComponentTests(unittest.TestCase): ) self.assertAlmostEqual(energy_flow, 0.0, places=10) + def test_port_2_outflow_enthalpy_remains_finite_through_flow_reversal(self) -> None: + node = AmesimP4Node2("p4_1") + connected_h = { + "port_1": 300.0, + "port_2": 300.0, + "port_3": 200.0, + "port_4": 100.0, + } + + outflow_enthalpies: list[float] = [] + for reference_flow in (-1.0e-6, 0.0, 1.0e-6): + node.port_1.m_flow = -1.0 + node.port_2.m_flow = reference_flow + node.port_3.m_flow = -reference_flow + node.port_4.m_flow = 1.0 + + node.update_stream_outflows(connected_h) + + self.assertAlmostEqual( + sum(port.m_flow for port in node.ports.values()), + 0.0, + places=15, + ) + outflow_enthalpies.append(node.port_2.h_outflow) + + self.assertTrue(all(abs(value) < 1.0e4 for value in outflow_enthalpies)) + self.assertLess( + max(outflow_enthalpies) - min(outflow_enthalpies), + 1.0, + ) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_mechanical_solver_causalization.py b/tests/test_mechanical_solver_causalization.py index 73ec7d3..91d69c5 100644 --- a/tests/test_mechanical_solver_causalization.py +++ b/tests/test_mechanical_solver_causalization.py @@ -1,6 +1,6 @@ from __future__ import annotations -from math import exp +from math import exp, isfinite import unittest from app.simulation.components.amesim.mechanical.translational import ( @@ -130,7 +130,93 @@ def _single_mass_system( return GenericFluidSystem(network), mass +def _high_stiffness_contact_impact_system() -> GenericFluidSystem: + """Return the reduced mechanical core of the four-branch stop impact. + + The production model has four 50 kg piston groups coupled through LSTP00A + contacts to one 170000 kg support group. When the support reaches its + ideal upper stop, its velocity is reset while each contact still carries + the finite, extremely fast LSTP damping transient. One branch is enough + to retain that stiffness and state-transition interaction in a fast test. + """ + + medium = IdealGasMedium() + zero_left = AmesimF000("zero_left") + moving_mass = AmesimMecmas21( + "moving_mass", + medium, + mass=50.0, + useFriction=1.0, + stoptype=4.0, + x0=0.0090025, + v0=1.0, + ) + contact = AmesimLstp00a( + "contact", + medium, + gap0=0.0, + kcont=1.0e11, + rcont=1.0e11, + Pdis=1.0e-7, + discContactOption=1.0, + ) + support_mass = AmesimMecmas21( + "support_mass", + medium, + mass=170000.0, + useFriction=1.0, + stoptype=1.0, + xmin=0.0, + xmax=0.01, + x0=0.009, + v0=1.0, + ) + zero_right = AmesimF000("zero_right") + + network = SimulationNetwork("high-stiffness-contact-impact") + for component in ( + zero_left, + moving_mass, + contact, + support_mass, + zero_right, + ): + network.add_component(component) + network.connect("zero_left", "port_1", "moving_mass", "port_1") + network.connect("moving_mass", "port_2", "contact", "port_1") + network.connect("contact", "port_2", "support_mass", "port_1") + network.connect("support_mass", "port_2", "zero_right", "port_1") + return GenericFluidSystem(network) + + class MechanicalSolverCausalizationTests(unittest.TestCase): + def test_high_stiffness_lstp_survives_ideal_stop_transition(self) -> None: + system = _high_stiffness_contact_impact_system() + + result = system.simulate( + SolveIVPConfig( + t_start=0.0, + t_stop=0.002, + max_step=0.002, + method="BDF", + rtol=1.0e-6, + ), + sample_step=0.00025, + ) + + self.assertTrue(result.success, result.message) + totals = result.diagnostics["integration"]["totals"] + self.assertEqual(totals["stateTransitionCount"], 1) + self.assertEqual(totals["solverStartCount"], 2) + self.assertEqual(totals["recoverableRetryCount"], 0) + + self.assertAlmostEqual(result.series["support_mass.x"][-1], 0.01) + self.assertAlmostEqual(result.series["support_mass.v"][-1], 0.0, delta=1.0e-9) + self.assertAlmostEqual(result.series["moving_mass.v"][-1], 0.0, delta=1.0e-4) + self.assertGreater(max(result.series["contact.force"]), 1.0e10) + for values in result.series.values(): + self.assertTrue(all(isfinite(value) for value in values)) + def test_lstp_contact_uses_exponential_damping_ramp_and_negative_force_option( self, ) -> None: