From a071896834216313325b5f1975488300a64b79a2 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Fri, 17 Jul 2026 09:38:32 +0000 Subject: [PATCH] =?UTF-8?q?=E5=AE=9E=E7=8E=B0test=5Fmql=20PNL0002=E7=AE=A1?= =?UTF-8?q?=E8=B7=AF=E5=8A=A8=E6=80=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim_pneumatic_line.py | 193 ++++++++++++++++++ PythonModels/systems/test_mql.py | 12 ++ .../systems/test_mql_line_parameters.py | 108 ++++++++++ .../systems/test_mql_pneumatic_lines.py | 39 ++++ .../test_amesim_pneumatic_line_components.py | 48 +++++ tests/test_test_mql_pnl0002.py | 92 +++++++++ 6 files changed, 492 insertions(+) create mode 100644 tests/test_test_mql_pnl0002.py diff --git a/PythonModels/components/amesim_pneumatic_line.py b/PythonModels/components/amesim_pneumatic_line.py index 83dcb5a..7ce54f8 100644 --- a/PythonModels/components/amesim_pneumatic_line.py +++ b/PythonModels/components/amesim_pneumatic_line.py @@ -477,6 +477,199 @@ class AmesimPnl0003Pipe(_DarcyPipeResistanceMixin, DynamicComponent): ) +class AmesimPnl0002Pipe(_DarcyPipeResistanceMixin, DynamicComponent): + """First-pass AMESim ``PNL0002`` (R-C-R) pipe. + + The center compliance owns the gas state. Positive connection mass flows + enter that center storage from each external port. + """ + + state_size = 2 + + def __init__( + self, + name: str, + *, + diameter_mm: float, + length_m: float, + relative_roughness: float, + polytropic_constant: float = 1.35, + heat_transfer_coefficient: float = 0.0, + external_temperature_k: float = 293.15, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + pctr_0: float = 101_325.0, + Tctr_0: float = 293.15, + ) -> None: + if diameter_mm <= 0.0: + raise ValueError("diameter_mm must be positive") + if length_m <= 0.0: + raise ValueError("length_m must be positive") + if relative_roughness < 0.0: + raise ValueError("relative_roughness must be non-negative") + if polytropic_constant <= 0.0: + raise ValueError("polytropic_constant must be positive") + if heat_transfer_coefficient < 0.0: + raise ValueError("heat_transfer_coefficient must be non-negative") + if external_temperature_k <= 0.0: + raise ValueError("external_temperature_k must be positive") + + super().__init__(name=name) + self.diameter = diameter_mm * 1.0e-3 + self.length = length_m + self.relative_roughness = relative_roughness + self.polytropic_constant = polytropic_constant + self.heat_transfer_coefficient = heat_transfer_coefficient + self.external_temperature = external_temperature_k + self.gas = gas + self.area = diameter_mm_to_area_m2(diameter_mm) + self.volume = self.area * self.length + self.heat_transfer_area = pi * self.diameter * self.length + self._resistance_length = self.length / 2.0 + + rho0 = gas.density(pctr_0, Tctr_0) + mass0 = rho0 * self.volume + self.state = VolumeState( + m=mass0, + U=mass0 * gas.specific_internal_energy(Tctr_0), + ) + self.port_1 = PortState() + self.port_2 = PortState() + + def get_state_vector(self) -> list[float]: + return self.state.as_vector() + + def set_state_vector(self, values: list[float]) -> None: + self.state = VolumeState.from_vector(values) + + def properties(self) -> ThermodynamicProperties: + if self.state.m <= 0.0: + raise ValueError("pipe mass must stay positive") + temperature = self.gas.temperature_from_internal_energy( + self.state.U / self.state.m + ) + density = self.state.m / self.volume + pressure = self.gas.pressure(density, temperature) + properties = ThermodynamicProperties( + p=pressure, + T=temperature, + rho=density, + u=self.state.U / self.state.m, + h=self.gas.specific_enthalpy(temperature), + ) + self.port_1.p = pressure + self.port_1.h_outflow = properties.h + self.port_2.p = pressure + self.port_2.h_outflow = properties.h + return properties + + def gas_mass_g(self) -> float: + return self.state.m * 1.0e3 + + def port_mass_flow( + self, + *, + port_pressure_pa: float, + port_temperature_k: float, + ) -> float: + """Return mass flow from an external port into the center storage.""" + if port_pressure_pa <= 0.0: + raise ValueError("port_pressure_pa must be positive") + if port_temperature_k <= 0.0: + raise ValueError("port_temperature_k must be positive") + + center = self.properties() + pressure_difference = port_pressure_pa - center.p + if pressure_difference == 0.0: + return 0.0 + upstream_pressure = max(port_pressure_pa, center.p) + upstream_temperature = ( + port_temperature_k if pressure_difference > 0.0 else center.T + ) + density = self.gas.density(upstream_pressure, upstream_temperature) + magnitude = self._mass_flow_for_resistance_pressure_drop( + abs(pressure_difference), + density=density, + temperature=upstream_temperature, + ) + return magnitude if pressure_difference > 0.0 else -magnitude + + def _mass_flow_for_resistance_pressure_drop( + self, + pressure_drop_pa: float, + *, + density: float, + temperature: float, + ) -> float: + original_length = self.length + self.length = self._resistance_length + try: + return self._mass_flow_for_pressure_drop( + pressure_drop_pa, + density=density, + temperature=temperature, + ) + finally: + self.length = original_length + + def diagnostics( + self, + *, + mass_flow_kg_s: float, + temperature_k: float | None = None, + ) -> AmesimPnl0001Diagnostics: + properties = self.properties() + temperature = temperature_k or properties.T + reynolds = self._reynolds_number(mass_flow_kg_s, temperature) + friction_factor = self._friction_factor(reynolds) + velocity = mass_flow_kg_s / (properties.rho * self.area) + original_length = self.length + self.length = self._resistance_length + try: + pressure_drop = self._darcy_pressure_drop( + mass_flow_kg_s, + density=properties.rho, + temperature=temperature, + ) + finally: + self.length = original_length + return AmesimPnl0001Diagnostics( + mass_flow_kg_s=mass_flow_kg_s, + reynolds_number=reynolds, + gas_velocity_m_s=velocity, + friction_factor=friction_factor, + pressure_drop_pa=pressure_drop, + ) + + def derivatives_from_connections( + self, + *, + port_1_m_flow: float, + connected_h_1: float, + port_2_m_flow: float, + connected_h_2: float, + ) -> VolumeState: + center = self.properties() + inlet_h_1 = self.connection_inlet_enthalpy( + port_m_flow=port_1_m_flow, + connected_h=connected_h_1, + internal_h=center.h, + ) + inlet_h_2 = self.connection_inlet_enthalpy( + port_m_flow=port_2_m_flow, + connected_h=connected_h_2, + internal_h=center.h, + ) + heat_flow = ( + self.heat_transfer_coefficient + * self.heat_transfer_area + * (self.external_temperature - center.T) + ) + 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, + ) + + class AmesimPnl00rPipe(_DarcyPipeResistanceMixin, AlgebraicComponent): """First-pass AMESim ``PNL00R`` (R) pipe resistance.""" diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 35d28de..ffff67a 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3860,6 +3860,7 @@ class TestMqlSystem: self.archive_path = archive_path or Path(__file__).resolve().parents[2] / AMESIM_ARCHIVE_RELATIVE_PATH self.network = SimulationNetwork(name=MODEL_NAME) self.pnl0001_assembly = self._build_pnl0001_assembly() + self.pnl0002_assembly = self._build_pnl0002_assembly() self.pnl0003_assembly = self._build_pnl0003_assembly() self.pnl00r_assembly = self._build_pnl00r_assembly() self.node3_assembly = self._build_node3_assembly() @@ -3893,6 +3894,13 @@ class TestMqlSystem: return build_test_mql_pnl0001_assembly(self.archive_path) + def _build_pnl0002_assembly(self): + from PythonModels.systems.test_mql_pneumatic_lines import ( + build_test_mql_pnl0002_assembly, + ) + + return build_test_mql_pnl0002_assembly(self.archive_path) + def _build_pnl0003_assembly(self): from PythonModels.systems.test_mql_pneumatic_lines import ( build_test_mql_pnl0003_assembly, @@ -3935,6 +3943,10 @@ class TestMqlSystem: def typed_pnl0001_line_count(self) -> int: return len(self.pnl0001_assembly.lines) + @property + def typed_pnl0002_line_count(self) -> int: + return len(self.pnl0002_assembly.lines) + @property def typed_pnl0003_line_count(self) -> int: return len(self.pnl0003_assembly.lines) diff --git a/PythonModels/systems/test_mql_line_parameters.py b/PythonModels/systems/test_mql_line_parameters.py index ff4c87d..e79298f 100644 --- a/PythonModels/systems/test_mql_line_parameters.py +++ b/PythonModels/systems/test_mql_line_parameters.py @@ -35,6 +35,29 @@ class TestMqlPnl0001Spec: return self.initial_gauge_pressure_pa + AMESIM_REFERENCE_PRESSURE_PA +@dataclass(frozen=True) +class TestMqlPnl0002Spec: + alias: str + source_component: str + source_port: str + target_component: str + target_port: str + diameter_mm: float + length_m: float + relative_roughness: float + polytropic_constant: float + heat_transfer_coefficient: float + external_temperature_k: float + gas_type_index: int + mode: int + initial_center_temperature_k: float + initial_center_gauge_pressure_pa: float + + @property + def initial_center_absolute_pressure_pa(self) -> float: + return self.initial_center_gauge_pressure_pa + AMESIM_REFERENCE_PRESSURE_PA + + @dataclass(frozen=True) class TestMqlPnl0003Spec: alias: str @@ -153,6 +176,82 @@ def load_test_mql_pnl0001_specs( return tuple(specs) +def load_test_mql_pnl0002_specs( + archive_path: str | Path, + *, + cir_member: str = "test_mql_.cir", +) -> tuple[TestMqlPnl0002Spec, ...]: + """Load resolved PNL0002 geometry and center compliance initial state.""" + with tarfile.open(archive_path) as archive: + cir_file = archive.extractfile(cir_member) + if cir_file is None: + raise ValueError(f"Missing AMESim circuit member: {cir_member}") + cir_text = cir_file.read().decode("latin1") + + numeric_globals = { + name: value + for name, expression in GLOBAL_PARAMETERS.items() + if (value := resolve_numeric_expression(expression, {})) is not None + } + connections = { + str(connection["alias"]): connection + for connection in CONNECTION_SPECS + if connection["submodel"] == "PNL0002" + } + specs = [] + for block in re.findall(r".*?", cir_text, flags=re.DOTALL): + if _optional_text(block, "SUB_NAME") != "PNL0002": + continue + alias = _required_text(block, "ALIAS") + connection = connections.get(alias) + if connection is None: + raise ValueError(f"PNL0002 line {alias!r} is absent from CONNECTION_SPECS") + real_parameters = _parameter_expressions(block, "RPARAM") + integer_parameters = _parameter_expressions(block, "IPARAM") + state_values = _ivar_values(block) + specs.append( + TestMqlPnl0002Spec( + alias=alias, + source_component=str(connection["source_component"]), + source_port=str(connection["source_port"]), + target_component=str(connection["target_component"]), + target_port=str(connection["target_port"]), + diameter_mm=_required_numeric( + alias, "diam", real_parameters, numeric_globals + ), + length_m=_required_numeric(alias, "le", real_parameters, numeric_globals), + relative_roughness=_required_numeric( + alias, "rr", real_parameters, numeric_globals + ), + polytropic_constant=_required_numeric( + alias, "k", real_parameters, numeric_globals + ), + heat_transfer_coefficient=_required_numeric( + alias, "kth", real_parameters, numeric_globals + ), + external_temperature_k=_required_numeric( + alias, "extemp", real_parameters, numeric_globals + ), + gas_type_index=int( + _required_numeric(alias, "gi", integer_parameters, numeric_globals) + ), + mode=int( + _required_numeric(alias, "mode", integer_parameters, numeric_globals) + ), + initial_center_temperature_k=_required_numeric( + alias, "tctr", state_values, numeric_globals + ), + initial_center_gauge_pressure_pa=_required_numeric( + alias, "pctr", state_values, numeric_globals + ), + ) + ) + if set(connections) != {spec.alias for spec in specs}: + missing = sorted(set(connections) - {spec.alias for spec in specs}) + raise ValueError(f"Missing PNL0002 parameter blocks: {missing}") + return tuple(specs) + + def load_test_mql_pnl0003_specs( archive_path: str | Path, *, @@ -306,6 +405,15 @@ def _parameter_expressions(block: str, tag_name: str) -> dict[str, str]: return parameters +def _ivar_values(block: str) -> dict[str, str]: + values = {} + for variable_block in re.findall(r".*?", block, flags=re.DOTALL): + value = _optional_text(variable_block, "VALUE") + if value: + values[_required_text(variable_block, "VARNAME")] = value + return values + + def _evar_values(block: str) -> dict[str, str]: values = {} for variable_block in re.findall(r".*?", block, flags=re.DOTALL): diff --git a/PythonModels/systems/test_mql_pneumatic_lines.py b/PythonModels/systems/test_mql_pneumatic_lines.py index 50ebc54..18b42ed 100644 --- a/PythonModels/systems/test_mql_pneumatic_lines.py +++ b/PythonModels/systems/test_mql_pneumatic_lines.py @@ -9,14 +9,17 @@ from PythonModels.components.amesim_pneumatic import ( ) from PythonModels.components.amesim_pneumatic_line import ( AmesimPnl0001Pipe, + AmesimPnl0002Pipe, AmesimPnl0003Pipe, AmesimPnl00rPipe, ) from PythonModels.systems.test_mql_line_parameters import ( TestMqlPnl0001Spec, + TestMqlPnl0002Spec, TestMqlPnl0003Spec, TestMqlPnl00rSpec, load_test_mql_pnl0001_specs, + load_test_mql_pnl0002_specs, load_test_mql_pnl0003_specs, load_test_mql_pnl00r_specs, ) @@ -34,6 +37,18 @@ class TestMqlPnl0001Assembly: raise KeyError(alias) +@dataclass(frozen=True) +class TestMqlPnl0002Assembly: + specs: tuple[TestMqlPnl0002Spec, ...] + lines: dict[str, AmesimPnl0002Pipe] + + def spec(self, alias: str) -> TestMqlPnl0002Spec: + for spec in self.specs: + if spec.alias == alias: + return spec + raise KeyError(alias) + + @dataclass(frozen=True) class TestMqlPnl0003Assembly: specs: tuple[TestMqlPnl0003Spec, ...] @@ -82,6 +97,30 @@ def build_test_mql_pnl0001_assembly( return TestMqlPnl0001Assembly(specs=specs, lines=lines) +def build_test_mql_pnl0002_assembly( + archive_path: str | Path, + *, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, +) -> TestMqlPnl0002Assembly: + specs = load_test_mql_pnl0002_specs(archive_path) + lines = { + spec.alias: AmesimPnl0002Pipe( + name=spec.alias, + diameter_mm=spec.diameter_mm, + length_m=spec.length_m, + relative_roughness=spec.relative_roughness, + polytropic_constant=spec.polytropic_constant, + heat_transfer_coefficient=spec.heat_transfer_coefficient, + external_temperature_k=spec.external_temperature_k, + gas=gas, + pctr_0=spec.initial_center_absolute_pressure_pa, + Tctr_0=spec.initial_center_temperature_k, + ) + for spec in specs + } + return TestMqlPnl0002Assembly(specs=specs, lines=lines) + + def build_test_mql_pnl0003_assembly( archive_path: str | Path, *, diff --git a/tests/test_amesim_pneumatic_line_components.py b/tests/test_amesim_pneumatic_line_components.py index 66a101b..4f73603 100644 --- a/tests/test_amesim_pneumatic_line_components.py +++ b/tests/test_amesim_pneumatic_line_components.py @@ -3,11 +3,59 @@ from __future__ import annotations import unittest from PythonModels.components.amesim_pneumatic_line import ( + AmesimPnl0002Pipe, AmesimPnl0003Pipe, AmesimPnl00rPipe, ) +class AmesimPnl0002PipeTests(unittest.TestCase): + def setUp(self) -> None: + self.pipe = AmesimPnl0002Pipe( + name="pneumatic_86", + diameter_mm=20.0, + length_m=2.0, + relative_roughness=0.045 / 20.0, + pctr_0=100_000.0, + Tctr_0=293.15, + ) + + def test_initial_state_uses_center_compliance(self) -> None: + properties = self.pipe.properties() + + self.assertAlmostEqual(self.pipe.volume, 6.283185307179586e-4) + self.assertEqual(len(self.pipe.get_state_vector()), 2) + self.assertAlmostEqual(properties.p, 100_000.0, delta=1.0e-6) + self.assertAlmostEqual(properties.T, 293.15) + self.assertAlmostEqual(self.pipe.port_1.p, properties.p) + self.assertAlmostEqual(self.pipe.port_2.p, properties.p) + + def test_port_flow_enters_center_from_higher_external_pressure(self) -> None: + forward = self.pipe.port_mass_flow( + port_pressure_pa=101_000.0, + port_temperature_k=293.15, + ) + reverse = self.pipe.port_mass_flow( + port_pressure_pa=99_000.0, + port_temperature_k=293.15, + ) + + self.assertGreater(forward, 0.0) + self.assertLess(reverse, 0.0) + + def test_connection_derivatives_conserve_two_external_port_flows(self) -> None: + properties = self.pipe.properties() + + derivative = self.pipe.derivatives_from_connections( + port_1_m_flow=0.2, + connected_h_1=properties.h + 1000.0, + port_2_m_flow=-0.1, + connected_h_2=properties.h - 1000.0, + ) + + self.assertAlmostEqual(derivative.m, 0.1) + + class AmesimPnl0003PipeTests(unittest.TestCase): def setUp(self) -> None: self.pipe = AmesimPnl0003Pipe( diff --git a/tests/test_test_mql_pnl0002.py b/tests/test_test_mql_pnl0002.py new file mode 100644 index 0000000..4139cb3 --- /dev/null +++ b/tests/test_test_mql_pnl0002.py @@ -0,0 +1,92 @@ +from __future__ import annotations + +import unittest +from pathlib import Path + +from PythonModels.reporting.amesim_results import load_test_mql_amesim_results +from PythonModels.systems.test_mql import TestMqlSystem +from PythonModels.systems.test_mql_line_parameters import load_test_mql_pnl0002_specs +from PythonModels.systems.test_mql_pneumatic_lines import build_test_mql_pnl0002_assembly + + +REPO_ROOT = Path(__file__).resolve().parents[1] +TEST_MQL_AME = REPO_ROOT / "AmesimModels" / "test_mql.ame" + + +class TestMqlPnl0002Tests(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + cls.pnl0002_specs = load_test_mql_pnl0002_specs(TEST_MQL_AME) + cls.pnl0002_assembly = build_test_mql_pnl0002_assembly(TEST_MQL_AME) + cls.results = load_test_mql_amesim_results(TEST_MQL_AME) + + def test_loads_all_real_pnl0002_parameters_from_cir(self) -> None: + self.assertEqual(len(self.pnl0002_specs), 8) + spec = self.pnl0002_assembly.spec("pneumatic_86") + + self.assertEqual(spec.source_component, "pnnode4_17") + self.assertEqual(spec.source_port, "port_1") + self.assertEqual(spec.target_component, "pnnode4_18") + self.assertEqual(spec.target_port, "port_3") + self.assertEqual(spec.diameter_mm, 20.0) + self.assertEqual(spec.length_m, 2.0) + self.assertAlmostEqual(spec.relative_roughness, 0.045 / 20.0) + self.assertEqual(spec.gas_type_index, 1) + self.assertEqual(spec.mode, 2) + self.assertAlmostEqual(spec.initial_center_temperature_k, 293.15) + self.assertAlmostEqual(spec.initial_center_gauge_pressure_pa, -1300.0) + self.assertAlmostEqual(spec.initial_center_absolute_pressure_pa, 100_000.0) + + def test_builds_pnl0002_center_compliance_line_components(self) -> None: + self.assertEqual(len(self.pnl0002_assembly.lines), 8) + self.assertTrue( + all(len(line.get_state_vector()) == 2 for line in self.pnl0002_assembly.lines.values()) + ) + + def test_pneumatic_86_initial_observables_match_amesim_baseline(self) -> None: + pipe = self.pnl0002_assembly.lines["pneumatic_86"] + properties = pipe.properties() + + self.assertAlmostEqual( + properties.p - 101_300.0, + self.results.series("pctr@pneumatic_86")[0], + delta=1.0e-6, + ) + self.assertAlmostEqual( + properties.T, + self.results.series("tctr@pneumatic_86")[0], + delta=1.0e-12, + ) + self.assertAlmostEqual( + pipe.gas_mass_g(), + self.results.series("mgas@pneumatic_86")[0], + delta=0.001, + ) + + def test_amesim_mass_balance_uses_opposite_saved_port_flow_signs(self) -> None: + index = 500 + times = self.results.times + mass = self.results.series("mgas@pneumatic_86") + finite_difference_g_s = ( + mass[index + 1] - mass[index - 1] + ) / (times[index + 1] - times[index - 1]) + saved_port_flow_g_s = ( + self.results.series("dm1@pneumatic_86")[index] + + self.results.series("dm2@pneumatic_86")[index] + ) + + self.assertAlmostEqual( + finite_difference_g_s, + -saved_port_flow_g_s, + delta=1.0e-5, + ) + + def test_system_exposes_real_pnl0002_assembly(self) -> None: + system = TestMqlSystem() + + self.assertEqual(system.typed_pnl0002_line_count, 8) + self.assertIn("pneumatic_86", system.pnl0002_assembly.lines) + + +if __name__ == "__main__": + unittest.main()