diff --git a/AmesimModels/test_mql/README.md b/AmesimModels/test_mql/README.md index 5a263b3..81723c8 100644 --- a/AmesimModels/test_mql/README.md +++ b/AmesimModels/test_mql/README.md @@ -42,6 +42,11 @@ - 当前覆盖压缩因子、摩尔体积、密度和由密度反算压力。 - 已内置 `HELIUM_PR`，供 `test_mql` 默认使用。 +- `PythonModels/components/amesim_pneumatic.py` + - 提供 `test_mql` 后续会用到的 AMESim 气动组件原语。 + - 当前包含氦气 PR 压力闭合的气动容腔、标准可压缩孔口流量、单位换算 helper。 + - 这层仍是首版近似原语，后续必须通过 AMESim CSV 对齐修正。 + ## 物性约定 AMESim 模型中 `test_mql` 使用氦气，Python 侧当前通过 `HELIUM_PR` 使用 Peng-Robinson 状态方程计算气体压缩因子和密度。当前物性层先覆盖状态方程相关量，完整焓/内能偏差函数后续在接气室能量方程时再补。 @@ -88,7 +93,7 @@ AMESim 模型中 `test_mql` 使用氦气，Python 侧当前通过 `HELIUM_PR` ## 验证方式 ```bash -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 py_compile PythonModels/components/amesim_pneumatic.py 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/components/amesim_pneumatic.py b/PythonModels/components/amesim_pneumatic.py new file mode 100644 index 0000000..e49741e --- /dev/null +++ b/PythonModels/components/amesim_pneumatic.py @@ -0,0 +1,216 @@ +from __future__ import annotations + +from dataclasses import dataclass +from math import pi, sqrt + +from PythonModels.core.base import AlgebraicComponent, DynamicComponent +from PythonModels.core.medium import ThermodynamicProperties +from PythonModels.core.peng_robinson import HELIUM_PR, PengRobinsonFluid +from PythonModels.core.ports import PortState +from PythonModels.core.state import VolumeState + + +@dataclass(frozen=True) +class AmesimPneumaticGas: + """Caloric constants plus Peng-Robinson EOS for AMESim pneumatic components.""" + + fluid: PengRobinsonFluid = HELIUM_PR + cp: float = 5193.0 + cv: float = 3116.0 + + @property + def gamma(self) -> float: + return self.cp / self.cv + + @property + def R_gas(self) -> float: + return self.fluid.specific_gas_constant + + def density(self, pressure: float, temperature: float) -> float: + return self.fluid.density(pressure, temperature) + + def pressure(self, density: float, temperature: float) -> float: + return self.fluid.pressure_from_density(temperature, density) + + def specific_internal_energy(self, temperature: float) -> float: + return self.cv * temperature + + def specific_enthalpy(self, temperature: float) -> float: + return self.cp * temperature + + def temperature_from_internal_energy(self, specific_internal_energy: float) -> float: + if self.cv <= 0.0: + raise ValueError("cv must be positive.") + return specific_internal_energy / self.cv + + +HELIUM_PNEUMATIC_GAS = AmesimPneumaticGas() + + +def liters_to_m3(value: float) -> float: + return value * 1.0e-3 + + +def mm2_to_m2(value: float) -> float: + return value * 1.0e-6 + + +def diameter_mm_to_area_m2(diameter_mm: float) -> float: + diameter_m = diameter_mm * 1.0e-3 + return pi * diameter_m * diameter_m / 4.0 + + +class AmesimPneumaticVolume(DynamicComponent): + """First-pass AMESim pneumatic control volume using helium PR pressure closure.""" + + def __init__( + self, + name: str, + volume: float, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + p0: float = 101_325.0, + T0: float = 293.15, + ) -> None: + if volume <= 0.0: + raise ValueError("volume must be positive.") + super().__init__(name=name) + self.volume = volume + self.gas = gas + rho0 = gas.density(p0, T0) + m0 = rho0 * volume + U0 = m0 * gas.specific_internal_energy(T0) + self.state = VolumeState(m=m0, U=U0) + self.port_a = PortState() + + @classmethod + def from_liters( + cls, + name: str, + volume_liters: float, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + p0: float = 101_325.0, + T0: float = 293.15, + ) -> "AmesimPneumaticVolume": + return cls(name=name, volume=liters_to_m3(volume_liters), gas=gas, p0=p0, T0=T0) + + 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("volume mass must stay positive.") + T = self.gas.temperature_from_internal_energy(self.state.U / self.state.m) + rho = self.state.m / self.volume + p = self.gas.pressure(rho, T) + u = self.state.U / self.state.m + h = self.gas.specific_enthalpy(T) + self.port_a.p = p + self.port_a.h_outflow = h + return ThermodynamicProperties(p=p, T=T, rho=rho, u=u, h=h) + + def derivatives(self, inlet_h: float, m_flow: float) -> VolumeState: + return VolumeState(m=m_flow, U=m_flow * inlet_h) + + +class AmesimPneumaticOrifice(AlgebraicComponent): + """First-pass PNOR001/PNVO001-style compressible helium orifice. + + This is a calibrated placeholder boundary for the Python port. It preserves + AMESim-style area and coefficient inputs, but final parity must be checked + against AMESim CSV results before treating it as numerically equivalent. + """ + + def __init__( + self, + name: str, + area: float, + flow_coefficient: float = 1.0, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + opening: float = 1.0, + ) -> None: + if area < 0.0: + raise ValueError("area must be non-negative.") + if flow_coefficient < 0.0: + raise ValueError("flow_coefficient must be non-negative.") + super().__init__(name=name) + self.area = area + self.flow_coefficient = flow_coefficient + self.gas = gas + self.opening = opening + self.port_a = PortState() + self.port_b = PortState() + + @classmethod + def from_mm2( + cls, + name: str, + area_mm2: float, + flow_coefficient: float = 1.0, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, + opening: float = 1.0, + ) -> "AmesimPneumaticOrifice": + return cls( + name=name, + area=mm2_to_m2(area_mm2), + flow_coefficient=flow_coefficient, + gas=gas, + opening=opening, + ) + + @property + def effective_area(self) -> float: + return self.area * max(self.opening, 0.0) + + def mass_flow(self, p_a: float, p_b: float, upstream_temperature: float) -> float: + if p_a == p_b or self.effective_area == 0.0 or self.flow_coefficient == 0.0: + return 0.0 + if p_a > p_b: + return compressible_orifice_mass_flow( + upstream_pressure=p_a, + downstream_pressure=p_b, + upstream_temperature=upstream_temperature, + area=self.effective_area, + flow_coefficient=self.flow_coefficient, + gas=self.gas, + ) + return -compressible_orifice_mass_flow( + upstream_pressure=p_b, + downstream_pressure=p_a, + upstream_temperature=upstream_temperature, + area=self.effective_area, + flow_coefficient=self.flow_coefficient, + gas=self.gas, + ) + + +def compressible_orifice_mass_flow( + *, + upstream_pressure: float, + downstream_pressure: float, + upstream_temperature: float, + area: float, + flow_coefficient: float, + gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS, +) -> float: + if upstream_pressure <= 0.0 or downstream_pressure < 0.0: + raise ValueError("pressures must be non-negative and upstream pressure must be positive.") + if upstream_temperature <= 0.0: + raise ValueError("upstream_temperature must be positive.") + if area < 0.0 or flow_coefficient < 0.0: + raise ValueError("area and flow_coefficient must be non-negative.") + if downstream_pressure >= upstream_pressure or area == 0.0 or flow_coefficient == 0.0: + return 0.0 + + gamma = gas.gamma + pressure_ratio = max(downstream_pressure / upstream_pressure, 0.0) + critical_ratio = (2.0 / (gamma + 1.0)) ** (gamma / (gamma - 1.0)) + coefficient = flow_coefficient * area * upstream_pressure / sqrt(gas.R_gas * upstream_temperature) + if pressure_ratio <= critical_ratio: + flow_function = sqrt(gamma) * (2.0 / (gamma + 1.0)) ** ((gamma + 1.0) / (2.0 * (gamma - 1.0))) + else: + term = pressure_ratio ** (2.0 / gamma) - pressure_ratio ** ((gamma + 1.0) / gamma) + flow_function = sqrt((2.0 * gamma / (gamma - 1.0)) * max(term, 0.0)) + return coefficient * flow_function diff --git a/tests/test_amesim_pneumatic_components.py b/tests/test_amesim_pneumatic_components.py new file mode 100644 index 0000000..5c9243f --- /dev/null +++ b/tests/test_amesim_pneumatic_components.py @@ -0,0 +1,92 @@ +from __future__ import annotations + +import unittest + +from PythonModels.components.amesim_pneumatic import ( + HELIUM_PNEUMATIC_GAS, + AmesimPneumaticOrifice, + AmesimPneumaticVolume, + compressible_orifice_mass_flow, + diameter_mm_to_area_m2, + liters_to_m3, + mm2_to_m2, +) + + +class AmesimPneumaticComponentsTest(unittest.TestCase): + def test_unit_conversions(self) -> None: + self.assertAlmostEqual(liters_to_m3(15.0), 0.015) + self.assertAlmostEqual(mm2_to_m2(78.5), 78.5e-6) + self.assertAlmostEqual(diameter_mm_to_area_m2(10.0), 7.853981633974483e-5) + + def test_volume_initial_state_matches_requested_pressure_temperature(self) -> None: + volume = AmesimPneumaticVolume.from_liters( + name="pn_general_chamber", + volume_liters=57.0, + p0=15.3e6, + T0=293.15, + ) + props = volume.properties() + + self.assertAlmostEqual(props.p, 15.3e6, delta=15.3e6 * 1.0e-12) + self.assertAlmostEqual(props.T, 293.15) + self.assertGreater(props.rho, 20.0) + self.assertLess(props.rho, 30.0) + + def test_orifice_returns_signed_mass_flow(self) -> None: + orifice = AmesimPneumaticOrifice.from_mm2( + name="pn_orifice_18", + area_mm2=78.5, + flow_coefficient=0.9, + ) + + forward = orifice.mass_flow(15.3e6, 1.0e6, 293.15) + reverse = orifice.mass_flow(1.0e6, 15.3e6, 293.15) + + self.assertGreater(forward, 0.0) + self.assertAlmostEqual(reverse, -forward) + self.assertEqual(orifice.mass_flow(1.0e6, 1.0e6, 293.15), 0.0) + + def test_choked_flow_is_independent_of_lower_downstream_pressure(self) -> None: + base = compressible_orifice_mass_flow( + upstream_pressure=15.3e6, + downstream_pressure=1.0e6, + upstream_temperature=293.15, + area=78.5e-6, + flow_coefficient=0.9, + gas=HELIUM_PNEUMATIC_GAS, + ) + lower_back_pressure = compressible_orifice_mass_flow( + upstream_pressure=15.3e6, + downstream_pressure=0.1e6, + upstream_temperature=293.15, + area=78.5e-6, + flow_coefficient=0.9, + gas=HELIUM_PNEUMATIC_GAS, + ) + + self.assertGreater(base, 0.0) + self.assertAlmostEqual(lower_back_pressure, base) + + def test_subcritical_flow_decreases_as_back_pressure_rises(self) -> None: + low_back_pressure = compressible_orifice_mass_flow( + upstream_pressure=1.0e6, + downstream_pressure=0.6e6, + upstream_temperature=293.15, + area=78.5e-6, + flow_coefficient=0.9, + ) + high_back_pressure = compressible_orifice_mass_flow( + upstream_pressure=1.0e6, + downstream_pressure=0.9e6, + upstream_temperature=293.15, + area=78.5e-6, + flow_coefficient=0.9, + ) + + self.assertGreater(low_back_pressure, high_back_pressure) + self.assertGreater(high_back_pressure, 0.0) + + +if __name__ == "__main__": + unittest.main()