Merge model-development into main
This commit is contained in:
huojiarong committed
commit
127ec36a55
218 files changed
+65243
-47
No files matched your search
@@ -0,0 +1,237 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from math import acos, cos, isfinite, log, pi, sqrt
|
||||
|
||||
UNIVERSAL_GAS_CONSTANT = 8.31446261815324
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PengRobinsonFluid:
|
||||
"""Pure-fluid Peng-Robinson equation-of-state helper.
|
||||
|
||||
The class covers the equation-of-state layer plus the enthalpy departure
|
||||
needed to compare AMESim pneumatic ``pn2hpti`` reference enthalpy flows.
|
||||
"""
|
||||
|
||||
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 alpha_temperature_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 -(
|
||||
alpha_base
|
||||
* self.kappa
|
||||
/ (self.critical_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 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)
|
||||
|
||||
def residual_specific_enthalpy(
|
||||
self,
|
||||
pressure: float,
|
||||
temperature: float,
|
||||
phase: str = "vapor",
|
||||
) -> float:
|
||||
"""Return Peng-Robinson enthalpy departure from ideal gas, J/kg."""
|
||||
self._validate_pressure_temperature(pressure, temperature)
|
||||
z = self.compressibility_factor(pressure, temperature, phase=phase)
|
||||
_, B = self.reduced_parameters(pressure, temperature)
|
||||
b = self.b_parameter
|
||||
attractive = self.attractive_parameter(temperature)
|
||||
d_attractive_d_temperature = (
|
||||
self.attractive_parameter_temperature_derivative(temperature)
|
||||
)
|
||||
log_argument = (z + (1.0 + sqrt(2.0)) * B) / (
|
||||
z + (1.0 - sqrt(2.0)) * B
|
||||
)
|
||||
residual_molar_enthalpy = (
|
||||
UNIVERSAL_GAS_CONSTANT * temperature * (z - 1.0)
|
||||
+ (
|
||||
temperature * d_attractive_d_temperature
|
||||
- attractive
|
||||
)
|
||||
* log(log_argument)
|
||||
/ (2.0 * sqrt(2.0) * b)
|
||||
)
|
||||
return residual_molar_enthalpy / self.molar_mass
|
||||
|
||||
@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)
|
||||
Reference in new issue
Block a user