diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py index 01001fd..bde49fc 100644 --- a/app/simulation/components/amesim/flow/orifices.py +++ b/app/simulation/components/amesim/flow/orifices.py @@ -35,6 +35,7 @@ _FLOWSET_USES_CV = (ParameterCondition("flowset", (2.0,)),) _FLOWSET_USES_KV = (ParameterCondition("flowset", (3.0,)),) _PN_PRESSURE_RATIO_ACCURACY = 0.9999 _PN_LAMINAR_SMOOTHING_GAIN = 12.0 +_PNVO001_CLOSED_OPENING_ABS_TOL = 1.0e-12 _PNOR001_FLOW_COEFFICIENT_GROUP = ParameterGroupDisplaySpec( id="flow_coefficient", label="流量系数", @@ -835,6 +836,16 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): downstream_pressure=downstream_pressure, upstream_temperature=upstream_temperature, ) + # AMESim reports no vena-contracta velocity while the valve is closed. + # Signal propagation around a step can leave a round-off-sized opening, + # so apply the same numerical-zero convention to this diagnostic only. + if isclose( + self.opening, + 0.0, + rel_tol=0.0, + abs_tol=_PNVO001_CLOSED_OPENING_ABS_TOL, + ): + gas_velocity = 0.0 return { "xv": self.opening, "cm": mass_flow_parameter, diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index d546a67..12d4b5b 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -75,6 +75,31 @@ def _reported_friction_factor(value: float) -> float: return min(float(value), _MAX_REPORTED_FRICTION_FACTOR) +def _bare_pipe_mass_flow_parameter( + mass_flow: float, + *, + area: float, + diameter: float, + upstream_pressure: float, + upstream_temperature: float, + resistance_length: float, + friction_factor: float, +) -> float: + """Return Amesim's ``Cm`` without the pipe flow coefficient ``Cq``.""" + + flow_coefficient = sqrt( + diameter / (resistance_length * friction_factor) + ) + return ( + abs(mass_flow) + * sqrt(upstream_temperature) + / max( + flow_coefficient * area * upstream_pressure, + 1.0e-18, + ) + ) + + _DYNAMIC_PIPE_POLYTROPIC_MODE = ParameterCondition("mode", (1.0,)) _DYNAMIC_PIPE_HEAT_EXCHANGE_MODE = ParameterCondition("mode", (2.0,)) _DYNAMIC_PIPE_PARAMETER_GROUPS = ( @@ -1014,6 +1039,7 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): upstream_pressure = max(self.port_1.p, props.p, 1.0) density = max(self.medium.density(upstream_pressure, props.T), 1.0e-12) reynolds = self.reynolds_number(flow, props.T) + friction = self.friction_factor(reynolds) return { "m": self.state.m, "U": self.state.U, @@ -1023,13 +1049,17 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): "u": props.u, "h": props.h, "re": reynolds, - "cm": ( - abs(flow) - * sqrt(props.T) - / max(self.area * upstream_pressure, 1.0e-18) + "cm": _bare_pipe_mass_flow_parameter( + flow, + area=self.area, + diameter=self.diam, + upstream_pressure=upstream_pressure, + upstream_temperature=props.T, + resistance_length=self.le, + friction_factor=friction, ), "v": flow / (density * self.area), - "ff": _reported_friction_factor(self.friction_factor(reynolds)), + "ff": _reported_friction_factor(friction), } def pressure_flow_equation_values(self) -> tuple[float, ...]: @@ -1344,16 +1374,21 @@ class AmesimPnl0002(AmesimPnl0001): 1.0e-12, ) reynolds = self.reynolds_number(flow, upstream_temperature) + friction = self.friction_factor(reynolds) resistance_diagnostics.append( ( reynolds, - ( - abs(flow) - * sqrt(upstream_temperature) - / max(self.area * upstream_pressure, 1.0e-18) + _bare_pipe_mass_flow_parameter( + flow, + area=self.area, + diameter=self.diam, + upstream_pressure=upstream_pressure, + upstream_temperature=upstream_temperature, + resistance_length=self.resistance_length, + friction_factor=friction, ), abs(flow) / (density * self.area), - self.friction_factor(reynolds), + friction, ) ) reynolds, cm, velocity, friction = ( @@ -1699,6 +1734,7 @@ class AmesimPnl0003(DynamicComponent): center_flow = self.resistance_mass_flow() upstream = port_1 if center_flow >= 0.0 else port_2 reynolds = self.reynolds_number(center_flow, upstream.T) + friction = self.friction_factor(reynolds) return { "m1": self.state_1.m, "U1": self.state_1.U, @@ -1716,13 +1752,17 @@ class AmesimPnl0003(DynamicComponent): "h2": port_2.h, "dmctr": center_flow, "re": reynolds, - "cm": ( - abs(center_flow) - * sqrt(upstream.T) - / max(self.area * max(port_1.p, port_2.p, 1.0), 1.0e-18) + "cm": _bare_pipe_mass_flow_parameter( + center_flow, + area=self.area, + diameter=self.diam, + upstream_pressure=max(port_1.p, port_2.p, 1.0), + upstream_temperature=upstream.T, + resistance_length=self.le, + friction_factor=friction, ), "v": center_flow / (max(upstream.rho, 1.0e-12) * self.area), - "ff": _reported_friction_factor(self.friction_factor(reynolds)), + "ff": _reported_friction_factor(friction), } def pressure_flow_equation_values(self) -> tuple[float, ...]: diff --git a/tests/test_amesim_pnl0001_component.py b/tests/test_amesim_pnl0001_component.py index e6bdd43..8a7c514 100644 --- a/tests/test_amesim_pnl0001_component.py +++ b/tests/test_amesim_pnl0001_component.py @@ -1,6 +1,7 @@ from __future__ import annotations import unittest +from math import sqrt from app.simulation.components.amesim.flow.pipes import AmesimPnl0001 from app.simulation.core.medium import IdealGasMedium @@ -148,6 +149,33 @@ class AmesimPnl0001ComponentTests(unittest.TestCase): {"m", "U", "p", "T", "rho", "u", "h", "re", "cm", "v", "ff"}, ) + def test_result_cm_excludes_pipe_flow_coefficient(self) -> None: + pipe = AmesimPnl0001( + "pnl_1", + self.medium, + diam=0.014, + le=0.3, + rr=0.045 / 14.0, + p0=2.5e6, + T0=290.0, + ) + pipe.port_1.p = 1.6e6 + props = pipe.properties() + flow = pipe.mass_flow(pipe.port_1.p, props.p, props.T) + reynolds = pipe.reynolds_number(flow, props.T) + friction = pipe.friction_factor(reynolds) + raw_cq_cm = ( + abs(flow) + * sqrt(props.T) + / (pipe.area * max(pipe.port_1.p, props.p)) + ) + flow_coefficient = sqrt(pipe.diam / (pipe.le * friction)) + + actual = pipe.component_result_values()["cm"] + + self.assertAlmostEqual(actual, raw_cq_cm / flow_coefficient) + self.assertNotAlmostEqual(actual, raw_cq_cm) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_amesim_pnl0002_pnl0003_component.py b/tests/test_amesim_pnl0002_pnl0003_component.py index 0aefc27..e470f1b 100644 --- a/tests/test_amesim_pnl0002_pnl0003_component.py +++ b/tests/test_amesim_pnl0002_pnl0003_component.py @@ -1,7 +1,7 @@ from __future__ import annotations import unittest -from math import log10 +from math import log10, sqrt from app.simulation.components.amesim.flow.pipes import AmesimPnl0002, AmesimPnl0003 from app.simulation.core.medium import IdealGasMedium @@ -200,6 +200,79 @@ class AmesimPnl0002ComponentTests(unittest.TestCase): self.assertAlmostEqual(derivative[0], 0.3) self.assertAlmostEqual(derivative[1], 0.3 * props.h) + def test_result_cm_averages_bare_parameter_for_each_resistance(self) -> None: + pipe = AmesimPnl0002( + "pnl_2", + self.medium, + diam=0.02, + le=2.0, + rr=0.045 / 20.0, + p0=100000.0, + T0=350.0, + ) + center = pipe.properties() + pipe.port_1.p = 130000.0 + pipe.port_2.p = 80000.0 + pipe.update_flow_temperature_references( + { + "port_1": self.medium.specific_enthalpy_at_pressure( + pipe.port_1.p, + 250.0, + ), + "port_2": self.medium.specific_enthalpy_at_pressure( + pipe.port_2.p, + 450.0, + ), + } + ) + + raw_parameters: list[float] = [] + bare_parameters: list[float] = [] + frictions: list[float] = [] + for port_name, port in ( + ("port_1", pipe.port_1), + ("port_2", pipe.port_2), + ): + flow = pipe.port_mass_flow( + port.p, + center.p, + center.T, + port_name=port_name, + ) + if flow >= 0.0: + upstream_pressure = port.p + upstream_temperature = pipe.medium.temperature_from_pressure_enthalpy( + port.p, + pipe._connected_h[port_name], + ) + else: + upstream_pressure = center.p + upstream_temperature = center.T + reynolds = pipe.reynolds_number(flow, upstream_temperature) + friction = pipe.friction_factor(reynolds) + raw = ( + abs(flow) + * sqrt(upstream_temperature) + / (pipe.area * upstream_pressure) + ) + flow_coefficient = sqrt( + pipe.diam / (pipe.resistance_length * friction) + ) + raw_parameters.append(raw) + bare_parameters.append(raw / flow_coefficient) + frictions.append(friction) + + expected = sum(bare_parameters) / 2.0 + collapsed = (sum(raw_parameters) / 2.0) / sqrt( + pipe.diam + / (pipe.resistance_length * (sum(frictions) / 2.0)) + ) + + actual = pipe.component_result_values()["cm"] + + self.assertAlmostEqual(actual, expected) + self.assertGreater(abs(actual - collapsed), 1.0e-8) + class AmesimPnl0003ComponentTests(unittest.TestCase): def setUp(self) -> None: @@ -389,6 +462,36 @@ class AmesimPnl0003ComponentTests(unittest.TestCase): self.assertIn("dmctr", values) self.assertIn("re", values) + def test_result_cm_excludes_center_resistance_flow_coefficient(self) -> None: + pipe = AmesimPnl0003( + "pnl_3", + self.medium, + diam=0.02, + le=0.3, + rr=0.045 / 20.0, + p1_0=10.7e6, + T1_0=263.4, + p2_0=10.68e6, + T2_0=263.42, + ) + port_1 = pipe.properties_1() + port_2 = pipe.properties_2() + flow = pipe.resistance_mass_flow() + upstream = port_1 if flow >= 0.0 else port_2 + reynolds = pipe.reynolds_number(flow, upstream.T) + friction = pipe.friction_factor(reynolds) + raw_cq_cm = ( + abs(flow) + * sqrt(upstream.T) + / (pipe.area * max(port_1.p, port_2.p)) + ) + flow_coefficient = sqrt(pipe.diam / (pipe.le * friction)) + + actual = pipe.component_result_values()["cm"] + + self.assertAlmostEqual(actual, raw_cq_cm / flow_coefficient) + self.assertNotAlmostEqual(actual, raw_cq_cm) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_amesim_pnvo001_fixed_component.py b/tests/test_amesim_pnvo001_fixed_component.py index 9e78ea8..69ceee5 100644 --- a/tests/test_amesim_pnvo001_fixed_component.py +++ b/tests/test_amesim_pnvo001_fixed_component.py @@ -54,6 +54,25 @@ class AmesimPnvo001FixedOpeningComponentTests(unittest.TestCase): self.assertEqual(valve.mass_flow(500000.0, 100000.0), 0.0) + def test_closed_opening_only_zeroes_reported_gas_velocity(self) -> None: + valve = AmesimPnvo001FixedOpening("valve_1", self.medium, opening=0.0) + valve.port_2.p = 500000.0 + valve.port_3.p = 100000.0 + + closed = valve.component_result_values() + self.assertGreater(closed["cm"], 0.0) + self.assertEqual(closed["gasvel"], 0.0) + + valve.opening = 7.0e-15 + roundoff_closed = valve.component_result_values() + self.assertEqual(roundoff_closed["gasvel"], 0.0) + self.assertEqual(roundoff_closed["cm"], closed["cm"]) + + valve.opening = 1.0e-9 + genuinely_open = valve.component_result_values() + self.assertNotEqual(genuinely_open["gasvel"], 0.0) + self.assertEqual(genuinely_open["cm"], closed["cm"]) + def test_mass_flow_is_bidirectional_by_pressure_order(self) -> None: valve = AmesimPnvo001FixedOpening("valve_1", self.medium, opening=0.5) diff --git a/tests/test_amesim_pnvo001_signal_xml.py b/tests/test_amesim_pnvo001_signal_xml.py index 03d09ed..eb28648 100644 --- a/tests/test_amesim_pnvo001_signal_xml.py +++ b/tests/test_amesim_pnvo001_signal_xml.py @@ -244,6 +244,8 @@ class AmesimPnvo001SignalXmlTests(unittest.TestCase): self.assertEqual(result["series"]["step_1.out.signal"][at_event], 1.0) self.assertEqual(result["series"]["valve_1.xv"][before_event], 0.0) self.assertEqual(result["series"]["valve_1.xv"][at_event], 1.0) + self.assertEqual(result["series"]["valve_1.gasvel"][before_event], 0.0) + self.assertNotEqual(result["series"]["valve_1.gasvel"][at_event], 0.0) self.assertGreater( result["series"]["valve_1.port_2.m_flow"][at_event], 0.45,