from __future__ import annotations from dataclasses import dataclass from math import log10, pi, sqrt from PythonModels.components.amesim_pneumatic import ( HELIUM_PNEUMATIC_GAS, AmesimPneumaticGas, compressible_orifice_mass_flow, diameter_mm_to_area_m2, ) from PythonModels.core.base import AlgebraicComponent, 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 _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 PythonModels 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() # PNL0001 is a fixed-volume distributed line store. Its transported # energy variable therefore follows specific internal energy, not the # chamber-style stagnation enthalpy contract. For this ideal gas, # h = gamma * u. Outflow always carries the local u. 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, ) 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) )