补充氦气Peng-Robinson物性库

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 .
 ```

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)

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,
         )
 

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()

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}