From 40c72422ff6919118d969d263d75d9b96ede8f43 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Wed, 5 Aug 2026 12:59:05 +0000 Subject: [PATCH] =?UTF-8?q?=E5=AF=B9=E9=BD=90PNVO=E8=BF=91=E7=AD=89?= =?UTF-8?q?=E5=8E=8B=E5=B1=82=E6=B5=81=E5=B9=B3=E6=BB=91?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim/flow/orifices.py | 74 +++++++++++---- tests/test_amesim_pnvo001_fixed_component.py | 90 +++++++++++++++++++ 2 files changed, 149 insertions(+), 15 deletions(-) diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py index c1c0f60..3b5b58e 100644 --- a/app/simulation/components/amesim/flow/orifices.py +++ b/app/simulation/components/amesim/flow/orifices.py @@ -1,7 +1,7 @@ from __future__ import annotations from collections.abc import Mapping -from math import isclose, sqrt +from math import isclose, log, sqrt, tanh from app.simulation.components.amesim.gases import ( AMESIM_GAS_INDEX_PARAMETER, @@ -32,6 +32,8 @@ _FLOW_COEFFICIENT_OPTIONS = ( _FLOWSET_USES_CQ = (ParameterCondition("flowset", (1.0,)),) _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 _PNOR001_FLOW_COEFFICIENT_GROUP = ParameterGroupDisplaySpec( id="flow_coefficient", label="流量系数", @@ -565,6 +567,31 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): 1.0, ) + @staticmethod + def _subsonic_mass_flow_parameter( + *, + pressure_ratio: float, + gamma_s: float, + density: float, + upstream_temperature: float, + upstream_pressure: float, + ) -> float: + expansion = ( + pressure_ratio ** (2.0 * gamma_s) + - pressure_ratio ** (1.0 + gamma_s) + ) + return sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * density + * upstream_temperature + / upstream_pressure + * expansion, + 0.0, + ) + ) + def mass_flow(self, p_2: float, p_3: float) -> float: if p_2 == p_3 or self.effective_area == 0.0: return 0.0 @@ -602,6 +629,7 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): 1.0 / (1.0 - gamma_s) ) if pressure_ratio <= critical_ratio: + effective_pressure_ratio = critical_ratio mass_flow_parameter = ( sqrt(2.0 / (1.0 + gamma_s) * density * T_up / p_up) * (2.0 * gamma_s / (gamma_s + 1.0)) @@ -611,20 +639,13 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): 2.0 / (1.0 + gamma_s) * p_up / density ) else: - expansion = ( - pressure_ratio ** (2.0 * gamma_s) - - pressure_ratio ** (1.0 + gamma_s) - ) - mass_flow_parameter = sqrt( - max( - 2.0 - / (1.0 - gamma_s) - * density - * T_up - / p_up - * expansion, - 0.0, - ) + effective_pressure_ratio = pressure_ratio + mass_flow_parameter = self._subsonic_mass_flow_parameter( + pressure_ratio=pressure_ratio, + gamma_s=gamma_s, + density=density, + upstream_temperature=T_up, + upstream_pressure=p_up, ) gas_velocity = sqrt( max( @@ -636,6 +657,29 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): 0.0, ) ) + + # AMESim's gas_cm_prc_ applies this factor continuously over the + # complete pressure-ratio range. It is effectively one outside the + # near-equal-pressure region and makes Cm (and vena-contracta + # velocity) approach zero quadratically as the pressure ratio tends + # to one. The reference Cm intentionally reuses the current gamma_s. + reference_mass_flow_parameter = self._subsonic_mass_flow_parameter( + pressure_ratio=_PN_PRESSURE_RATIO_ACCURACY, + gamma_s=gamma_s, + density=density, + upstream_temperature=T_up, + upstream_pressure=p_up, + ) + if mass_flow_parameter > 0.0 and reference_mass_flow_parameter > 0.0: + smoothing_argument = ( + _PN_LAMINAR_SMOOTHING_GAIN + * abs(mass_flow_parameter / reference_mass_flow_parameter) + * log(effective_pressure_ratio) + / log(_PN_PRESSURE_RATIO_ACCURACY) + ) + smoothing_factor = tanh(max(smoothing_argument, 0.0)) + mass_flow_parameter *= smoothing_factor + gas_velocity *= smoothing_factor return mass_flow_parameter, gas_velocity def _one_way_mass_flow( diff --git a/tests/test_amesim_pnvo001_fixed_component.py b/tests/test_amesim_pnvo001_fixed_component.py index b69a003..c3a4740 100644 --- a/tests/test_amesim_pnvo001_fixed_component.py +++ b/tests/test_amesim_pnvo001_fixed_component.py @@ -144,6 +144,96 @@ class AmesimPnvo001FixedOpeningComponentTests(unittest.TestCase): delta=1.0e-8, ) + def test_helium_laminar_transition_matches_amesim_baseline(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + valve = AmesimPnvo001FixedOpening( + "valve_1", + medium, + cq=0.45, + area0=78.5e-6, + opening=1.0, + ) + cases = ( + # p_up [PaA], p_down [PaA], T_up [K], dm [g/s], Cm, velocity [m/s] + ( + 10_837_357.913884956, + 10_835_428.362241976, + 255.45507181057303, + 9.770711291393273, + 4.07922685863417e-4, + 13.951480313568323, + ), + ( + 10_735_827.842429036, + 10_735_543.135671863, + 254.49635939397584, + 3.4711877823059916, + 1.4601624857015259e-4, + 4.982778680076626, + ), + ( + 10_703_773.135050302, + 10_703_621.857082237, + 254.1925398502598, + 1.53560321377487, + 6.475023461883112e-5, + 2.208065833961227, + ), + ( + 10_682_359.261700785, + 10_682_301.022976216, + 253.98926908210996, + 0.25626141624058846, + 1.0822847918204546e-5, + 0.36890236539556787, + ), + ( + 10_678_100.351449525, + 10_678_093.704329137, + 253.94881208705195, + 0.0033658306146710985, + 1.421965896377586e-7, + 0.004846389527047451, + ), + ) + + for ( + p_up, + p_down, + T_up, + expected_dm, + expected_cm, + expected_velocity, + ) in cases: + with self.subTest(pressure_difference=p_up - p_down): + mass_flow_parameter, gas_velocity = ( + valve._one_way_flow_characteristics( + upstream_pressure=p_up, + downstream_pressure=p_down, + upstream_temperature=T_up, + ) + ) + mass_flow_g_s = ( + valve._one_way_mass_flow( + upstream_pressure=p_up, + downstream_pressure=p_down, + upstream_temperature=T_up, + ) + * 1.0e3 + ) + + self.assertAlmostEqual(mass_flow_g_s, expected_dm, delta=1.0e-11) + self.assertAlmostEqual( + mass_flow_parameter, + expected_cm, + delta=1.0e-15, + ) + self.assertAlmostEqual( + gas_velocity, + expected_velocity, + delta=5.0e-10, + ) + if __name__ == "__main__": unittest.main()