349 lines
12 KiB
Python
349 lines
12 KiB
Python
from __future__ import annotations
|
||
|
||
from collections.abc import Callable
|
||
from dataclasses import dataclass
|
||
from typing import ClassVar
|
||
|
||
from app.simulation.core.errors import RecoverableTrialStateError
|
||
from app.simulation.core.medium import (
|
||
GasMedium,
|
||
IdealGasMedium,
|
||
ThermodynamicProperties,
|
||
)
|
||
from app.simulation.core.peng_robinson import HELIUM_PR, PengRobinsonFluid
|
||
from app.simulation.performance import profile_property, record_property_iterations
|
||
from app.simulation.property_cache import cache_property_calculation
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class AmesimIdealAirMedium(IdealGasMedium):
|
||
"""AMESim air properties evaluated with the ideal-gas method.
|
||
|
||
Substance identity and property method are part of the concrete Python
|
||
type. A future air correlation or helium Peng-Robinson implementation can
|
||
therefore coexist as a sibling type without turning ``gi`` into a fluid
|
||
enumeration.
|
||
"""
|
||
|
||
SUBSTANCE_ID: ClassVar[str] = "air"
|
||
PROPERTY_METHOD_ID: ClassVar[str] = "ideal_gas"
|
||
|
||
name: str = "AMESimAirIdealGas"
|
||
R_gas: float = 287.0
|
||
cp_ref: float = 1005.0
|
||
T_ref: float = 300.0
|
||
cp_slope: float = 0.0
|
||
viscosity_ref: float = 1.82e-5
|
||
viscosity_T_ref: float = 293.15
|
||
sutherland_constant: float = 110.4
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
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 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 = 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
|
||
viscosity_T_ref: float = 293.15
|
||
sutherland_constant: float = 79.4
|
||
|
||
@property
|
||
def cv(self) -> float:
|
||
return (self.nasa_cp_over_R - 1.0) * self.R_gas
|
||
|
||
def cv_at_temperature(self, T: float) -> float:
|
||
del T
|
||
return self.cv
|
||
|
||
@profile_property("density")
|
||
@cache_property_calculation("density")
|
||
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
|
||
|
||
@profile_property("isentropic_density_pressure_factor")
|
||
@cache_property_calculation("isentropic_density_pressure_factor")
|
||
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)
|
||
|
||
@profile_property("specific_internal_energy")
|
||
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
|
||
)
|
||
|
||
@profile_property("specific_internal_energy_at_pressure")
|
||
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)
|
||
)
|
||
|
||
@profile_property("specific_enthalpy")
|
||
def specific_enthalpy(self, T: float) -> float:
|
||
return self.R_gas * (
|
||
self.nasa_cp_over_R * T
|
||
+ self.nasa_enthalpy_constant_K
|
||
)
|
||
|
||
@profile_property("specific_enthalpy_at_pressure")
|
||
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
|
||
|
||
@profile_property("temperature_from_pressure_enthalpy")
|
||
@cache_property_calculation("temperature_from_pressure_enthalpy")
|
||
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,
|
||
):
|
||
record_property_iterations(
|
||
"temperature_from_pressure_enthalpy",
|
||
_iteration + 1,
|
||
True,
|
||
)
|
||
return next_temperature
|
||
temperature = next_temperature
|
||
record_property_iterations(
|
||
"temperature_from_pressure_enthalpy",
|
||
16,
|
||
False,
|
||
)
|
||
return temperature
|
||
|
||
def temperature_from_mass_internal_energy(self, m: float, U: float) -> float:
|
||
if m <= 0.0:
|
||
raise RecoverableTrialStateError(
|
||
"Mass must stay positive when recovering temperature."
|
||
)
|
||
return self.temperature_from_internal_energy(U / m)
|
||
|
||
@profile_property("properties_from_mU")
|
||
@cache_property_calculation("properties_from_mU")
|
||
def properties_from_mU(
|
||
self,
|
||
m: float,
|
||
U: float,
|
||
V: float,
|
||
) -> ThermodynamicProperties:
|
||
"""Recover a real-gas state, reusing exact repeated evaluations.
|
||
|
||
Implicit integration asks several component interfaces for the same
|
||
``(m, U, V)`` state while closing one RHS evaluation and while building
|
||
finite-difference Jacobians. The calculation is pure and its result is
|
||
immutable, so an exact-key bounded cache avoids repeating the
|
||
Peng-Robinson temperature iteration without changing model semantics.
|
||
"""
|
||
if m <= 0.0:
|
||
raise RecoverableTrialStateError(
|
||
"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,
|
||
)
|
||
converged = False
|
||
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
|
||
converged = True
|
||
break
|
||
temperature = next_temperature
|
||
record_property_iterations(
|
||
"properties_from_mU",
|
||
_iteration + 1,
|
||
converged,
|
||
)
|
||
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:
|
||
"""A selectable calculation method for one AMESim gas substance."""
|
||
|
||
value: int
|
||
label: str
|
||
method_id: str
|
||
factory: Callable[[], GasMedium]
|
||
eos_type: int
|
||
|
||
def build_medium(self) -> GasMedium:
|
||
return self.factory()
|
||
|
||
|
||
AMESIM_AIR_IDEAL_GAS_PROPERTY_MODEL = 0
|
||
AMESIM_AIR_PROPERTY_MODELS = (
|
||
AmesimGasPropertyModelSpec(
|
||
value=AMESIM_AIR_IDEAL_GAS_PROPERTY_MODEL,
|
||
label="理想气体",
|
||
method_id=AmesimIdealAirMedium.PROPERTY_METHOD_ID,
|
||
factory=AmesimIdealAirMedium,
|
||
eos_type=1,
|
||
),
|
||
)
|
||
|
||
|
||
AMESIM_HELIUM_PENG_ROBINSON_PROPERTY_MODEL = 0
|
||
AMESIM_HELIUM_PROPERTY_MODELS = (
|
||
AmesimGasPropertyModelSpec(
|
||
value=AMESIM_HELIUM_PENG_ROBINSON_PROPERTY_MODEL,
|
||
label="Peng–Robinson",
|
||
method_id=AmesimHeliumPengRobinsonMedium.PROPERTY_METHOD_ID,
|
||
factory=AmesimHeliumPengRobinsonMedium,
|
||
eos_type=6,
|
||
),
|
||
)
|