From 046aa49814177baad789f03d59d16ef87cbac356 Mon Sep 17 00:00:00 2001 From: huojiarong Date: Mon, 3 Aug 2026 15:34:33 +0000 Subject: [PATCH] =?UTF-8?q?=E5=AF=B9=E9=BD=90Amesim=E6=B0=A6=E6=B0=94PR?= =?UTF-8?q?=E7=89=A9=E6=80=A7=E4=B8=8EPNVO=E6=B5=81=E9=87=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../components/amesim/flow/orifices.py | 143 +++++++++--- .../components/amesim/flow/pipes.py | 26 ++- .../components/amesim/media/mediums.py | 207 +++++++++++++++++- .../components/amesim/storage/chambers.py | 8 +- app/simulation/core/medium.py | 37 ++++ app/simulation/core/peng_robinson.py | 138 +++++++++++- app/simulation/systems/generic.py | 7 +- tests/test_amesim_helium_medium.py | 74 ++++++- tests/test_amesim_pneumatic_components.py | 4 +- tests/test_amesim_pnl0001_pipe.py | 6 +- tests/test_amesim_pnvo001_fixed_component.py | 83 +++++++ tests/test_peng_robinson.py | 37 +++- ...est_pressure_flow_solver_initialization.py | 33 ++- 13 files changed, 730 insertions(+), 73 deletions(-) diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py index 689c73c..bdfe46b 100644 --- a/app/simulation/components/amesim/flow/orifices.py +++ b/app/simulation/components/amesim/flow/orifices.py @@ -206,9 +206,13 @@ class AmesimPnor001(AlgebraicComponent): def _upstream_temperature(self, port_name: str) -> float: port = self.get_port(port_name) - if port.h_outflow > 0.0: - return max(self.medium.temperature_from_enthalpy(port.h_outflow), 1.0) - return self.medium.T_ref + return max( + self.medium.temperature_from_pressure_enthalpy( + max(port.p, 1.0), + port.h_outflow, + ), + 1.0, + ) def mass_flow(self, p_1: float, p_2: float) -> float: if p_1 == p_2 or self.effective_area == 0.0: @@ -461,6 +465,7 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): self.port_2.h_outflow = initial_h self.port_3 = self.register_declared_port("port_3") self.port_3.h_outflow = initial_h + self._connected_h: dict[str, float] = {} @staticmethod def _integer_parameter(name: str, value: float) -> int: @@ -507,9 +512,16 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): def _upstream_temperature(self, port_name: str) -> float: port = self.get_port(port_name) - if port.h_outflow > 0.0: - return max(self.medium.temperature_from_enthalpy(port.h_outflow), 1.0) - return self.medium.T_ref + # A component port's h_outflow describes fluid leaving the valve; the + # upstream state comes from the connection on that same physical side. + inlet_h = self._connected_h.get(port_name, port.h_outflow) + return max( + self.medium.temperature_from_pressure_enthalpy( + max(port.p, 1.0), + inlet_h, + ), + 1.0, + ) def mass_flow(self, p_2: float, p_3: float) -> float: if p_2 == p_3 or self.effective_area == 0.0: @@ -526,6 +538,64 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): upstream_temperature=self._upstream_temperature("port_3"), ) + def _one_way_flow_characteristics( + self, + *, + upstream_pressure: float, + downstream_pressure: float, + upstream_temperature: float, + ) -> tuple[float, float]: + p_up = max(upstream_pressure, 1.0) + p_down = max(min(downstream_pressure, p_up), 0.0) + T_up = max(upstream_temperature, 1.0) + gamma_s = self.medium.isentropic_density_pressure_factor( + p_up, + T_up, + p_down, + ) + gamma_s = min(max(gamma_s, 1.0e-9), 1.0 - 1.0e-9) + density = max(self.medium.density(p_up, T_up), 1.0e-12) + pressure_ratio = max(p_down / p_up, 0.0) + critical_ratio = (2.0 * gamma_s / (gamma_s + 1.0)) ** ( + 1.0 / (1.0 - gamma_s) + ) + if 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)) + ** (gamma_s / (1.0 - gamma_s)) + ) + gas_velocity = sqrt( + 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, + ) + ) + gas_velocity = sqrt( + max( + 2.0 + / (1.0 - gamma_s) + * p_up + / density + * (1.0 - pressure_ratio ** (1.0 - gamma_s)), + 0.0, + ) + ) + return mass_flow_parameter, gas_velocity + def _one_way_mass_flow( self, *, @@ -534,44 +604,45 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): upstream_temperature: float, ) -> float: p_up = max(upstream_pressure, 1.0) - p_down = max(min(downstream_pressure, p_up), 0.0) T_up = max(upstream_temperature, 1.0) - gamma = max(self.medium.gamma, 1.000001) - pressure_ratio = max(p_down / p_up, 0.0) - critical_ratio = (2.0 / (gamma + 1.0)) ** (gamma / (gamma - 1.0)) - if pressure_ratio <= critical_ratio: - flow_factor = sqrt(gamma / (self.medium.R_gas * T_up)) * ( - 2.0 / (gamma + 1.0) - ) ** ((gamma + 1.0) / (2.0 * (gamma - 1.0))) - else: - expansion = pressure_ratio ** (2.0 / gamma) - pressure_ratio ** ( - (gamma + 1.0) / gamma - ) - flow_factor = sqrt( - max( - 2.0 - * gamma - * expansion - / (self.medium.R_gas * T_up * (gamma - 1.0)), - 0.0, - ) - ) - return self.effective_cq * self.effective_area * p_up * flow_factor + mass_flow_parameter, _gas_velocity = self._one_way_flow_characteristics( + upstream_pressure=p_up, + downstream_pressure=downstream_pressure, + upstream_temperature=T_up, + ) + return ( + self.effective_cq + * self.effective_area + * p_up + * mass_flow_parameter + / sqrt(T_up) + ) def component_result_values(self) -> Mapping[str, float]: p_2 = max(self.port_2.p, 1.0) p_3 = max(self.port_3.p, 1.0) - m_flow = abs(self.mass_flow(self.port_2.p, self.port_3.p)) - upstream_pressure = max(p_2, p_3) + if p_2 >= p_3: + upstream_port_name = "port_2" + upstream_pressure = p_2 + downstream_pressure = p_3 + flow_direction = 1.0 + else: + upstream_port_name = "port_3" + upstream_pressure = p_3 + downstream_pressure = p_2 + flow_direction = -1.0 upstream_temperature = self._upstream_temperature( - "port_2" if p_2 >= p_3 else "port_3" + upstream_port_name + ) + mass_flow_parameter, gas_velocity = self._one_way_flow_characteristics( + upstream_pressure=upstream_pressure, + downstream_pressure=downstream_pressure, + upstream_temperature=upstream_temperature, ) - density = max(self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12) - area = max(self.effective_area, 1.0e-18) return { "xv": self.opening, - "cm": m_flow / (self.effective_cq * area * upstream_pressure), - "gasvel": m_flow / (density * area), + "cm": mass_flow_parameter, + "gasvel": flow_direction * gas_velocity, } def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: @@ -605,6 +676,7 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent): ) def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: + self._connected_h = dict(connected_h) self.port_2.h_outflow = connected_h["port_3"] self.port_3.h_outflow = connected_h["port_2"] @@ -684,6 +756,7 @@ class AmesimPnvo001SignalOpening(AmesimPnvo001FixedOpening): self.port_2.h_outflow = initial_h self.port_3 = self.register_declared_port("port_3") self.port_3.h_outflow = initial_h + self._connected_h: dict[str, float] = {} @classmethod def create( diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index b445a72..625b320 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -164,9 +164,13 @@ class AmesimPnl00r(AlgebraicComponent): def _port_temperature(self, port_name: str) -> float: port = self.get_port(port_name) - if port.h_outflow > 0.0: - return max(self.medium.temperature_from_enthalpy(port.h_outflow), 1.0) - return self.medium.T_ref + return max( + self.medium.temperature_from_pressure_enthalpy( + max(port.p, 1.0), + port.h_outflow, + ), + 1.0, + ) def _dynamic_viscosity(self, temperature_k: float) -> float: return self.medium.dynamic_viscosity(temperature_k) @@ -494,9 +498,9 @@ class AmesimPnl0001(ThermodynamicVolumeComponent): self.volume = self.area * self.le self.exchange_area = pi * self.diam * self.le m0 = medium.density(self.p0, self.T0) * self.volume - U0 = m0 * medium.specific_internal_energy(self.T0) + U0 = m0 * medium.specific_internal_energy_at_pressure(self.p0, self.T0) self.state = VolumeState(m=m0, U=U0) - initial_h = medium.specific_enthalpy(self.T0) + initial_h = medium.specific_enthalpy_at_pressure(self.p0, self.T0) self.port_1 = self.register_declared_port("port_1") self.port_1.p = self.p0 self.port_1.h_outflow = initial_h @@ -978,8 +982,8 @@ class AmesimPnl0003(DynamicComponent): self.exchange_area = pi * self.diam * self.le self.state_1 = self._initial_state(float(p1_0), float(T1_0)) self.state_2 = self._initial_state(float(p2_0), float(T2_0)) - h1 = medium.specific_enthalpy(float(T1_0)) - h2 = medium.specific_enthalpy(float(T2_0)) + h1 = medium.specific_enthalpy_at_pressure(float(p1_0), float(T1_0)) + h2 = medium.specific_enthalpy_at_pressure(float(p2_0), float(T2_0)) self.port_1 = self.register_declared_port("port_1") self.port_1.p = float(p1_0) self.port_1.h_outflow = h1 @@ -999,7 +1003,13 @@ class AmesimPnl0003(DynamicComponent): def _initial_state(self, pressure: float, temperature: float) -> VolumeState: mass = self.medium.density(pressure, temperature) * self.compliance_volume - return VolumeState(m=mass, U=mass * self.medium.specific_internal_energy(temperature)) + return VolumeState( + m=mass, + U=mass * self.medium.specific_internal_energy_at_pressure( + pressure, + temperature, + ), + ) def get_state_vector(self) -> list[float]: return [*self.state_1.as_vector(), *self.state_2.as_vector()] diff --git a/app/simulation/components/amesim/media/mediums.py b/app/simulation/components/amesim/media/mediums.py index 69bc002..7ef9290 100644 --- a/app/simulation/components/amesim/media/mediums.py +++ b/app/simulation/components/amesim/media/mediums.py @@ -4,7 +4,11 @@ from collections.abc import Callable from dataclasses import dataclass from typing import ClassVar -from app.simulation.core.medium import GasMedium, IdealGasMedium +from app.simulation.core.medium import ( + GasMedium, + IdealGasMedium, + ThermodynamicProperties, +) from app.simulation.core.peng_robinson import HELIUM_PR, PengRobinsonFluid @@ -36,18 +40,19 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): """AMESim helium with a Peng-Robinson mechanical equation of state. The pressure-density-temperature relation is evaluated by the shared - ``HELIUM_PR`` fluid. The first public AMESim port keeps the committed - constant-heat-capacity caloric model so it can be consumed through the - same :class:`GasMedium` contract as ideal-gas air. + ``HELIUM_PR`` fluid. The caloric reference follows the constant NASA + polynomial from Simcenter Amesim 2404 ``helium_cp_h_s.data``. """ SUBSTANCE_ID: ClassVar[str] = "helium" PROPERTY_METHOD_ID: ClassVar[str] = "peng_robinson" fluid: ClassVar[PengRobinsonFluid] = HELIUM_PR + nasa_cp_over_R: ClassVar[float] = 2.5 + nasa_enthalpy_constant_K: ClassVar[float] = -745.375 name: str = "AMESimHeliumPengRobinson" R_gas: float = HELIUM_PR.specific_gas_constant - cp_ref: float = 5193.0 + cp_ref: float = nasa_cp_over_R * HELIUM_PR.specific_gas_constant T_ref: float = 293.15 cp_slope: float = 0.0 viscosity_ref: float = 1.96e-5 @@ -56,7 +61,7 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): @property def cv(self) -> float: - return 3116.0 + return (self.nasa_cp_over_R - 1.0) * self.R_gas def cv_at_temperature(self, T: float) -> float: del T @@ -65,11 +70,201 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): def density(self, p: float, T: float) -> float: return self.fluid.density(p, T) + def _real_heat_capacities( + self, + p: float, + T: float, + ) -> tuple[float, float, float, float, float]: + density = self.density(p, T) + pressure_density_derivative = ( + self.fluid.pressure_density_derivative_at_temperature( + T, + density, + ) + ) + pressure_temperature_derivative = ( + self.fluid.pressure_temperature_derivative_at_density( + T, + density, + ) + ) + cv = ( + self.cv_at_temperature(T) + + self.fluid.residual_isochoric_heat_capacity_at_density(T, density) + ) + cp = ( + cv + + T + * pressure_temperature_derivative + * pressure_temperature_derivative + / (density * density * pressure_density_derivative) + ) + if cp <= 0.0 or cv <= 0.0: + raise ValueError("Real-gas heat capacities must be positive.") + return ( + cp, + cv, + density, + pressure_density_derivative, + pressure_temperature_derivative, + ) + + def _local_isentropic_density_pressure_factor( + self, + p: float, + T: float, + ) -> tuple[float, float]: + cp, cv, density, pressure_density_derivative, pressure_temperature_derivative = ( + self._real_heat_capacities(p, T) + ) + heat_capacity_ratio = cp / cv + factor = p / ( + density * pressure_density_derivative * heat_capacity_ratio + ) + exponent = ( + p + * (heat_capacity_ratio - 1.0) + / ( + heat_capacity_ratio + * T + * pressure_temperature_derivative + ) + ) + return factor, exponent + + def isentropic_density_pressure_factor( + self, + p: float, + T: float, + downstream_pressure: float | None = None, + ) -> float: + upstream_factor, isentropic_temperature_exponent = ( + self._local_isentropic_density_pressure_factor(p, T) + ) + if downstream_pressure is None or downstream_pressure >= p: + return upstream_factor + + pressure_ratio = max(downstream_pressure / p, 1.0e-12) + isentropic_temperature = max( + T * pressure_ratio**isentropic_temperature_exponent, + 2.2, + ) + downstream_factor, _unused_exponent = ( + self._local_isentropic_density_pressure_factor( + max(downstream_pressure, 1.0), + isentropic_temperature, + ) + ) + # AMESim 2404 saggs_ evaluates the local factor at the upstream + # state and at an approximate isentropic downstream state. + return 0.5 * (upstream_factor + downstream_factor) + def pressure(self, m: float, T: float, V: float) -> float: if V <= 0.0: raise ValueError("Volume must stay positive.") return self.fluid.pressure_from_density(T, m / V) + def specific_internal_energy(self, T: float) -> float: + return self.R_gas * ( + (self.nasa_cp_over_R - 1.0) * T + + self.nasa_enthalpy_constant_K + ) + + def specific_internal_energy_at_pressure(self, p: float, T: float) -> float: + density = self.density(p, T) + return ( + self.specific_internal_energy(T) + + self.fluid.residual_specific_internal_energy_at_density(T, density) + ) + + def specific_enthalpy(self, T: float) -> float: + return self.R_gas * ( + self.nasa_cp_over_R * T + + self.nasa_enthalpy_constant_K + ) + + def specific_enthalpy_at_pressure(self, p: float, T: float) -> float: + return self.specific_enthalpy(T) + self.fluid.residual_specific_enthalpy(p, T) + + def temperature_from_internal_energy(self, u: float) -> float: + return ( + u / self.R_gas - self.nasa_enthalpy_constant_K + ) / (self.nasa_cp_over_R - 1.0) + + def temperature_from_enthalpy(self, h: float) -> float: + return ( + h / self.R_gas - self.nasa_enthalpy_constant_K + ) / self.nasa_cp_over_R + + def temperature_from_pressure_enthalpy(self, p: float, h: float) -> float: + temperature = max(self.temperature_from_enthalpy(h), 2.2) + for _iteration in range(16): + residual_enthalpy = self.fluid.residual_specific_enthalpy(p, temperature) + next_temperature = max( + self.temperature_from_enthalpy(h - residual_enthalpy), + 2.2, + ) + if abs(next_temperature - temperature) <= 1.0e-10 * max( + temperature, + 1.0, + ): + return next_temperature + temperature = next_temperature + return temperature + + def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: + if m <= 0.0: + raise ValueError("Mass must stay positive when recovering temperature.") + return self.temperature_from_internal_energy(U / m) + + def properties_from_mU( + self, + m: float, + U: float, + V: float, + ) -> ThermodynamicProperties: + if m <= 0.0: + raise ValueError("Mass must stay positive when recovering temperature.") + if V <= 0.0: + raise ValueError("Volume must stay positive.") + density = m / V + target_internal_energy = U / m + temperature = max( + self.temperature_from_internal_energy(target_internal_energy), + 2.2, + ) + for _iteration in range(16): + residual_internal_energy = ( + self.fluid.residual_specific_internal_energy_at_density( + temperature, + density, + ) + ) + next_temperature = max( + self.temperature_from_internal_energy( + target_internal_energy - residual_internal_energy + ), + 2.2, + ) + if abs(next_temperature - temperature) <= 1.0e-10 * max( + temperature, + 1.0, + ): + temperature = next_temperature + break + temperature = next_temperature + pressure = self.fluid.pressure_from_density(temperature, density) + return ThermodynamicProperties( + p=pressure, + T=temperature, + rho=density, + u=target_internal_energy, + h=self.specific_enthalpy_at_pressure( + pressure, + temperature, + ), + ) + @dataclass(frozen=True) class AmesimGasPropertyModelSpec: diff --git a/app/simulation/components/amesim/storage/chambers.py b/app/simulation/components/amesim/storage/chambers.py index 1bcf5c7..0be73d0 100644 --- a/app/simulation/components/amesim/storage/chambers.py +++ b/app/simulation/components/amesim/storage/chambers.py @@ -142,9 +142,9 @@ class AmesimPnch023(ThermodynamicVolumeComponent): self.p0 = float(p0) self.T0 = float(T0) m0 = medium.density(self.p0, self.T0) * self.cvol - U0 = m0 * medium.specific_internal_energy(self.T0) + U0 = m0 * medium.specific_internal_energy_at_pressure(self.p0, self.T0) self.state = VolumeState(m=m0, U=U0) - initial_h = medium.specific_enthalpy(self.T0) + initial_h = medium.specific_enthalpy_at_pressure(self.p0, self.T0) self.port_1 = self.register_declared_port("port_1") self.port_1.p = self.p0 self.port_1.h_outflow = initial_h @@ -413,9 +413,9 @@ class AmesimPnch012(ThermodynamicVolumeComponent): if self.total_volume() <= 0.0: raise ValueError("PNCH012 total volume must be positive.") m0 = medium.density(self.p0, self.T0) * self.total_volume() - U0 = m0 * medium.specific_internal_energy(self.T0) + U0 = m0 * medium.specific_internal_energy_at_pressure(self.p0, self.T0) self.state = VolumeState(m=m0, U=U0) - initial_h = medium.specific_enthalpy(self.T0) + initial_h = medium.specific_enthalpy_at_pressure(self.p0, self.T0) for port_name in ("port_1", "port_2", "port_3", "port_4"): port = self.register_declared_port(port_name) port.p = self.p0 diff --git a/app/simulation/core/medium.py b/app/simulation/core/medium.py index c18541c..9108ec4 100644 --- a/app/simulation/core/medium.py +++ b/app/simulation/core/medium.py @@ -38,16 +38,29 @@ class GasMedium(Protocol): def density(self, p: float, T: float) -> float: ... + def isentropic_density_pressure_factor( + self, + p: float, + T: float, + downstream_pressure: float | None = None, + ) -> float: ... + def dynamic_viscosity(self, T: float) -> float: ... def specific_internal_energy(self, T: float) -> float: ... + def specific_internal_energy_at_pressure(self, p: float, T: float) -> float: ... + def specific_enthalpy(self, T: float) -> float: ... + def specific_enthalpy_at_pressure(self, p: float, T: float) -> float: ... + def temperature_from_internal_energy(self, u: float) -> float: ... def temperature_from_enthalpy(self, h: float) -> float: ... + def temperature_from_pressure_enthalpy(self, p: float, h: float) -> float: ... + def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: ... def pressure(self, m: float, T: float, V: float) -> float: ... @@ -96,6 +109,18 @@ class IdealGasMedium: def density(self, p: float, T: float) -> float: return p / (self.R_gas * T) + def isentropic_density_pressure_factor( + self, + p: float, + T: float, + downstream_pressure: float | None = None, + ) -> float: + del p + del downstream_pressure + cp = self.cp_at_temperature(T) + cv = self.cv_at_temperature(T) + return cv / cp + def dynamic_viscosity(self, T: float) -> float: """Return dynamic viscosity using the default air Sutherland law.""" @@ -116,6 +141,10 @@ class IdealGasMedium: + 0.5 * self.cp_slope * delta_T * delta_T ) + def specific_internal_energy_at_pressure(self, p: float, T: float) -> float: + del p + return self.specific_internal_energy(T) + def specific_enthalpy(self, T: float) -> float: delta_T = T - self.T_ref return ( @@ -124,6 +153,10 @@ class IdealGasMedium: + 0.5 * self.cp_slope * delta_T * delta_T ) + def specific_enthalpy_at_pressure(self, p: float, T: float) -> float: + del p + return self.specific_enthalpy(T) + def temperature_from_internal_energy(self, u: float) -> float: reference_internal_energy = self.cv * self.T_ref delta_u = u - reference_internal_energy @@ -156,6 +189,10 @@ class IdealGasMedium: delta_T = positive_root if abs(positive_root) <= abs(negative_root) else negative_root return self.T_ref + delta_T + def temperature_from_pressure_enthalpy(self, p: float, h: float) -> float: + del p + return self.temperature_from_enthalpy(h) + def temperature_from_mass_internal_energy(self, m: float, U: float) -> float: if m <= 0.0: raise ValueError("Mass must stay positive when recovering temperature.") diff --git a/app/simulation/core/peng_robinson.py b/app/simulation/core/peng_robinson.py index 51ea44f..2365a0f 100644 --- a/app/simulation/core/peng_robinson.py +++ b/app/simulation/core/peng_robinson.py @@ -6,6 +6,10 @@ from dataclasses import dataclass from math import acos, cos, isfinite, log, pi, sqrt UNIVERSAL_GAS_CONSTANT = 8.31446261815324 +# Simcenter Amesim 2404 ``sag_reinit_eos_`` keeps more digits than the +# commonly printed Peng-Robinson constants 0.45724 and 0.07780. +PENG_ROBINSON_A_COEFFICIENT = 0.457235583 +PENG_ROBINSON_B_COEFFICIENT = 0.07779607 @dataclass(frozen=True) @@ -29,7 +33,7 @@ class PengRobinsonFluid: @property def a_parameter(self) -> float: return ( - 0.45724 + PENG_ROBINSON_A_COEFFICIENT * UNIVERSAL_GAS_CONSTANT * UNIVERSAL_GAS_CONSTANT * self.critical_temperature @@ -39,7 +43,12 @@ class PengRobinsonFluid: @property def b_parameter(self) -> float: - return 0.07780 * UNIVERSAL_GAS_CONSTANT * self.critical_temperature / self.critical_pressure + return ( + PENG_ROBINSON_B_COEFFICIENT + * UNIVERSAL_GAS_CONSTANT + * self.critical_temperature + / self.critical_pressure + ) @property def kappa(self) -> float: @@ -62,12 +71,32 @@ class PengRobinsonFluid: / (self.critical_temperature * sqrt_reduced_temperature) ) + def alpha_temperature_second_derivative(self, temperature: float) -> float: + self._validate_temperature(temperature) + reduced_temperature = temperature / self.critical_temperature + sqrt_reduced_temperature = sqrt(reduced_temperature) + alpha_base = 1.0 + self.kappa * (1.0 - sqrt_reduced_temperature) + return ( + self.kappa + / (2.0 * self.critical_temperature * self.critical_temperature) + * ( + self.kappa / reduced_temperature + + alpha_base / (reduced_temperature * sqrt_reduced_temperature) + ) + ) + def attractive_parameter(self, temperature: float) -> float: return self.a_parameter * self.alpha(temperature) def attractive_parameter_temperature_derivative(self, temperature: float) -> float: return self.a_parameter * self.alpha_temperature_derivative(temperature) + def attractive_parameter_temperature_second_derivative( + self, + temperature: float, + ) -> float: + return self.a_parameter * self.alpha_temperature_second_derivative(temperature) + def pressure_from_molar_volume(self, temperature: float, molar_volume: float) -> float: self._validate_temperature(temperature) if molar_volume <= self.b_parameter: @@ -83,6 +112,51 @@ class PengRobinsonFluid: raise ValueError("Density must be positive.") return self.pressure_from_molar_volume(temperature, self.molar_mass / density) + def pressure_temperature_derivative_at_density( + self, + temperature: float, + density: float, + ) -> float: + self._validate_temperature(temperature) + if density <= 0.0: + raise ValueError("Density must be positive.") + molar_volume = self.molar_mass / density + if molar_volume <= self.b_parameter: + raise RecoverableTrialStateError( + "Molar volume must be larger than Peng-Robinson b parameter." + ) + b = self.b_parameter + denominator = molar_volume * (molar_volume + b) + b * (molar_volume - b) + return ( + UNIVERSAL_GAS_CONSTANT / (molar_volume - b) + - self.attractive_parameter_temperature_derivative(temperature) / denominator + ) + + def pressure_density_derivative_at_temperature( + self, + temperature: float, + density: float, + ) -> float: + self._validate_temperature(temperature) + if density <= 0.0: + raise ValueError("Density must be positive.") + molar_volume = self.molar_mass / density + if molar_volume <= self.b_parameter: + raise RecoverableTrialStateError( + "Molar volume must be larger than Peng-Robinson b parameter." + ) + b = self.b_parameter + denominator = molar_volume * (molar_volume + b) + b * (molar_volume - b) + pressure_molar_volume_derivative = ( + -UNIVERSAL_GAS_CONSTANT * temperature / (molar_volume - b) ** 2 + + self.attractive_parameter(temperature) + * 2.0 + * (molar_volume + b) + / denominator**2 + ) + molar_volume_density_derivative = -self.molar_mass / (density * density) + return pressure_molar_volume_derivative * molar_volume_density_derivative + def reduced_parameters(self, pressure: float, temperature: float) -> tuple[float, float]: self._validate_pressure_temperature(pressure, temperature) a_alpha = self.attractive_parameter(temperature) @@ -165,6 +239,63 @@ class PengRobinsonFluid: ) return residual_molar_enthalpy / self.molar_mass + def residual_specific_internal_energy_at_density( + self, + temperature: float, + density: float, + ) -> float: + """Return Peng-Robinson internal-energy departure, J/kg.""" + self._validate_temperature(temperature) + if density <= 0.0: + raise ValueError("Density must be positive.") + molar_volume = self.molar_mass / density + b = self.b_parameter + if molar_volume <= b: + raise RecoverableTrialStateError( + "Molar volume must be larger than Peng-Robinson b parameter." + ) + attractive = self.attractive_parameter(temperature) + d_attractive_d_temperature = ( + self.attractive_parameter_temperature_derivative(temperature) + ) + log_argument = ( + molar_volume + (1.0 + sqrt(2.0)) * b + ) / ( + molar_volume + (1.0 - sqrt(2.0)) * b + ) + residual_molar_internal_energy = ( + temperature * d_attractive_d_temperature - attractive + ) * log(log_argument) / (2.0 * sqrt(2.0) * b) + return residual_molar_internal_energy / self.molar_mass + + def residual_isochoric_heat_capacity_at_density( + self, + temperature: float, + density: float, + ) -> float: + """Return the constant-volume heat-capacity departure, J/kg/K.""" + self._validate_temperature(temperature) + if density <= 0.0: + raise ValueError("Density must be positive.") + molar_volume = self.molar_mass / density + b = self.b_parameter + if molar_volume <= b: + raise RecoverableTrialStateError( + "Molar volume must be larger than Peng-Robinson b parameter." + ) + log_argument = ( + molar_volume + (1.0 + sqrt(2.0)) * b + ) / ( + molar_volume + (1.0 - sqrt(2.0)) * b + ) + residual_molar_cv = ( + temperature + * self.attractive_parameter_temperature_second_derivative(temperature) + * log(log_argument) + / (2.0 * sqrt(2.0) * b) + ) + return residual_molar_cv / self.molar_mass + @staticmethod def _validate_temperature(temperature: float) -> None: if temperature <= 0.0: @@ -181,7 +312,8 @@ HELIUM_PR = PengRobinsonFluid( molar_mass=0.004002602, critical_temperature=5.1953, critical_pressure=227_460.0, - acentric_factor=-0.385, + # Simcenter Amesim 2404 helium_eos.data. + acentric_factor=-0.382, ) NITROGEN_PR = PengRobinsonFluid( diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py index 47a02e1..0871ca7 100644 --- a/app/simulation/systems/generic.py +++ b/app/simulation/systems/generic.py @@ -268,8 +268,13 @@ class GenericFluidSystem: component.refresh_thermodynamic_ports() algebraic = self.pressure_flow_solver.solve() stream, connected_h = self.stream_resolver.solve() + # Some constitutive flow laws recover their upstream temperature from + # the connected stream enthalpy. Stream propagation updates that + # cache after the first pressure-flow pass, so refresh explicit flows + # once more before evaluating state derivatives and result variables. + algebraic = self.pressure_flow_solver.solve() self.mechanical_state_reducer.update_constraint_accelerations() - self.algebraic_solve_count += 1 + int(bool(pneumatic_volume.propagated)) + self.algebraic_solve_count += 2 + int(bool(pneumatic_volume.propagated)) self.max_algebraic_residual = max( self.max_algebraic_residual, algebraic.max_scaled_residual, diff --git a/tests/test_amesim_helium_medium.py b/tests/test_amesim_helium_medium.py index 0770bf8..f909368 100644 --- a/tests/test_amesim_helium_medium.py +++ b/tests/test_amesim_helium_medium.py @@ -18,7 +18,7 @@ from tests.test_system_xml_protocol import physical_port class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): - def test_uses_shared_peng_robinson_eos_and_committed_caloric_constants( + def test_uses_amesim_2404_peng_robinson_and_nasa_constants( self, ) -> None: medium = AmesimHeliumPengRobinsonMedium() @@ -32,11 +32,11 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): self.assertEqual(medium.SUBSTANCE_ID, "helium") self.assertEqual(medium.PROPERTY_METHOD_ID, "peng_robinson") self.assertAlmostEqual(medium.R_gas, HELIUM_PR.specific_gas_constant) - self.assertEqual(medium.cp_ref, 5193.0) - self.assertEqual(medium.cv, 3116.0) - self.assertEqual(medium.cp_at_temperature(400.0), 5193.0) - self.assertEqual(medium.cv_at_temperature(400.0), 3116.0) - self.assertAlmostEqual(medium.gamma, 5193.0 / 3116.0) + self.assertEqual(medium.cp_ref, 2.5 * medium.R_gas) + self.assertEqual(medium.cv, 1.5 * medium.R_gas) + self.assertEqual(medium.cp_at_temperature(400.0), 2.5 * medium.R_gas) + self.assertEqual(medium.cv_at_temperature(400.0), 1.5 * medium.R_gas) + self.assertAlmostEqual(medium.gamma, 5.0 / 3.0) self.assertAlmostEqual(density, HELIUM_PR.density(pressure, temperature)) self.assertAlmostEqual( medium.pressure(density * volume, temperature, volume), @@ -44,7 +44,7 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): delta=pressure * 1.0e-12, ) - def test_constant_caloric_model_and_sutherland_viscosity_complete_contract( + def test_nasa_caloric_reference_and_sutherland_viscosity_complete_contract( self, ) -> None: medium = AmesimHeliumPengRobinsonMedium() @@ -53,8 +53,14 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): internal_energy = medium.specific_internal_energy(temperature) enthalpy = medium.specific_enthalpy(temperature) - self.assertEqual(internal_energy, 3116.0 * temperature) - self.assertEqual(enthalpy, 5193.0 * temperature) + self.assertEqual( + internal_energy, + medium.R_gas * (1.5 * temperature - 745.375), + ) + self.assertEqual( + enthalpy, + medium.R_gas * (2.5 * temperature - 745.375), + ) self.assertEqual( medium.temperature_from_internal_energy(internal_energy), temperature, @@ -66,6 +72,56 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): self.assertAlmostEqual(medium.dynamic_viscosity(293.15), 1.96e-5) self.assertEqual(medium.sutherland_constant, 79.4) + def test_pressure_transport_enthalpy_includes_peng_robinson_departure( + self, + ) -> None: + medium = AmesimHeliumPengRobinsonMedium() + pressure = 13_839_965.0 + temperature = 287.7322 + + ideal_enthalpy = medium.specific_enthalpy(temperature) + transport_enthalpy = medium.specific_enthalpy_at_pressure(pressure, temperature) + + self.assertGreater( + transport_enthalpy, + ideal_enthalpy, + ) + self.assertAlmostEqual( + medium.temperature_from_pressure_enthalpy(pressure, transport_enthalpy), + temperature, + delta=1.0e-8, + ) + + density = medium.density(pressure, temperature) + volume = 0.01 + properties = medium.properties_from_mU( + density * volume, + density + * volume + * medium.specific_internal_energy_at_pressure(pressure, temperature), + volume, + ) + self.assertAlmostEqual(properties.p, pressure, delta=pressure * 1.0e-10) + self.assertAlmostEqual(properties.T, temperature, delta=1.0e-8) + self.assertAlmostEqual( + properties.h, + transport_enthalpy, + delta=2.0e-8, + ) + + def test_amesim_2404_real_gas_isentropic_factor_reference(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + + self.assertAlmostEqual( + medium.isentropic_density_pressure_factor( + 15_201_996.778497815, + 292.3997890954483, + 405_072.8123335606, + ), + 0.5807873871273339, + delta=1.0e-12, + ) + def test_definition_maps_local_selector_to_amesim_fluid_and_eos_codes( self, ) -> None: diff --git a/tests/test_amesim_pneumatic_components.py b/tests/test_amesim_pneumatic_components.py index 08cd201..3d0c748 100644 --- a/tests/test_amesim_pneumatic_components.py +++ b/tests/test_amesim_pneumatic_components.py @@ -51,8 +51,8 @@ class AmesimPneumaticComponentsTest(unittest.TestCase): self.assertAlmostEqual(ideal_reference_h, -54099.6354, delta=0.001) self.assertAlmostEqual( pressure_reference_h - ideal_reference_h, - 11936.1, - delta=0.1, + 12762.690087540526, + delta=1.0e-4, ) def test_pressure_transport_enthalpy_restores_absolute_energy_offset(self) -> None: diff --git a/tests/test_amesim_pnl0001_pipe.py b/tests/test_amesim_pnl0001_pipe.py index 0243547..01a9a8a 100644 --- a/tests/test_amesim_pnl0001_pipe.py +++ b/tests/test_amesim_pnl0001_pipe.py @@ -23,7 +23,11 @@ class AmesimPnl0001PipeTests(unittest.TestCase): self.assertEqual(len(self.pipe.get_state_vector()), 2) self.assertAlmostEqual(properties.p, 15.3e6, delta=1.0e-5) self.assertAlmostEqual(properties.T, 293.15) - self.assertAlmostEqual(self.pipe.gas_mass_g(), 3.716965219188, places=10) + self.assertAlmostEqual( + self.pipe.gas_mass_g(), + self.pipe.gas.density(15.3e6, 293.15) * self.pipe.volume * 1000.0, + places=10, + ) def test_resistance_flow_follows_pressure_gradient(self) -> None: forward = self.pipe.resistance_mass_flow( diff --git a/tests/test_amesim_pnvo001_fixed_component.py b/tests/test_amesim_pnvo001_fixed_component.py index 6f5d0c5..3b51c34 100644 --- a/tests/test_amesim_pnvo001_fixed_component.py +++ b/tests/test_amesim_pnvo001_fixed_component.py @@ -3,6 +3,9 @@ from __future__ import annotations import unittest from app.simulation.components.amesim.flow.orifices import AmesimPnvo001FixedOpening +from app.simulation.components.amesim.media.mediums import ( + AmesimHeliumPengRobinsonMedium, +) from app.simulation.core.medium import IdealGasMedium from app.simulation.registry import COMPONENT_MODEL_REGISTRY @@ -61,6 +64,86 @@ class AmesimPnvo001FixedOpeningComponentTests(unittest.TestCase): self.assertLess(reverse, 0.0) self.assertAlmostEqual(forward, -reverse) + def test_mass_flow_uses_connected_enthalpy_from_the_upstream_side(self) -> None: + valve = AmesimPnvo001FixedOpening("valve_1", self.medium, opening=0.5) + hot_h = self.medium.specific_enthalpy(600.0) + cold_h = self.medium.specific_enthalpy(200.0) + + valve.update_stream_outflows({"port_2": hot_h, "port_3": cold_h}) + + self.assertEqual(valve.port_2.h_outflow, cold_h) + self.assertEqual(valve.port_3.h_outflow, hot_h) + self.assertAlmostEqual( + valve.mass_flow(500000.0, 100000.0), + valve._one_way_mass_flow( + upstream_pressure=500000.0, + downstream_pressure=100000.0, + upstream_temperature=600.0, + ), + ) + self.assertAlmostEqual( + valve.mass_flow(100000.0, 500000.0), + -valve._one_way_mass_flow( + upstream_pressure=500000.0, + downstream_pressure=100000.0, + upstream_temperature=200.0, + ), + ) + + + def test_helium_flow_uses_pressure_enthalpy_and_real_gas_factor(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + valve = AmesimPnvo001FixedOpening( + "valve_1", + medium, + cq=0.45, + area0=78.5e-6, + opening=1.0, + ) + upstream_pressure = 15_201_996.778497815 + downstream_pressure = 405_072.8123335606 + upstream_temperature = 292.3997890954483 + valve.port_2.p = upstream_pressure + valve.port_3.p = downstream_pressure + upstream_enthalpy = medium.specific_enthalpy_at_pressure( + upstream_pressure, + upstream_temperature, + ) + valve.update_stream_outflows( + { + "port_2": upstream_enthalpy, + "port_3": medium.specific_enthalpy_at_pressure( + downstream_pressure, + 393.46713105173205, + ), + } + ) + + self.assertLess(upstream_enthalpy, 0.0) + self.assertAlmostEqual( + valve._upstream_temperature("port_2"), + upstream_temperature, + delta=1.0e-8, + ) + + self.assertAlmostEqual( + valve.mass_flow(upstream_pressure, downstream_pressure), + 0.49551308906777447, + delta=1.0e-9, + ) + results = valve.component_result_values() + + self.assertAlmostEqual( + results["cm"], + 0.015778323343746598, + delta=1.0e-11, + ) + self.assertAlmostEqual( + results["gasvel"], + 894.7011595213406, + delta=1.0e-8, + ) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_peng_robinson.py b/tests/test_peng_robinson.py index 3d85626..4312219 100644 --- a/tests/test_peng_robinson.py +++ b/tests/test_peng_robinson.py @@ -19,6 +19,11 @@ class PengRobinsonTest(unittest.TestCase): density = HELIUM_PR.density(pressure, temperature) self.assertGreater(HELIUM_PR.compressibility_factor(pressure, temperature), 1.0) + self.assertAlmostEqual( + density, + 24.114225444477153, + delta=1.0e-12, + ) self.assertAlmostEqual( HELIUM_PR.pressure_from_density(temperature, density), pressure, @@ -28,17 +33,43 @@ class PengRobinsonTest(unittest.TestCase): def test_helium_residual_enthalpy_is_small_at_atmosphere(self) -> None: self.assertAlmostEqual( HELIUM_PR.residual_specific_enthalpy(101_300.0, 298.15), - 17.834, + 25.3556466754, delta=0.01, ) def test_helium_residual_enthalpy_captures_high_pressure_departure(self) -> None: self.assertAlmostEqual( HELIUM_PR.residual_specific_enthalpy(13_839_965.0, 287.7322), - 11953.9, + 12788.0457342, delta=0.1, ) + def test_residual_isochoric_heat_capacity_is_internal_energy_derivative( + self, + ) -> None: + temperature = 292.3997890954483 + density = 24.02697515640726 + temperature_step = 1.0e-3 + numerical_derivative = ( + HELIUM_PR.residual_specific_internal_energy_at_density( + temperature + temperature_step, + density, + ) + - HELIUM_PR.residual_specific_internal_energy_at_density( + temperature - temperature_step, + density, + ) + ) / (2.0 * temperature_step) + + self.assertAlmostEqual( + HELIUM_PR.residual_isochoric_heat_capacity_at_density( + temperature, + density, + ), + numerical_derivative, + delta=1.0e-6, + ) + 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) @@ -60,7 +91,7 @@ class PengRobinsonTest(unittest.TestCase): molar_mass=0.004002602, critical_temperature=5.1953, critical_pressure=227_460.0, - acentric_factor=-0.385, + acentric_factor=-0.382, ) self.assertAlmostEqual(fluid.specific_gas_constant, 2077.3, delta=0.5) diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py index b20b4ed..bc844d6 100644 --- a/tests/test_pressure_flow_solver_initialization.py +++ b/tests/test_pressure_flow_solver_initialization.py @@ -16,6 +16,7 @@ from app.simulation.components.experimental.storage.tank import Tank from app.simulation.core.medium import IdealGasMedium from app.simulation.core.state import VolumeState from app.simulation.solvers.algebraic import AlgebraicSolveError, PressureFlowSolver +from app.simulation.systems.generic import GenericFluidSystem from app.simulation.systems.network import SimulationNetwork @@ -59,10 +60,40 @@ class PressureFlowSolverInitializationTests(unittest.TestCase): mass = medium.density(pressure, temperature) * component.V component.state = VolumeState( m=mass, - U=mass * medium.specific_internal_energy(temperature), + U=mass + * medium.specific_internal_energy_at_pressure( + pressure, + temperature, + ), ) component.refresh_thermodynamic_ports() + def test_stream_enthalpy_refreshes_explicit_orifice_flow(self) -> None: + network, medium, high, low, valve = self._near_equal_pressure_network() + self._set_pressure_temperature( + high, + medium, + 15_201_996.778497815, + 292.3997890954483, + ) + self._set_pressure_temperature( + low, + medium, + 405_072.8123335606, + 393.46713105173205, + ) + + system = GenericFluidSystem(network) + system.consistent_initial_state_vector() + + self.assertLess(valve._connected_h["port_2"], 0.0) + self.assertAlmostEqual( + valve.port_2.m_flow, + valve.mass_flow(valve.port_2.p, valve.port_3.p), + places=12, + ) + self.assertAlmostEqual(valve.port_3.m_flow, -valve.port_2.m_flow, places=12) + def test_current_storage_pressure_reseeds_stale_orifice_ports_and_flow(self) -> None: network, medium, high, low, valve = self._near_equal_pressure_network() solver = PressureFlowSolver(network, max_evaluations=10)