feat: integrate AMESim media models and editor UI
This commit is contained in:
ljz committed
1 parent
3e4f466114
commit
e7177ab03e
181 files changed
+7042
-25270
No files matched your search
@@ -0,0 +1 @@
|
||||
"""Calibrated component primitives used only by the ``test_mql`` example."""
|
||||
@@ -0,0 +1,168 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from math import pi
|
||||
|
||||
|
||||
MM_TO_M = 1.0e-3
|
||||
M_TO_MM = 1.0e3
|
||||
M3_TO_CM3 = 1.0e6
|
||||
M3_PER_S_TO_L_PER_MIN = 60_000.0
|
||||
|
||||
|
||||
def circular_area(diameter_m: float) -> float:
|
||||
if diameter_m < 0.0:
|
||||
raise ValueError("diameter_m must be non-negative.")
|
||||
return pi * diameter_m * diameter_m / 4.0
|
||||
|
||||
|
||||
def mm_to_m(value: float) -> float:
|
||||
return value * MM_TO_M
|
||||
|
||||
|
||||
def m_to_mm(value: float) -> float:
|
||||
return value * M_TO_MM
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AmesimPistonGeometry:
|
||||
"""Geometry relations used by AMESim PNRP17 pneumatic piston variables."""
|
||||
|
||||
piston_diameter_m: float
|
||||
rod_diameter_m: float = 0.0
|
||||
zero_length_m: float = 0.0
|
||||
|
||||
@property
|
||||
def piston_area_m2(self) -> float:
|
||||
return circular_area(self.piston_diameter_m)
|
||||
|
||||
@property
|
||||
def rod_area_m2(self) -> float:
|
||||
return circular_area(self.rod_diameter_m)
|
||||
|
||||
@property
|
||||
def annulus_area_m2(self) -> float:
|
||||
return self.piston_area_m2 - self.rod_area_m2
|
||||
|
||||
def chamber_length_m(self, port4_displacement_m: float, port5_displacement_m: float) -> float:
|
||||
return self.zero_length_m + port5_displacement_m - port4_displacement_m
|
||||
|
||||
def chamber_length_mm(self, port4_displacement_m: float, port5_displacement_m: float) -> float:
|
||||
return m_to_mm(self.chamber_length_m(port4_displacement_m, port5_displacement_m))
|
||||
|
||||
@property
|
||||
def chamber_area_m2(self) -> float:
|
||||
return self.annulus_area_m2
|
||||
|
||||
def chamber_volume_m3(self, port4_displacement_m: float, port5_displacement_m: float) -> float:
|
||||
return self.chamber_area_m2 * self.chamber_length_m(
|
||||
port4_displacement_m,
|
||||
port5_displacement_m,
|
||||
)
|
||||
|
||||
def chamber_volume_cm3(self, port4_displacement_m: float, port5_displacement_m: float) -> float:
|
||||
return self.chamber_volume_m3(port4_displacement_m, port5_displacement_m) * M3_TO_CM3
|
||||
|
||||
def chamber_volume_rate_m3_s(self, port4_velocity_m_s: float, port5_velocity_m_s: float) -> float:
|
||||
return self.chamber_area_m2 * (port5_velocity_m_s - port4_velocity_m_s)
|
||||
|
||||
def chamber_volume_rate_l_min(self, port4_velocity_m_s: float, port5_velocity_m_s: float) -> float:
|
||||
return self.chamber_volume_rate_m3_s(
|
||||
port4_velocity_m_s,
|
||||
port5_velocity_m_s,
|
||||
) * M3_PER_S_TO_L_PER_MIN
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AmesimElasticEndstop:
|
||||
"""Contact force part of AMESim LSTP00A elastic endstop."""
|
||||
|
||||
contact_stiffness_n_per_m: float
|
||||
contact_damping_n_per_m_per_s: float = 0.0
|
||||
gap0_m: float = 0.0
|
||||
|
||||
def penetration_m_from_gap_mm(self, gap_mm: float) -> float:
|
||||
return max(-(mm_to_m(gap_mm) - self.gap0_m), 0.0)
|
||||
|
||||
def static_contact_force(self, gap_mm: float) -> float:
|
||||
return self.contact_stiffness_n_per_m * self.penetration_m_from_gap_mm(gap_mm)
|
||||
|
||||
def contact_force(self, gap_mm: float, penetration_velocity_m_s: float = 0.0) -> float:
|
||||
if self.penetration_m_from_gap_mm(gap_mm) <= 0.0:
|
||||
return 0.0
|
||||
damping_force = self.contact_damping_n_per_m_per_s * penetration_velocity_m_s
|
||||
return max(self.static_contact_force(gap_mm) + damping_force, 0.0)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AmesimMassFrictionEndstops:
|
||||
"""Parameter and observable helpers for AMESim MECMAS21 translation masses."""
|
||||
|
||||
mass_kg: float
|
||||
lower_limit_m: float
|
||||
upper_limit_m: float
|
||||
lower_stiffness_n_per_m: float
|
||||
upper_stiffness_n_per_m: float
|
||||
lower_damping_n_per_m_per_s: float = 0.0
|
||||
upper_damping_n_per_m_per_s: float = 0.0
|
||||
viscous_friction_n_per_m_per_s: float = 0.0
|
||||
coulomb_friction_n: float = 0.0
|
||||
stiction_force_n: float = 0.0
|
||||
windage_n_per_m2_per_s2: float = 0.0
|
||||
|
||||
def lower_penetration_m(self, displacement_m: float) -> float:
|
||||
return max(self.lower_limit_m - displacement_m, 0.0)
|
||||
|
||||
def upper_penetration_m(self, displacement_m: float) -> float:
|
||||
return max(displacement_m - self.upper_limit_m, 0.0)
|
||||
|
||||
def lower_static_force_magnitude(self, displacement_m: float) -> float:
|
||||
return self.lower_stiffness_n_per_m * self.lower_penetration_m(displacement_m)
|
||||
|
||||
def upper_static_force_magnitude(self, displacement_m: float) -> float:
|
||||
return self.upper_stiffness_n_per_m * self.upper_penetration_m(displacement_m)
|
||||
|
||||
def viscous_friction_force(self, velocity_m_s: float) -> float:
|
||||
return -self.viscous_friction_n_per_m_per_s * velocity_m_s
|
||||
|
||||
def windage_force(self, velocity_m_s: float) -> float:
|
||||
return -self.windage_n_per_m2_per_s2 * velocity_m_s * abs(velocity_m_s)
|
||||
|
||||
def dry_friction_force(self, velocity_m_s: float) -> float:
|
||||
if velocity_m_s > 0.0:
|
||||
return -self.coulomb_friction_n
|
||||
if velocity_m_s < 0.0:
|
||||
return self.coulomb_friction_n
|
||||
return 0.0
|
||||
|
||||
def limit_contact_force(self, displacement_m: float, velocity_m_s: float) -> float:
|
||||
lower_force = self.lower_static_force_magnitude(displacement_m)
|
||||
if lower_force > 0.0:
|
||||
lower_force += max(-self.lower_damping_n_per_m_per_s * velocity_m_s, 0.0)
|
||||
|
||||
upper_force = self.upper_static_force_magnitude(displacement_m)
|
||||
if upper_force > 0.0:
|
||||
upper_force += max(self.upper_damping_n_per_m_per_s * velocity_m_s, 0.0)
|
||||
|
||||
return lower_force - upper_force
|
||||
|
||||
def derivatives(
|
||||
self,
|
||||
*,
|
||||
velocity_m_s: float,
|
||||
displacement_m: float,
|
||||
port_1_force_n: float = 0.0,
|
||||
port_2_force_n: float = 0.0,
|
||||
external_force_n: float = 0.0,
|
||||
) -> tuple[float, float]:
|
||||
total_force = (
|
||||
port_1_force_n
|
||||
+ port_2_force_n
|
||||
+ external_force_n
|
||||
+ self.viscous_friction_force(velocity_m_s)
|
||||
+ self.windage_force(velocity_m_s)
|
||||
+ self.dry_friction_force(velocity_m_s)
|
||||
+ self.limit_contact_force(displacement_m, velocity_m_s)
|
||||
)
|
||||
return total_force / self.mass_kg, velocity_m_s
|
||||
|
||||
@@ -0,0 +1,437 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from math import pi, sqrt
|
||||
|
||||
from app.simulation.core.base import AlgebraicComponent, DynamicComponent
|
||||
from app.simulation.core.medium import ThermodynamicProperties
|
||||
from app.simulation.core.peng_robinson import HELIUM_PR, PengRobinsonFluid
|
||||
from app.simulation.core.ports import PortState
|
||||
from app.simulation.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 specific_reference_enthalpy(
|
||||
self,
|
||||
temperature: float,
|
||||
reference_temperature: float = 298.15,
|
||||
) -> float:
|
||||
return self.cp * (temperature - reference_temperature)
|
||||
|
||||
def reference_temperature_from_specific_enthalpy(
|
||||
self,
|
||||
specific_enthalpy: float,
|
||||
reference_temperature: float = 298.15,
|
||||
) -> float:
|
||||
if self.cp <= 0.0:
|
||||
raise ValueError("cp must be positive.")
|
||||
return reference_temperature + specific_enthalpy / self.cp
|
||||
|
||||
def pressure_reference_enthalpy(
|
||||
self,
|
||||
pressure: float,
|
||||
temperature: float,
|
||||
reference_pressure: float = 101_300.0,
|
||||
reference_temperature: float = 298.15,
|
||||
) -> float:
|
||||
return (
|
||||
self.specific_reference_enthalpy(temperature, reference_temperature)
|
||||
+ self.fluid.residual_specific_enthalpy(pressure, temperature)
|
||||
- self.fluid.residual_specific_enthalpy(
|
||||
reference_pressure,
|
||||
reference_temperature,
|
||||
)
|
||||
)
|
||||
|
||||
def pressure_transport_enthalpy(
|
||||
self,
|
||||
pressure: float,
|
||||
temperature: float,
|
||||
reference_pressure: float = 101_300.0,
|
||||
reference_temperature: float = 298.15,
|
||||
) -> float:
|
||||
"""Convert AMESim reference enthalpy to the absolute-energy state basis."""
|
||||
return (
|
||||
self.pressure_reference_enthalpy(
|
||||
pressure,
|
||||
temperature,
|
||||
reference_pressure,
|
||||
reference_temperature,
|
||||
)
|
||||
+ self.cp * reference_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 m3_to_cm3(value: float) -> float:
|
||||
return value * 1.0e6
|
||||
|
||||
|
||||
def cm3_to_m3(value: float) -> float:
|
||||
return value * 1.0e-6
|
||||
|
||||
|
||||
def kg_to_g(value: float) -> float:
|
||||
return value * 1.0e3
|
||||
|
||||
|
||||
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,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
heat_transfer_area: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
) -> None:
|
||||
if volume <= 0.0:
|
||||
raise ValueError("volume must be positive.")
|
||||
if heat_transfer_coefficient < 0.0:
|
||||
raise ValueError("heat_transfer_coefficient must be non-negative.")
|
||||
if heat_transfer_area < 0.0:
|
||||
raise ValueError("heat_transfer_area must be non-negative.")
|
||||
if external_temperature_k <= 0.0:
|
||||
raise ValueError("external_temperature_k must be positive.")
|
||||
super().__init__(name=name)
|
||||
self.volume = volume
|
||||
self.gas = gas
|
||||
self.heat_transfer_coefficient = heat_transfer_coefficient
|
||||
self.heat_transfer_area = heat_transfer_area
|
||||
self.external_temperature = external_temperature_k
|
||||
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()
|
||||
self.port_b = 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,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
heat_transfer_area: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
) -> "AmesimPneumaticVolume":
|
||||
return cls(
|
||||
name=name,
|
||||
volume=liters_to_m3(volume_liters),
|
||||
gas=gas,
|
||||
p0=p0,
|
||||
T0=T0,
|
||||
heat_transfer_coefficient=heat_transfer_coefficient,
|
||||
heat_transfer_area=heat_transfer_area,
|
||||
external_temperature_k=external_temperature_k,
|
||||
)
|
||||
|
||||
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 volume_cm3(self) -> float:
|
||||
return m3_to_cm3(self.volume)
|
||||
|
||||
def volume_rate_m3_s(self) -> float:
|
||||
return 0.0
|
||||
|
||||
def thermal_energy_flow_w(self, temperature_k: float | None = None) -> float:
|
||||
temperature = self.properties().T if temperature_k is None else temperature_k
|
||||
return (
|
||||
self.heat_transfer_coefficient
|
||||
* self.heat_transfer_area
|
||||
* (self.external_temperature - temperature)
|
||||
)
|
||||
|
||||
def gas_mass_g(self) -> float:
|
||||
return kg_to_g(self.state.m)
|
||||
|
||||
def pressure_gauge_pa(self, reference_pressure_pa: float = 101_300.0) -> float:
|
||||
return self.properties().p - reference_pressure_pa
|
||||
|
||||
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
|
||||
self.port_b.p = p
|
||||
self.port_b.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 + self.thermal_energy_flow_w(),
|
||||
)
|
||||
|
||||
def derivatives_from_two_connections(
|
||||
self,
|
||||
*,
|
||||
port_a_m_flow: float,
|
||||
connected_h_a: float,
|
||||
port_b_m_flow: float,
|
||||
connected_h_b: float,
|
||||
internal_h: float,
|
||||
volume_rate_m3_s: float | None = None,
|
||||
) -> VolumeState:
|
||||
properties = self.properties()
|
||||
inlet_h_a = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_a_m_flow,
|
||||
connected_h=connected_h_a,
|
||||
internal_h=internal_h,
|
||||
)
|
||||
inlet_h_b = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_b_m_flow,
|
||||
connected_h=connected_h_b,
|
||||
internal_h=internal_h,
|
||||
)
|
||||
return VolumeState(
|
||||
m=port_a_m_flow + port_b_m_flow,
|
||||
U=(
|
||||
port_a_m_flow * inlet_h_a
|
||||
+ port_b_m_flow * inlet_h_b
|
||||
+ self.thermal_energy_flow_w(properties.T)
|
||||
- properties.p * (
|
||||
self.volume_rate_m3_s()
|
||||
if volume_rate_m3_s is None
|
||||
else volume_rate_m3_s
|
||||
)
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
class AmesimVariablePneumaticVolume(AmesimPneumaticVolume):
|
||||
"""PNCH012-style volume with a dead volume plus an external moving volume."""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
name: str,
|
||||
dead_volume: float,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
p0: float = 101_325.0,
|
||||
T0: float = 293.15,
|
||||
external_volume: float = 0.0,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
heat_transfer_area: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
) -> None:
|
||||
if dead_volume <= 0.0:
|
||||
raise ValueError("dead_volume must be positive.")
|
||||
if dead_volume + external_volume <= 0.0:
|
||||
raise ValueError("total volume must be positive.")
|
||||
self.dead_volume = dead_volume
|
||||
self.external_volume = external_volume
|
||||
self.external_volume_rate = 0.0
|
||||
super().__init__(
|
||||
name=name,
|
||||
volume=dead_volume + external_volume,
|
||||
gas=gas,
|
||||
p0=p0,
|
||||
T0=T0,
|
||||
heat_transfer_coefficient=heat_transfer_coefficient,
|
||||
heat_transfer_area=heat_transfer_area,
|
||||
external_temperature_k=external_temperature_k,
|
||||
)
|
||||
|
||||
@classmethod
|
||||
def from_liters(
|
||||
cls,
|
||||
name: str,
|
||||
dead_volume_liters: float,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
p0: float = 101_325.0,
|
||||
T0: float = 293.15,
|
||||
external_volume_liters: float = 0.0,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
heat_transfer_area: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
) -> "AmesimVariablePneumaticVolume":
|
||||
return cls(
|
||||
name=name,
|
||||
dead_volume=liters_to_m3(dead_volume_liters),
|
||||
gas=gas,
|
||||
p0=p0,
|
||||
T0=T0,
|
||||
external_volume=liters_to_m3(external_volume_liters),
|
||||
heat_transfer_coefficient=heat_transfer_coefficient,
|
||||
heat_transfer_area=heat_transfer_area,
|
||||
external_temperature_k=external_temperature_k,
|
||||
)
|
||||
|
||||
def volume_rate_m3_s(self) -> float:
|
||||
return self.external_volume_rate
|
||||
|
||||
def set_external_volume_m3(
|
||||
self,
|
||||
external_volume: float,
|
||||
external_volume_rate_m3_s: float = 0.0,
|
||||
) -> None:
|
||||
if self.dead_volume + external_volume <= 0.0:
|
||||
raise ValueError("total volume must be positive.")
|
||||
self.external_volume = external_volume
|
||||
self.external_volume_rate = external_volume_rate_m3_s
|
||||
self.volume = self.dead_volume + self.external_volume
|
||||
|
||||
|
||||
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:
|
||||
opening = min(max(self.opening, 0.0), 1.0)
|
||||
return self.area * opening
|
||||
|
||||
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
|
||||
@@ -0,0 +1,881 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from math import log10, pi, sqrt
|
||||
|
||||
from app.simulation.examples.test_mql.primitives.pneumatic import (
|
||||
HELIUM_PNEUMATIC_GAS,
|
||||
AmesimPneumaticGas,
|
||||
compressible_orifice_mass_flow,
|
||||
diameter_mm_to_area_m2,
|
||||
)
|
||||
from app.simulation.core.base import AlgebraicComponent, DynamicComponent
|
||||
from app.simulation.core.medium import ThermodynamicProperties
|
||||
from app.simulation.core.ports import PortState
|
||||
from app.simulation.core.state import VolumeState
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AmesimPnl0001Diagnostics:
|
||||
mass_flow_kg_s: float
|
||||
reynolds_number: float
|
||||
gas_velocity_m_s: float
|
||||
friction_factor: float
|
||||
pressure_drop_pa: float
|
||||
|
||||
|
||||
class _DarcyPipeResistanceMixin:
|
||||
diameter: float
|
||||
length: float
|
||||
relative_roughness: float
|
||||
area: float
|
||||
|
||||
def _mass_flow_for_pressure_drop(
|
||||
self,
|
||||
pressure_drop_pa: float,
|
||||
*,
|
||||
density: float,
|
||||
temperature: float,
|
||||
) -> float:
|
||||
if pressure_drop_pa <= 0.0:
|
||||
return 0.0
|
||||
upper = 1.0e-9
|
||||
while self._darcy_pressure_drop(
|
||||
upper,
|
||||
density=density,
|
||||
temperature=temperature,
|
||||
) < pressure_drop_pa:
|
||||
upper *= 10.0
|
||||
if upper > 1.0e3:
|
||||
raise ValueError("unable to bracket pneumatic pipe resistance flow")
|
||||
lower = 0.0
|
||||
for _ in range(48):
|
||||
middle = 0.5 * (lower + upper)
|
||||
if self._darcy_pressure_drop(
|
||||
middle,
|
||||
density=density,
|
||||
temperature=temperature,
|
||||
) < pressure_drop_pa:
|
||||
lower = middle
|
||||
else:
|
||||
upper = middle
|
||||
return 0.5 * (lower + upper)
|
||||
|
||||
def pn2pipefr_mass_flow(
|
||||
self,
|
||||
*,
|
||||
port_1_pressure_pa: float,
|
||||
port_1_temperature_k: float,
|
||||
port_2_pressure_pa: float,
|
||||
port_2_temperature_k: float,
|
||||
length: float | None = None,
|
||||
) -> float:
|
||||
pressure_difference = port_1_pressure_pa - port_2_pressure_pa
|
||||
if pressure_difference == 0.0:
|
||||
return 0.0
|
||||
upstream_pressure = max(port_1_pressure_pa, port_2_pressure_pa)
|
||||
downstream_pressure = min(port_1_pressure_pa, port_2_pressure_pa)
|
||||
upstream_temperature = (
|
||||
port_1_temperature_k
|
||||
if pressure_difference > 0.0
|
||||
else port_2_temperature_k
|
||||
)
|
||||
resistance_length = self.length if length is None else length
|
||||
if resistance_length <= 0.0:
|
||||
raise ValueError("length must be positive")
|
||||
|
||||
def target_flow(mass_flow_kg_s: float) -> float:
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, upstream_temperature)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
flow_coefficient = sqrt(
|
||||
self.diameter / (resistance_length * friction_factor)
|
||||
)
|
||||
return compressible_orifice_mass_flow(
|
||||
upstream_pressure=upstream_pressure,
|
||||
downstream_pressure=downstream_pressure,
|
||||
upstream_temperature=upstream_temperature,
|
||||
area=self.area,
|
||||
flow_coefficient=flow_coefficient,
|
||||
gas=self.gas,
|
||||
)
|
||||
|
||||
flow_coefficient = sqrt(self.diameter / (resistance_length * 0.02))
|
||||
magnitude = compressible_orifice_mass_flow(
|
||||
upstream_pressure=upstream_pressure,
|
||||
downstream_pressure=downstream_pressure,
|
||||
upstream_temperature=upstream_temperature,
|
||||
area=self.area,
|
||||
flow_coefficient=flow_coefficient,
|
||||
gas=self.gas,
|
||||
)
|
||||
for _ in range(12):
|
||||
next_magnitude = target_flow(magnitude)
|
||||
if abs(next_magnitude - magnitude) <= max(1.0e-12, abs(magnitude) * 1.0e-9):
|
||||
magnitude = next_magnitude
|
||||
break
|
||||
magnitude = 0.5 * (magnitude + next_magnitude)
|
||||
return magnitude if pressure_difference > 0.0 else -magnitude
|
||||
|
||||
def _darcy_pressure_drop(
|
||||
self,
|
||||
mass_flow_kg_s: float,
|
||||
*,
|
||||
density: float,
|
||||
temperature: float,
|
||||
) -> float:
|
||||
if mass_flow_kg_s == 0.0:
|
||||
return 0.0
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
velocity = mass_flow_kg_s / (density * self.area)
|
||||
magnitude = (
|
||||
friction_factor
|
||||
* (self.length / self.diameter)
|
||||
* density
|
||||
* velocity
|
||||
* velocity
|
||||
/ 2.0
|
||||
)
|
||||
return magnitude if mass_flow_kg_s > 0.0 else -magnitude
|
||||
|
||||
def _reynolds_number(self, mass_flow_kg_s: float, temperature: float) -> float:
|
||||
viscosity = helium_dynamic_viscosity(temperature)
|
||||
return 4.0 * abs(mass_flow_kg_s) / (pi * self.diameter * viscosity)
|
||||
|
||||
def _friction_factor(self, reynolds_number: float) -> float:
|
||||
if reynolds_number <= 0.0:
|
||||
return 64_000_000.0
|
||||
laminar = 64.0 / reynolds_number
|
||||
if reynolds_number <= 2_300.0:
|
||||
return laminar
|
||||
turbulent = 1.0 / (
|
||||
-1.8
|
||||
* log10(
|
||||
(self.relative_roughness / 3.7) ** 1.11
|
||||
+ 6.9 / reynolds_number
|
||||
)
|
||||
) ** 2
|
||||
if reynolds_number >= 4_000.0:
|
||||
return turbulent
|
||||
fraction = (reynolds_number - 2_300.0) / 1_700.0
|
||||
return laminar + fraction * (turbulent - laminar)
|
||||
|
||||
|
||||
class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
|
||||
"""Physical first-pass implementation of AMESim ``PNL0001`` (C-R).
|
||||
|
||||
Port 2 owns the lumped gas storage. Port 1 is connected through a Darcy
|
||||
resistance. Both connection mass flows use the simulation convention:
|
||||
positive values enter the pipe storage.
|
||||
|
||||
AMESim's proprietary ``pn2pipefr`` utility is represented by an
|
||||
optional calibrated linear conductance when a model-specific baseline
|
||||
supports it; otherwise the component falls back to an auditable
|
||||
Darcy-Weisbach law. Both paths preserve the real geometry, state count,
|
||||
mass/energy balance, heat-transfer parameter, and observable diagnostics.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
name: str,
|
||||
*,
|
||||
diameter_mm: float,
|
||||
length_m: float,
|
||||
relative_roughness: float,
|
||||
polytropic_constant: float = 1.35,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
calibrated_linear_conductance: float | None = None,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
p0: float = 101_325.0,
|
||||
T0: float = 293.15,
|
||||
) -> None:
|
||||
if diameter_mm <= 0.0:
|
||||
raise ValueError("diameter_mm must be positive")
|
||||
if length_m <= 0.0:
|
||||
raise ValueError("length_m must be positive")
|
||||
if relative_roughness < 0.0:
|
||||
raise ValueError("relative_roughness must be non-negative")
|
||||
if polytropic_constant <= 0.0:
|
||||
raise ValueError("polytropic_constant must be positive")
|
||||
if heat_transfer_coefficient < 0.0:
|
||||
raise ValueError("heat_transfer_coefficient must be non-negative")
|
||||
if external_temperature_k <= 0.0:
|
||||
raise ValueError("external_temperature_k must be positive")
|
||||
if (
|
||||
calibrated_linear_conductance is not None
|
||||
and calibrated_linear_conductance <= 0.0
|
||||
):
|
||||
raise ValueError("calibrated_linear_conductance must be positive")
|
||||
|
||||
super().__init__(name=name)
|
||||
self.diameter = diameter_mm * 1.0e-3
|
||||
self.length = length_m
|
||||
self.relative_roughness = relative_roughness
|
||||
self.polytropic_constant = polytropic_constant
|
||||
self.heat_transfer_coefficient = heat_transfer_coefficient
|
||||
self.external_temperature = external_temperature_k
|
||||
self.calibrated_linear_conductance = calibrated_linear_conductance
|
||||
self.gas = gas
|
||||
self.area = diameter_mm_to_area_m2(diameter_mm)
|
||||
self.volume = self.area * self.length
|
||||
self.heat_transfer_area = pi * self.diameter * self.length
|
||||
|
||||
rho0 = gas.density(p0, T0)
|
||||
mass0 = rho0 * self.volume
|
||||
self.state = VolumeState(
|
||||
m=mass0,
|
||||
U=mass0 * gas.specific_internal_energy(T0),
|
||||
)
|
||||
self.port_1 = PortState()
|
||||
self.port_2 = PortState()
|
||||
|
||||
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("pipe mass must stay positive")
|
||||
temperature = self.gas.temperature_from_internal_energy(
|
||||
self.state.U / self.state.m
|
||||
)
|
||||
density = self.state.m / self.volume
|
||||
pressure = self.gas.pressure(density, temperature)
|
||||
properties = ThermodynamicProperties(
|
||||
p=pressure,
|
||||
T=temperature,
|
||||
rho=density,
|
||||
u=self.state.U / self.state.m,
|
||||
h=self.gas.specific_enthalpy(temperature),
|
||||
)
|
||||
self.port_2.p = pressure
|
||||
self.port_2.h_outflow = properties.h
|
||||
return properties
|
||||
|
||||
def gas_mass_g(self) -> float:
|
||||
return self.state.m * 1.0e3
|
||||
|
||||
def resistance_mass_flow(
|
||||
self,
|
||||
*,
|
||||
port_1_pressure_pa: float,
|
||||
port_1_temperature_k: float,
|
||||
) -> float:
|
||||
"""Return mass flow from port 1 into the port-2 storage in kg/s."""
|
||||
if port_1_pressure_pa <= 0.0:
|
||||
raise ValueError("port_1_pressure_pa must be positive")
|
||||
if port_1_temperature_k <= 0.0:
|
||||
raise ValueError("port_1_temperature_k must be positive")
|
||||
|
||||
internal = self.properties()
|
||||
pressure_difference = port_1_pressure_pa - internal.p
|
||||
if pressure_difference == 0.0:
|
||||
return 0.0
|
||||
if self.calibrated_linear_conductance is not None:
|
||||
return (
|
||||
self.calibrated_linear_conductance
|
||||
* pressure_difference
|
||||
/ sqrt(internal.T)
|
||||
)
|
||||
upstream_pressure = max(port_1_pressure_pa, internal.p)
|
||||
upstream_temperature = (
|
||||
port_1_temperature_k if pressure_difference > 0.0 else internal.T
|
||||
)
|
||||
density = self.gas.density(upstream_pressure, upstream_temperature)
|
||||
magnitude = self._mass_flow_for_pressure_drop(
|
||||
abs(pressure_difference),
|
||||
density=density,
|
||||
temperature=upstream_temperature,
|
||||
)
|
||||
return magnitude if pressure_difference > 0.0 else -magnitude
|
||||
|
||||
def diagnostics(
|
||||
self,
|
||||
*,
|
||||
mass_flow_kg_s: float,
|
||||
temperature_k: float | None = None,
|
||||
) -> AmesimPnl0001Diagnostics:
|
||||
properties = self.properties()
|
||||
temperature = temperature_k or properties.T
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
velocity = mass_flow_kg_s / (properties.rho * self.area)
|
||||
pressure_drop = self._darcy_pressure_drop(
|
||||
mass_flow_kg_s,
|
||||
density=properties.rho,
|
||||
temperature=temperature,
|
||||
)
|
||||
return AmesimPnl0001Diagnostics(
|
||||
mass_flow_kg_s=mass_flow_kg_s,
|
||||
reynolds_number=reynolds,
|
||||
gas_velocity_m_s=velocity,
|
||||
friction_factor=friction_factor,
|
||||
pressure_drop_pa=pressure_drop,
|
||||
)
|
||||
|
||||
def darcy_pressure_drop_for_state(
|
||||
self,
|
||||
*,
|
||||
mass_flow_kg_s: float,
|
||||
pressure_pa: float,
|
||||
temperature_k: float,
|
||||
) -> float:
|
||||
if pressure_pa <= 0.0:
|
||||
raise ValueError("pressure_pa must be positive")
|
||||
if temperature_k <= 0.0:
|
||||
raise ValueError("temperature_k must be positive")
|
||||
density = self.gas.density(pressure_pa, temperature_k)
|
||||
return self._darcy_pressure_drop(
|
||||
mass_flow_kg_s,
|
||||
density=density,
|
||||
temperature=temperature_k,
|
||||
)
|
||||
|
||||
def derivatives_from_connections(
|
||||
self,
|
||||
*,
|
||||
port_1_m_flow: float,
|
||||
connected_h_1: float,
|
||||
port_2_m_flow: float,
|
||||
connected_h_2: float,
|
||||
) -> VolumeState:
|
||||
internal = self.properties()
|
||||
# Default first-pass PNL0001 behavior uses the historical internal-energy
|
||||
# approximation. AMESim-specific transport-enthalpy corrections are kept
|
||||
# behind derivatives_from_transport_enthalpy_connections so they can be
|
||||
# applied only where validated against baseline data.
|
||||
inlet_u_1 = (
|
||||
connected_h_1 / self.gas.gamma
|
||||
if port_1_m_flow > 0.0
|
||||
else internal.u
|
||||
)
|
||||
inlet_u_2 = (
|
||||
connected_h_2 / self.gas.gamma
|
||||
if port_2_m_flow > 0.0
|
||||
else internal.u
|
||||
)
|
||||
heat_flow = (
|
||||
self.heat_transfer_coefficient
|
||||
* self.heat_transfer_area
|
||||
* (self.external_temperature - internal.T)
|
||||
)
|
||||
return VolumeState(
|
||||
m=port_1_m_flow + port_2_m_flow,
|
||||
U=port_1_m_flow * inlet_u_1 + port_2_m_flow * inlet_u_2 + heat_flow,
|
||||
)
|
||||
|
||||
def derivatives_from_transport_enthalpy_connections(
|
||||
self,
|
||||
*,
|
||||
port_1_m_flow: float,
|
||||
connected_h_1: float,
|
||||
port_2_m_flow: float,
|
||||
connected_h_2: float,
|
||||
) -> VolumeState:
|
||||
internal = self.properties()
|
||||
inlet_h_1 = connected_h_1 if port_1_m_flow > 0.0 else internal.h
|
||||
inlet_h_2 = connected_h_2 if port_2_m_flow > 0.0 else internal.h
|
||||
heat_flow = (
|
||||
self.heat_transfer_coefficient
|
||||
* self.heat_transfer_area
|
||||
* (self.external_temperature - internal.T)
|
||||
)
|
||||
return VolumeState(
|
||||
m=port_1_m_flow + port_2_m_flow,
|
||||
U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow,
|
||||
)
|
||||
|
||||
|
||||
class AmesimPnl0003Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
|
||||
"""First-pass AMESim ``PNL0003`` (C-R-C) pipe.
|
||||
|
||||
The two pipe-end compliances are represented as equal half-volume gas
|
||||
stores connected by the same auditable Darcy resistance used for PNL0001.
|
||||
Center flow is positive from port 1 storage to port 2 storage.
|
||||
"""
|
||||
|
||||
state_size = 4
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
name: str,
|
||||
*,
|
||||
diameter_mm: float,
|
||||
length_m: float,
|
||||
relative_roughness: float,
|
||||
polytropic_constant: float = 1.35,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
p1_0: float = 101_325.0,
|
||||
T1_0: float = 293.15,
|
||||
p2_0: float = 101_325.0,
|
||||
T2_0: float = 293.15,
|
||||
) -> None:
|
||||
if diameter_mm <= 0.0:
|
||||
raise ValueError("diameter_mm must be positive")
|
||||
if length_m <= 0.0:
|
||||
raise ValueError("length_m must be positive")
|
||||
if relative_roughness < 0.0:
|
||||
raise ValueError("relative_roughness must be non-negative")
|
||||
if polytropic_constant <= 0.0:
|
||||
raise ValueError("polytropic_constant must be positive")
|
||||
if heat_transfer_coefficient < 0.0:
|
||||
raise ValueError("heat_transfer_coefficient must be non-negative")
|
||||
if external_temperature_k <= 0.0:
|
||||
raise ValueError("external_temperature_k must be positive")
|
||||
|
||||
super().__init__(name=name)
|
||||
self.diameter = diameter_mm * 1.0e-3
|
||||
self.length = length_m
|
||||
self.relative_roughness = relative_roughness
|
||||
self.polytropic_constant = polytropic_constant
|
||||
self.heat_transfer_coefficient = heat_transfer_coefficient
|
||||
self.external_temperature = external_temperature_k
|
||||
self.gas = gas
|
||||
self.area = diameter_mm_to_area_m2(diameter_mm)
|
||||
self.volume = self.area * self.length
|
||||
self.compliance_volume = self.volume / 2.0
|
||||
self.heat_transfer_area = pi * self.diameter * self.length
|
||||
|
||||
self.state_1 = self._initial_state(p1_0, T1_0)
|
||||
self.state_2 = self._initial_state(p2_0, T2_0)
|
||||
self.port_1 = PortState()
|
||||
self.port_2 = PortState()
|
||||
|
||||
def _initial_state(self, pressure: float, temperature: float) -> VolumeState:
|
||||
rho = self.gas.density(pressure, temperature)
|
||||
mass = rho * self.compliance_volume
|
||||
return VolumeState(
|
||||
m=mass,
|
||||
U=mass * self.gas.specific_internal_energy(temperature),
|
||||
)
|
||||
|
||||
def get_state_vector(self) -> list[float]:
|
||||
return [*self.state_1.as_vector(), *self.state_2.as_vector()]
|
||||
|
||||
def set_state_vector(self, values: list[float]) -> None:
|
||||
if len(values) != 4:
|
||||
raise ValueError("PNL0003 state vector requires four values")
|
||||
self.state_1 = VolumeState.from_vector(values[:2])
|
||||
self.state_2 = VolumeState.from_vector(values[2:])
|
||||
|
||||
def properties_1(self) -> ThermodynamicProperties:
|
||||
properties = self._properties(self.state_1)
|
||||
self.port_1.p = properties.p
|
||||
self.port_1.h_outflow = properties.h
|
||||
return properties
|
||||
|
||||
def properties_2(self) -> ThermodynamicProperties:
|
||||
properties = self._properties(self.state_2)
|
||||
self.port_2.p = properties.p
|
||||
self.port_2.h_outflow = properties.h
|
||||
return properties
|
||||
|
||||
def _properties(self, state: VolumeState) -> ThermodynamicProperties:
|
||||
if state.m <= 0.0:
|
||||
raise ValueError("pipe mass must stay positive")
|
||||
temperature = self.gas.temperature_from_internal_energy(state.U / state.m)
|
||||
density = state.m / self.compliance_volume
|
||||
pressure = self.gas.pressure(density, temperature)
|
||||
return ThermodynamicProperties(
|
||||
p=pressure,
|
||||
T=temperature,
|
||||
rho=density,
|
||||
u=state.U / state.m,
|
||||
h=self.gas.specific_enthalpy(temperature),
|
||||
)
|
||||
|
||||
def gas_mass_g(self) -> float:
|
||||
return (self.state_1.m + self.state_2.m) * 1.0e3
|
||||
|
||||
def resistance_mass_flow(self) -> float:
|
||||
"""Return center mass flow from port 1 storage to port 2 storage."""
|
||||
port_1 = self.properties_1()
|
||||
port_2 = self.properties_2()
|
||||
pressure_difference = port_1.p - port_2.p
|
||||
if pressure_difference == 0.0:
|
||||
return 0.0
|
||||
upstream = port_1 if pressure_difference > 0.0 else port_2
|
||||
magnitude = self._mass_flow_for_pressure_drop(
|
||||
abs(pressure_difference),
|
||||
density=upstream.rho,
|
||||
temperature=upstream.T,
|
||||
)
|
||||
return magnitude if pressure_difference > 0.0 else -magnitude
|
||||
|
||||
def diagnostics(
|
||||
self,
|
||||
*,
|
||||
mass_flow_kg_s: float,
|
||||
temperature_k: float | None = None,
|
||||
) -> AmesimPnl0001Diagnostics:
|
||||
port_1 = self.properties_1()
|
||||
port_2 = self.properties_2()
|
||||
temperature = temperature_k or (port_1.T if mass_flow_kg_s >= 0.0 else port_2.T)
|
||||
density = port_1.rho if mass_flow_kg_s >= 0.0 else port_2.rho
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
velocity = mass_flow_kg_s / (density * self.area)
|
||||
pressure_drop = self._darcy_pressure_drop(
|
||||
mass_flow_kg_s,
|
||||
density=density,
|
||||
temperature=temperature,
|
||||
)
|
||||
return AmesimPnl0001Diagnostics(
|
||||
mass_flow_kg_s=mass_flow_kg_s,
|
||||
reynolds_number=reynolds,
|
||||
gas_velocity_m_s=velocity,
|
||||
friction_factor=friction_factor,
|
||||
pressure_drop_pa=pressure_drop,
|
||||
)
|
||||
|
||||
def derivatives_from_connections(
|
||||
self,
|
||||
*,
|
||||
port_1_m_flow: float,
|
||||
connected_h_1: float,
|
||||
port_2_m_flow: float,
|
||||
connected_h_2: float,
|
||||
) -> tuple[VolumeState, VolumeState]:
|
||||
port_1 = self.properties_1()
|
||||
port_2 = self.properties_2()
|
||||
center_flow = self.resistance_mass_flow()
|
||||
heat_flow_each = (
|
||||
self.heat_transfer_coefficient
|
||||
* self.heat_transfer_area
|
||||
* (self.external_temperature - 0.5 * (port_1.T + port_2.T))
|
||||
/ 2.0
|
||||
)
|
||||
port_1_external_h = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_1_m_flow,
|
||||
connected_h=connected_h_1,
|
||||
internal_h=port_1.h,
|
||||
)
|
||||
port_2_external_h = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_2_m_flow,
|
||||
connected_h=connected_h_2,
|
||||
internal_h=port_2.h,
|
||||
)
|
||||
port_1_center_h = self.connection_inlet_enthalpy(
|
||||
port_m_flow=-center_flow,
|
||||
connected_h=port_2.h,
|
||||
internal_h=port_1.h,
|
||||
)
|
||||
port_2_center_h = self.connection_inlet_enthalpy(
|
||||
port_m_flow=center_flow,
|
||||
connected_h=port_1.h,
|
||||
internal_h=port_2.h,
|
||||
)
|
||||
return (
|
||||
VolumeState(
|
||||
m=port_1_m_flow - center_flow,
|
||||
U=(
|
||||
port_1_m_flow * port_1_external_h
|
||||
- center_flow * port_1_center_h
|
||||
+ heat_flow_each
|
||||
),
|
||||
),
|
||||
VolumeState(
|
||||
m=port_2_m_flow + center_flow,
|
||||
U=(
|
||||
port_2_m_flow * port_2_external_h
|
||||
+ center_flow * port_2_center_h
|
||||
+ heat_flow_each
|
||||
),
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
class AmesimPnl0002Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
|
||||
"""First-pass AMESim ``PNL0002`` (R-C-R) pipe.
|
||||
|
||||
The center compliance owns the gas state. Positive connection mass flows
|
||||
enter that center storage from each external port.
|
||||
"""
|
||||
|
||||
state_size = 2
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
name: str,
|
||||
*,
|
||||
diameter_mm: float,
|
||||
length_m: float,
|
||||
relative_roughness: float,
|
||||
polytropic_constant: float = 1.35,
|
||||
heat_transfer_coefficient: float = 0.0,
|
||||
external_temperature_k: float = 293.15,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
pctr_0: float = 101_325.0,
|
||||
Tctr_0: float = 293.15,
|
||||
) -> None:
|
||||
if diameter_mm <= 0.0:
|
||||
raise ValueError("diameter_mm must be positive")
|
||||
if length_m <= 0.0:
|
||||
raise ValueError("length_m must be positive")
|
||||
if relative_roughness < 0.0:
|
||||
raise ValueError("relative_roughness must be non-negative")
|
||||
if polytropic_constant <= 0.0:
|
||||
raise ValueError("polytropic_constant must be positive")
|
||||
if heat_transfer_coefficient < 0.0:
|
||||
raise ValueError("heat_transfer_coefficient must be non-negative")
|
||||
if external_temperature_k <= 0.0:
|
||||
raise ValueError("external_temperature_k must be positive")
|
||||
|
||||
super().__init__(name=name)
|
||||
self.diameter = diameter_mm * 1.0e-3
|
||||
self.length = length_m
|
||||
self.relative_roughness = relative_roughness
|
||||
self.polytropic_constant = polytropic_constant
|
||||
self.heat_transfer_coefficient = heat_transfer_coefficient
|
||||
self.external_temperature = external_temperature_k
|
||||
self.gas = gas
|
||||
self.area = diameter_mm_to_area_m2(diameter_mm)
|
||||
self.volume = self.area * self.length
|
||||
self.heat_transfer_area = pi * self.diameter * self.length
|
||||
self._resistance_length = self.length / 2.0
|
||||
|
||||
rho0 = gas.density(pctr_0, Tctr_0)
|
||||
mass0 = rho0 * self.volume
|
||||
self.state = VolumeState(
|
||||
m=mass0,
|
||||
U=mass0 * gas.specific_internal_energy(Tctr_0),
|
||||
)
|
||||
self.port_1 = PortState()
|
||||
self.port_2 = PortState()
|
||||
|
||||
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("pipe mass must stay positive")
|
||||
temperature = self.gas.temperature_from_internal_energy(
|
||||
self.state.U / self.state.m
|
||||
)
|
||||
density = self.state.m / self.volume
|
||||
pressure = self.gas.pressure(density, temperature)
|
||||
properties = ThermodynamicProperties(
|
||||
p=pressure,
|
||||
T=temperature,
|
||||
rho=density,
|
||||
u=self.state.U / self.state.m,
|
||||
h=self.gas.specific_enthalpy(temperature),
|
||||
)
|
||||
self.port_1.p = pressure
|
||||
self.port_1.h_outflow = properties.h
|
||||
self.port_2.p = pressure
|
||||
self.port_2.h_outflow = properties.h
|
||||
return properties
|
||||
|
||||
def gas_mass_g(self) -> float:
|
||||
return self.state.m * 1.0e3
|
||||
|
||||
def port_mass_flow(
|
||||
self,
|
||||
*,
|
||||
port_pressure_pa: float,
|
||||
port_temperature_k: float,
|
||||
) -> float:
|
||||
"""Return mass flow from an external port into the center storage."""
|
||||
if port_pressure_pa <= 0.0:
|
||||
raise ValueError("port_pressure_pa must be positive")
|
||||
if port_temperature_k <= 0.0:
|
||||
raise ValueError("port_temperature_k must be positive")
|
||||
|
||||
center = self.properties()
|
||||
pressure_difference = port_pressure_pa - center.p
|
||||
if pressure_difference == 0.0:
|
||||
return 0.0
|
||||
upstream_pressure = max(port_pressure_pa, center.p)
|
||||
upstream_temperature = (
|
||||
port_temperature_k if pressure_difference > 0.0 else center.T
|
||||
)
|
||||
density = self.gas.density(upstream_pressure, upstream_temperature)
|
||||
magnitude = self._mass_flow_for_resistance_pressure_drop(
|
||||
abs(pressure_difference),
|
||||
density=density,
|
||||
temperature=upstream_temperature,
|
||||
)
|
||||
return magnitude if pressure_difference > 0.0 else -magnitude
|
||||
|
||||
def _mass_flow_for_resistance_pressure_drop(
|
||||
self,
|
||||
pressure_drop_pa: float,
|
||||
*,
|
||||
density: float,
|
||||
temperature: float,
|
||||
) -> float:
|
||||
original_length = self.length
|
||||
self.length = self._resistance_length
|
||||
try:
|
||||
return self._mass_flow_for_pressure_drop(
|
||||
pressure_drop_pa,
|
||||
density=density,
|
||||
temperature=temperature,
|
||||
)
|
||||
finally:
|
||||
self.length = original_length
|
||||
|
||||
def diagnostics(
|
||||
self,
|
||||
*,
|
||||
mass_flow_kg_s: float,
|
||||
temperature_k: float | None = None,
|
||||
) -> AmesimPnl0001Diagnostics:
|
||||
properties = self.properties()
|
||||
temperature = temperature_k or properties.T
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
velocity = mass_flow_kg_s / (properties.rho * self.area)
|
||||
original_length = self.length
|
||||
self.length = self._resistance_length
|
||||
try:
|
||||
pressure_drop = self._darcy_pressure_drop(
|
||||
mass_flow_kg_s,
|
||||
density=properties.rho,
|
||||
temperature=temperature,
|
||||
)
|
||||
finally:
|
||||
self.length = original_length
|
||||
return AmesimPnl0001Diagnostics(
|
||||
mass_flow_kg_s=mass_flow_kg_s,
|
||||
reynolds_number=reynolds,
|
||||
gas_velocity_m_s=velocity,
|
||||
friction_factor=friction_factor,
|
||||
pressure_drop_pa=pressure_drop,
|
||||
)
|
||||
|
||||
def derivatives_from_connections(
|
||||
self,
|
||||
*,
|
||||
port_1_m_flow: float,
|
||||
connected_h_1: float,
|
||||
port_2_m_flow: float,
|
||||
connected_h_2: float,
|
||||
) -> VolumeState:
|
||||
center = self.properties()
|
||||
inlet_h_1 = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_1_m_flow,
|
||||
connected_h=connected_h_1,
|
||||
internal_h=center.h,
|
||||
)
|
||||
inlet_h_2 = self.connection_inlet_enthalpy(
|
||||
port_m_flow=port_2_m_flow,
|
||||
connected_h=connected_h_2,
|
||||
internal_h=center.h,
|
||||
)
|
||||
heat_flow = (
|
||||
self.heat_transfer_coefficient
|
||||
* self.heat_transfer_area
|
||||
* (self.external_temperature - center.T)
|
||||
)
|
||||
return VolumeState(
|
||||
m=port_1_m_flow + port_2_m_flow,
|
||||
U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow,
|
||||
)
|
||||
|
||||
|
||||
class AmesimPnl00rPipe(_DarcyPipeResistanceMixin, AlgebraicComponent):
|
||||
"""First-pass AMESim ``PNL00R`` (R) pipe resistance."""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
name: str,
|
||||
*,
|
||||
diameter_mm: float,
|
||||
length_m: float,
|
||||
relative_roughness: float,
|
||||
gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
|
||||
) -> None:
|
||||
if diameter_mm <= 0.0:
|
||||
raise ValueError("diameter_mm must be positive")
|
||||
if length_m <= 0.0:
|
||||
raise ValueError("length_m must be positive")
|
||||
if relative_roughness < 0.0:
|
||||
raise ValueError("relative_roughness must be non-negative")
|
||||
|
||||
super().__init__(name=name)
|
||||
self.diameter = diameter_mm * 1.0e-3
|
||||
self.length = length_m
|
||||
self.relative_roughness = relative_roughness
|
||||
self.gas = gas
|
||||
self.area = diameter_mm_to_area_m2(diameter_mm)
|
||||
self.port_1 = PortState()
|
||||
self.port_2 = PortState()
|
||||
|
||||
def mass_flow(
|
||||
self,
|
||||
*,
|
||||
port_1_pressure_pa: float,
|
||||
port_1_temperature_k: float,
|
||||
port_2_pressure_pa: float,
|
||||
port_2_temperature_k: float,
|
||||
) -> float:
|
||||
"""Return mass flow from port 1 to port 2 in kg/s."""
|
||||
if port_1_pressure_pa <= 0.0 or port_2_pressure_pa <= 0.0:
|
||||
raise ValueError("port pressures must be positive")
|
||||
if port_1_temperature_k <= 0.0 or port_2_temperature_k <= 0.0:
|
||||
raise ValueError("port temperatures must be positive")
|
||||
pressure_difference = port_1_pressure_pa - port_2_pressure_pa
|
||||
if pressure_difference == 0.0:
|
||||
return 0.0
|
||||
upstream_pressure = max(port_1_pressure_pa, port_2_pressure_pa)
|
||||
upstream_temperature = (
|
||||
port_1_temperature_k
|
||||
if pressure_difference > 0.0
|
||||
else port_2_temperature_k
|
||||
)
|
||||
density = self.gas.density(upstream_pressure, upstream_temperature)
|
||||
magnitude = self._mass_flow_for_pressure_drop(
|
||||
abs(pressure_difference),
|
||||
density=density,
|
||||
temperature=upstream_temperature,
|
||||
)
|
||||
return magnitude if pressure_difference > 0.0 else -magnitude
|
||||
|
||||
def diagnostics(
|
||||
self,
|
||||
*,
|
||||
mass_flow_kg_s: float,
|
||||
pressure_pa: float,
|
||||
temperature_k: float,
|
||||
) -> AmesimPnl0001Diagnostics:
|
||||
density = self.gas.density(pressure_pa, temperature_k)
|
||||
reynolds = self._reynolds_number(mass_flow_kg_s, temperature_k)
|
||||
friction_factor = self._friction_factor(reynolds)
|
||||
velocity = mass_flow_kg_s / (density * self.area)
|
||||
pressure_drop = self._darcy_pressure_drop(
|
||||
mass_flow_kg_s,
|
||||
density=density,
|
||||
temperature=temperature_k,
|
||||
)
|
||||
return AmesimPnl0001Diagnostics(
|
||||
mass_flow_kg_s=mass_flow_kg_s,
|
||||
reynolds_number=reynolds,
|
||||
gas_velocity_m_s=velocity,
|
||||
friction_factor=friction_factor,
|
||||
pressure_drop_pa=pressure_drop,
|
||||
)
|
||||
|
||||
|
||||
def helium_dynamic_viscosity(temperature_k: float) -> float:
|
||||
"""Sutherland approximation centered on the test_mql initial condition."""
|
||||
if temperature_k <= 0.0:
|
||||
raise ValueError("temperature_k must be positive")
|
||||
reference_temperature = 293.15
|
||||
reference_viscosity = 2.0e-5
|
||||
sutherland_constant = 79.4
|
||||
return (
|
||||
reference_viscosity
|
||||
* (temperature_k / reference_temperature) ** 1.5
|
||||
* (reference_temperature + sutherland_constant)
|
||||
/ (temperature_k + sutherland_constant)
|
||||
)
|
||||
Reference in new issue
Block a user