from __future__ import annotations from dataclasses import dataclass from math import log10, pi from PythonModels.components.amesim_pneumatic import ( HELIUM_PNEUMATIC_GAS, AmesimPneumaticGas, diameter_mm_to_area_m2, ) from PythonModels.core.base import DynamicComponent from PythonModels.core.medium import ThermodynamicProperties from PythonModels.core.ports import PortState from PythonModels.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 AmesimPnl0001Pipe(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 PythonModels convention: positive values enter the pipe storage. AMESim's proprietary pressure-loss calibration is not available in the archive. This implementation therefore uses an explicit Darcy-Weisbach law while preserving 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, 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") 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 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 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 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() inlet_h_1 = self.connection_inlet_enthalpy( port_m_flow=port_1_m_flow, connected_h=connected_h_1, internal_h=internal.h, ) inlet_h_2 = self.connection_inlet_enthalpy( port_m_flow=port_2_m_flow, connected_h=connected_h_2, internal_h=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, ) 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 PNL0001 resistance flow") lower = 0.0 for _ in range(80): 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 _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) 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) )