From cdead7897798b9a2f57772f1b6a3bd4b5b241b08 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Wed, 15 Jul 2026 09:16:11 +0000 Subject: [PATCH] =?UTF-8?q?=E8=A1=A5=E5=85=85=E6=B0=A6=E6=B0=94Peng-Robins?= =?UTF-8?q?on=E7=89=A9=E6=80=A7=E5=BA=93?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- AmesimModels/test_mql/README.md | 26 +++- PythonModels/core/peng_robinson.py | 195 ++++++++++++++++++++++++ PythonModels/systems/test_mql_config.py | 3 + tests/test_peng_robinson.py | 65 ++++++++ tests/test_test_mql_config.py | 2 + 5 files changed, 289 insertions(+), 2 deletions(-) create mode 100644 PythonModels/core/peng_robinson.py create mode 100644 tests/test_peng_robinson.py diff --git a/AmesimModels/test_mql/README.md b/AmesimModels/test_mql/README.md index 7455546..5a263b3 100644 --- a/AmesimModels/test_mql/README.md +++ b/AmesimModels/test_mql/README.md @@ -23,6 +23,7 @@ - `PythonModels/systems/test_mql_config.py` - 提供 `TestMqlConfig.from_amesim_specs()`。 - 负责把 AMESim 全局参数和组件参数解析成后续可用的 Python 配置对象。 + - `test_mql` 默认工质已按 AMESim 模型确认设为氦气。 - 已支持常量、全局参数引用、四则运算、括号和 AMESim 风格的 `^` 指数表达式。 - 非数值文本会保留为 `None`，后续按具体子模型显式处理。 @@ -34,7 +35,16 @@ - 保护组件数、连接数、状态数、全局参数和关键子模型计数。 - `tests/test_test_mql_config.py` - - 保护全局参数解析、AMESim 表达式解析、按子模型分组和典型组件参数解析。 + - 保护全局参数解析、AMESim 表达式解析、按子模型分组、默认氦气工质和典型组件参数解析。 + +- `PythonModels/core/peng_robinson.py` + - 提供 Peng-Robinson 状态方程小物性库。 + - 当前覆盖压缩因子、摩尔体积、密度和由密度反算压力。 + - 已内置 `HELIUM_PR`，供 `test_mql` 默认使用。 + +## 物性约定 + +AMESim 模型中 `test_mql` 使用氦气，Python 侧当前通过 `HELIUM_PR` 使用 Peng-Robinson 状态方程计算气体压缩因子和密度。当前物性层先覆盖状态方程相关量，完整焓/内能偏差函数后续在接气室能量方程时再补。 ## 当前提取结果 @@ -63,10 +73,22 @@ - `UD00`: 2 - `PNGD00`: 1 + +## 结果对齐原则 + +后续 Python 仿真结果必须以 AMESim 为基准做一一对应校验： + +1. 组件别名优先沿用 AMESim 原名，例如 `pn_brp2_8`、`pn_c1_8`、`pn_orifice_18`。 +2. 子模型族优先沿用 AMESim 子模型名，例如 `PNRP17`、`PNCH012`、`PNOR001`、`PNVO001`。 +3. 输出变量命名要能追溯到 AMESim 的 `Data_Path`，避免 Python 侧改名后无法对比。 +4. Python 每完成一批物理方程，都应拿 AMESim 导出的 CSV 结果做误差对比。 +5. 当前 `.ame` 包内的 `test_mql_.results` 是 AMESim 二进制结果文件，暂不作为直接解析基准；建议后续从 AMESim 导出 CSV 结果放到本目录下。 +6. 在没有 AMESim CSV 基准前，不能把 Python 输出描述为“已经和 AMESim 几乎一致”，只能说明完成了结构、参数或某类组件方程。 + ## 验证方式 ```bash -python3 -m py_compile PythonModels/systems/test_mql.py PythonModels/systems/test_mql_config.py PythonModels/scripts/run_test_mql.py +python3 -m py_compile PythonModels/core/peng_robinson.py PythonModels/systems/test_mql.py PythonModels/systems/test_mql_config.py PythonModels/scripts/run_test_mql.py python3 -m PythonModels.scripts.run_test_mql python3 -m unittest discover -s tests -t . ``` diff --git a/PythonModels/core/peng_robinson.py b/PythonModels/core/peng_robinson.py new file mode 100644 index 0000000..488bec0 --- /dev/null +++ b/PythonModels/core/peng_robinson.py @@ -0,0 +1,195 @@ +from __future__ import annotations + +from dataclasses import dataclass +from math import acos, cos, isfinite, pi, sqrt + +UNIVERSAL_GAS_CONSTANT = 8.31446261815324 + + +@dataclass(frozen=True) +class PengRobinsonFluid: + """Pure-fluid Peng-Robinson equation-of-state helper. + + The class intentionally covers the equation-of-state layer first: pressure, + compressibility factor, molar volume, and density. Caloric departure + properties are left out until the test_mql energy equations need them. + """ + + name: str + molar_mass: float + critical_temperature: float + critical_pressure: float + acentric_factor: float + + @property + def specific_gas_constant(self) -> float: + return UNIVERSAL_GAS_CONSTANT / self.molar_mass + + @property + def a_parameter(self) -> float: + return ( + 0.45724 + * UNIVERSAL_GAS_CONSTANT + * UNIVERSAL_GAS_CONSTANT + * self.critical_temperature + * self.critical_temperature + / self.critical_pressure + ) + + @property + def b_parameter(self) -> float: + return 0.07780 * UNIVERSAL_GAS_CONSTANT * self.critical_temperature / self.critical_pressure + + @property + def kappa(self) -> float: + omega = self.acentric_factor + return 0.37464 + 1.54226 * omega - 0.26992 * omega * omega + + def alpha(self, temperature: float) -> float: + self._validate_temperature(temperature) + reduced_temperature = temperature / self.critical_temperature + return (1.0 + self.kappa * (1.0 - sqrt(reduced_temperature))) ** 2.0 + + def attractive_parameter(self, temperature: float) -> float: + return self.a_parameter * self.alpha(temperature) + + def pressure_from_molar_volume(self, temperature: float, molar_volume: float) -> float: + self._validate_temperature(temperature) + if molar_volume <= self.b_parameter: + raise ValueError("Molar volume must be larger than Peng-Robinson b parameter.") + a_alpha = self.attractive_parameter(temperature) + b = self.b_parameter + repulsive = UNIVERSAL_GAS_CONSTANT * temperature / (molar_volume - b) + attractive = a_alpha / (molar_volume * (molar_volume + b) + b * (molar_volume - b)) + return repulsive - attractive + + def pressure_from_density(self, temperature: float, density: float) -> float: + if density <= 0.0: + raise ValueError("Density must be positive.") + return self.pressure_from_molar_volume(temperature, self.molar_mass / density) + + def reduced_parameters(self, pressure: float, temperature: float) -> tuple[float, float]: + self._validate_pressure_temperature(pressure, temperature) + a_alpha = self.attractive_parameter(temperature) + b = self.b_parameter + A = a_alpha * pressure / (UNIVERSAL_GAS_CONSTANT * UNIVERSAL_GAS_CONSTANT * temperature * temperature) + B = b * pressure / (UNIVERSAL_GAS_CONSTANT * temperature) + return A, B + + def compressibility_roots(self, pressure: float, temperature: float) -> tuple[float, ...]: + A, B = self.reduced_parameters(pressure, temperature) + coefficients = ( + -(1.0 - B), + A - 3.0 * B * B - 2.0 * B, + -(A * B - B * B - B * B * B), + ) + roots = _real_cubic_roots(*coefficients) + physical_roots = tuple(sorted(root for root in roots if root > B and isfinite(root))) + if not physical_roots: + raise ValueError("Peng-Robinson cubic produced no physical compressibility root.") + return physical_roots + + def compressibility_factor( + self, + pressure: float, + temperature: float, + phase: str = "vapor", + ) -> float: + roots = self.compressibility_roots(pressure, temperature) + if phase == "vapor": + return roots[-1] + if phase == "liquid": + return roots[0] + if phase == "stable-single-root": + return roots[-1] + raise ValueError(f"Unsupported phase selector: {phase!r}") + + def molar_volume( + self, + pressure: float, + temperature: float, + phase: str = "vapor", + ) -> float: + z = self.compressibility_factor(pressure, temperature, phase=phase) + return z * UNIVERSAL_GAS_CONSTANT * temperature / pressure + + def density( + self, + pressure: float, + temperature: float, + phase: str = "vapor", + ) -> float: + return self.molar_mass / self.molar_volume(pressure, temperature, phase=phase) + + @staticmethod + def _validate_temperature(temperature: float) -> None: + if temperature <= 0.0: + raise ValueError("Temperature must be positive.") + + @classmethod + def _validate_pressure_temperature(cls, pressure: float, temperature: float) -> None: + if pressure <= 0.0: + raise ValueError("Pressure must be positive.") + cls._validate_temperature(temperature) + +HELIUM_PR = PengRobinsonFluid( + name="helium", + molar_mass=0.004002602, + critical_temperature=5.1953, + critical_pressure=227_460.0, + acentric_factor=-0.385, +) + +NITROGEN_PR = PengRobinsonFluid( + name="nitrogen", + molar_mass=0.0280134, + critical_temperature=126.192, + critical_pressure=3.3958e6, + acentric_factor=0.0372, +) + +AIR_PR = PengRobinsonFluid( + name="air", + molar_mass=0.02896513, + critical_temperature=132.5306, + critical_pressure=3.786e6, + acentric_factor=0.0335, +) + + +def _real_cubic_roots(a: float, b: float, c: float) -> tuple[float, ...]: + """Return real roots for x**3 + a*x**2 + b*x + c = 0.""" + + depressed_p = b - a * a / 3.0 + depressed_q = 2.0 * a * a * a / 27.0 - a * b / 3.0 + c + discriminant = (depressed_q / 2.0) ** 2.0 + (depressed_p / 3.0) ** 3.0 + offset = -a / 3.0 + tolerance = 1e-14 + + if discriminant > tolerance: + sqrt_discriminant = sqrt(discriminant) + u = _real_cube_root(-depressed_q / 2.0 + sqrt_discriminant) + v = _real_cube_root(-depressed_q / 2.0 - sqrt_discriminant) + return (u + v + offset,) + + if abs(discriminant) <= tolerance: + u = _real_cube_root(-depressed_q / 2.0) + return tuple(sorted({2.0 * u + offset, -u + offset})) + + if depressed_p >= 0.0: + raise ValueError("Unexpected cubic state with three real roots and non-negative p.") + radius = 2.0 * sqrt(-depressed_p / 3.0) + argument = (3.0 * depressed_q / (2.0 * depressed_p)) * sqrt(-3.0 / depressed_p) + argument = max(-1.0, min(1.0, argument)) + theta = acos(argument) / 3.0 + roots = [ + radius * cos(theta - 2.0 * pi * index / 3.0) + offset + for index in range(3) + ] + return tuple(sorted(roots)) + + +def _real_cube_root(value: float) -> float: + if value == 0.0: + return 0.0 + return (1.0 if value > 0.0 else -1.0) * abs(value) ** (1.0 / 3.0) diff --git a/PythonModels/systems/test_mql_config.py b/PythonModels/systems/test_mql_config.py index b227752..487936c 100644 --- a/PythonModels/systems/test_mql_config.py +++ b/PythonModels/systems/test_mql_config.py @@ -6,6 +6,7 @@ from dataclasses import dataclass from math import isfinite from typing import Any +from PythonModels.core.peng_robinson import HELIUM_PR, PengRobinsonFluid from PythonModels.systems.test_mql import COMPONENT_SPECS, GLOBAL_PARAMETERS @@ -58,6 +59,7 @@ class TestMqlResolvedComponent: class TestMqlConfig: raw_global_parameters: dict[str, str] global_parameters: dict[str, float] + fluid: PengRobinsonFluid components: tuple[TestMqlResolvedComponent, ...] @classmethod @@ -75,6 +77,7 @@ class TestMqlConfig: return cls( raw_global_parameters=raw_globals, global_parameters=numeric_globals, + fluid=HELIUM_PR, components=components, ) diff --git a/tests/test_peng_robinson.py b/tests/test_peng_robinson.py new file mode 100644 index 0000000..cc8590d --- /dev/null +++ b/tests/test_peng_robinson.py @@ -0,0 +1,65 @@ +from __future__ import annotations + +import unittest + +from PythonModels.core.peng_robinson import AIR_PR, HELIUM_PR, NITROGEN_PR, PengRobinsonFluid + + +class PengRobinsonTest(unittest.TestCase): + def test_helium_near_ideal_at_atmospheric_condition(self) -> None: + z = HELIUM_PR.compressibility_factor(101_325.0, 300.0) + density = HELIUM_PR.density(101_325.0, 300.0) + + self.assertAlmostEqual(z, 1.0, delta=1.0e-3) + self.assertAlmostEqual(density, 0.1625, delta=0.002) + + def test_helium_density_pressure_round_trip(self) -> None: + pressure = 15.3e6 + temperature = 293.15 + density = HELIUM_PR.density(pressure, temperature) + + self.assertGreater(HELIUM_PR.compressibility_factor(pressure, temperature), 1.0) + self.assertAlmostEqual( + HELIUM_PR.pressure_from_density(temperature, density), + pressure, + delta=pressure * 1.0e-12, + ) + + def test_air_reference_remains_available_for_other_models(self) -> None: + z = AIR_PR.compressibility_factor(101_325.0, 300.0) + density = AIR_PR.density(101_325.0, 300.0) + + self.assertAlmostEqual(z, 1.0, delta=1.0e-3) + self.assertAlmostEqual(density, 1.177, delta=0.01) + + def test_nitrogen_can_return_multiple_roots(self) -> None: + roots = NITROGEN_PR.compressibility_roots(1.0e6, 100.0) + + self.assertEqual(len(roots), 3) + self.assertLess(roots[0], roots[-1]) + self.assertEqual(NITROGEN_PR.compressibility_factor(1.0e6, 100.0, phase="liquid"), roots[0]) + self.assertEqual(NITROGEN_PR.compressibility_factor(1.0e6, 100.0, phase="vapor"), roots[-1]) + + def test_custom_fluid_can_be_defined_explicitly(self) -> None: + fluid = PengRobinsonFluid( + name="custom_helium", + molar_mass=0.004002602, + critical_temperature=5.1953, + critical_pressure=227_460.0, + acentric_factor=-0.385, + ) + + self.assertAlmostEqual(fluid.specific_gas_constant, 2077.3, delta=0.5) + self.assertAlmostEqual(fluid.compressibility_factor(1.0e6, 298.15), HELIUM_PR.compressibility_factor(1.0e6, 298.15)) + + def test_rejects_invalid_states(self) -> None: + with self.assertRaises(ValueError): + HELIUM_PR.compressibility_factor(0.0, 300.0) + with self.assertRaises(ValueError): + HELIUM_PR.density(101_325.0, 0.0) + with self.assertRaises(ValueError): + HELIUM_PR.pressure_from_density(300.0, -1.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_test_mql_config.py b/tests/test_test_mql_config.py index e3d61d9..d16ab58 100644 --- a/tests/test_test_mql_config.py +++ b/tests/test_test_mql_config.py @@ -16,6 +16,8 @@ class TestMqlConfigTest(unittest.TestCase): self.assertEqual(config.global_parameters["Pdq"], 1.0) self.assertEqual(config.global_parameters["V"], 15.0) self.assertEqual(config.global_parameters["cf"], 0.45) + self.assertEqual(config.fluid.name, "helium") + self.assertAlmostEqual(config.fluid.specific_gas_constant, 2077.3, delta=0.5) def test_resolves_amesim_parameter_expressions(self) -> None: variables = {"D2": 20.0, "P0": 153.0, "cf": 0.45}