from __future__ import annotations from collections.abc import Mapping from functools import lru_cache from math import isclose, log, log10, pi, sqrt, tanh from app.simulation.components.amesim.gases import ( AMESIM_GAS_INDEX_PARAMETER, normalize_amesim_gas_index, ) from app.simulation.core.base import AlgebraicComponent, DynamicComponent, ThermodynamicVolumeComponent from app.simulation.core.catalog import ( ComponentDisplaySpec, ParameterGroupDisplaySpec, PortDisplaySpec, ) from app.simulation.core.equations import EquationResidual from app.simulation.core.metadata import ( ParameterCondition, ParameterDefinition, ParameterOption, ResultVariableDefinition, THERMODYNAMIC_VOLUME_RESULT_VARIABLES, ) from app.simulation.core.medium import GasMedium, ThermodynamicProperties from app.simulation.core.ports import PortDefinition from app.simulation.core.state import VolumeState _MAX_REPORTED_FRICTION_FACTOR = 64_000_000.0 def _reported_friction_factor(value: float) -> float: return min(float(value), _MAX_REPORTED_FRICTION_FACTOR) _DYNAMIC_PIPE_POLYTROPIC_MODE = ParameterCondition("mode", (1.0,)) _DYNAMIC_PIPE_HEAT_EXCHANGE_MODE = ParameterCondition("mode", (2.0,)) _DYNAMIC_PIPE_PARAMETER_GROUPS = ( ParameterGroupDisplaySpec( id="thermodynamics", label="热力学", parameters=("k", "kth", "extemp"), order=10, ), ) class AmesimPnl00r(AlgebraicComponent): """AMESim PNL00R pneumatic pipe friction resistance. The public model exposes the AMESim PNL00R catalog/XML contract and uses an auditable Darcy-Weisbach resistance with Reynolds/roughness-dependent friction. Exact `pn2pipefr_` parity is left for the later model tuning pass. """ MODEL_TYPE = "amesim_pnl00r" MODEL_VERSION = "0.3.0" PRESSURE_FLOW_DEPENDS_ON_STREAM = True PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = ( AMESIM_GAS_INDEX_PARAMETER, ParameterDefinition( "diam", 0.01, label="管径", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, description="管路的有效内径，用于计算流通面积和摩擦压降。", ), ParameterDefinition( "le", 1.0, label="管长", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, description="参与摩擦压降计算的管路有效长度。", ), ParameterDefinition( "rr", 1.0e-5, label="相对粗糙度", quantity="dimensionless", unit="", minimum=0.0, maximum=0.1, description="管壁绝对粗糙度与管径之比，用于计算 Darcy 摩擦因子。", ), ) RESULT_VARIABLES = ( ResultVariableDefinition( "re", label="Reynolds 数", quantity="dimensionless", unit="", category="derived", order=10, ), ResultVariableDefinition( "cm", label="质量流量参数", quantity="dimensionless", unit="", category="derived", order=20, ), ResultVariableDefinition( "v", label="平均气体速度", quantity="velocity", unit="m/s", category="derived", order=30, ), ResultVariableDefinition( "ff", label="摩擦因子", quantity="dimensionless", unit="", category="derived", order=40, ), ) DISPLAY = ComponentDisplaySpec( label="PNL00R 气动管路阻力", library_id="amesim", category_id="flow", symbol="amesim_pnl00r", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=20, ) def __init__( self, name: str, medium: GasMedium, *, diam: float = 0.01, le: float = 1.0, rr: float = 1.0e-5, gi: float = 1.0, ) -> None: super().__init__(name=name) self.set_parameter_values({"diam": diam, "le": le, "rr": rr, "gi": gi}) self.medium = medium self.diam = float(diam) self.le = float(le) self.rr = float(rr) self.gi = normalize_amesim_gas_index(gi) self.area = pi * self.diam * self.diam / 4.0 initial_h = medium.specific_enthalpy(medium.T_ref) self.port_1 = self.register_declared_port("port_1") self.port_1.h_outflow = initial_h self.port_2 = self.register_declared_port("port_2") self.port_2.h_outflow = initial_h @staticmethod def _integer_parameter(name: str, value: float) -> int: rounded = round(value) if not isclose(value, rounded, rel_tol=0.0, abs_tol=1.0e-12): raise ValueError(f"PNL00R parameter {name} must be an integer value.") return int(rounded) @classmethod def create( cls, *, name: str, medium: GasMedium, parameters: Mapping[str, float], ) -> AmesimPnl00r: return cls( name=name, medium=medium, diam=parameters["diam"], le=parameters["le"], rr=parameters["rr"], gi=parameters["gi"], ) def _port_temperature(self, port_name: str) -> float: port = self.get_port(port_name) return max( self.medium.temperature_from_pressure_enthalpy( max(port.p, 1.0), port.h_outflow, ), 1.0, ) def _dynamic_viscosity(self, temperature_k: float) -> float: return self.medium.dynamic_viscosity(temperature_k) def reynolds_number(self, mass_flow: float, temperature: float) -> float: viscosity = self._dynamic_viscosity(temperature) return 4.0 * abs(mass_flow) / (pi * self.diam * 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 <= 2300.0: return laminar # pn2pipefr does not apply the fully rough correction at every # turbulent Reynolds number. Its saved ff curves first follow the # hydraulically smooth law and approach the rough asymptote as Re*rr # grows. Keeping those two limits separate reproduces the AMESim # curves for both 14 mm and 20 mm test_mql pipes; putting both terms # directly inside one Haaland logarithm over-predicts PNL0002 friction # by about 23 percent near Re=57,000. smooth_turbulent = 1.0 / ( -1.8 * log10(6.9 / reynolds_number) ) ** 2 if self.rr <= 0.0: turbulent = smooth_turbulent else: # Nikuradse's fully rough asymptote is the Re-independent limit # of Colebrook. Haaland's rounded all-regime approximation is # about 0.2003% high at the test_mql roughness values, enough to # bias its long high-Re PNL0001 filling transient. fully_rough = 1.0 / ( -2.0 * log10(self.rr / 3.7) ) ** 2 roughness_reynolds = reynolds_number * self.rr roughness_weight = roughness_reynolds * roughness_reynolds / ( roughness_reynolds * roughness_reynolds + 180.0 * 180.0 ) turbulent = smooth_turbulent + roughness_weight * ( fully_rough - smooth_turbulent ) if reynolds_number >= 4000.0: return turbulent fraction = (reynolds_number - 2300.0) / 1700.0 return laminar + fraction**0.58 * (turbulent - laminar) def darcy_pressure_drop( self, mass_flow: float, *, density: float, temperature: float, ) -> float: if mass_flow == 0.0: return 0.0 reynolds = self.reynolds_number(mass_flow, temperature) friction = self.friction_factor(reynolds) velocity = mass_flow / (density * self.area) magnitude = ( friction * (self.le / self.diam) * density * velocity * velocity / 2.0 ) return magnitude if mass_flow > 0.0 else -magnitude @lru_cache(maxsize=32768) def _mass_flow_for_pressure_drop( self, pressure_drop: float, *, density: float, temperature: float, ) -> float: if pressure_drop <= 0.0: return 0.0 upper = 1.0e-9 while self.darcy_pressure_drop(upper, density=density, temperature=temperature) < pressure_drop: upper *= 10.0 if upper > 1.0e3: return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) if self.darcy_pressure_drop(middle, density=density, temperature=temperature) < pressure_drop: lower = middle else: upper = middle return 0.5 * (lower + upper) def mass_flow(self, p_1: float, p_2: float) -> float: if isclose(p_1, p_2, rel_tol=0.0, abs_tol=1.0e-8): return 0.0 pressure_difference = p_1 - p_2 upstream_pressure = max(p_1, p_2, 1.0) upstream_temperature = self._port_temperature("port_1" if pressure_difference > 0.0 else "port_2") density = max(self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12) 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 component_result_values(self) -> Mapping[str, float]: m_flow = self.mass_flow(self.port_1.p, self.port_2.p) upstream_pressure = max(self.port_1.p, self.port_2.p, 1.0) upstream_temperature = self._port_temperature( "port_1" if self.port_1.p >= self.port_2.p else "port_2" ) density = max(self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12) reynolds = self.reynolds_number(m_flow, upstream_temperature) velocity = m_flow / (density * self.area) cm = ( abs(m_flow) * sqrt(upstream_temperature) / max(self.area * upstream_pressure, 1.0e-18) ) return { "re": reynolds, "cm": cm, "v": velocity, "ff": _reported_friction_factor(self.friction_factor(reynolds)), } def pressure_flow_equation_values(self) -> tuple[float, ...]: return ( self.port_1.m_flow + self.port_2.m_flow, self.port_1.m_flow - self.mass_flow(self.port_1.p, self.port_2.p), ) def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: return ( EquationResidual( id=f"{self.name}:mass_flow_balance", owner="component", owner_id=self.name, relation="sumToZero", variables=( f"{self.name}.port_1.m_flow", f"{self.name}.port_2.m_flow", ), role="flow", value=self.port_1.m_flow + self.port_2.m_flow, ), EquationResidual( id=f"{self.name}:pressure_flow_relation", owner="component", owner_id=self.name, relation="constitutive", variables=( f"{self.name}.port_1.p", f"{self.name}.port_2.p", f"{self.name}.port_1.m_flow", ), role="flow", value=self.port_1.m_flow - self.mass_flow(self.port_1.p, self.port_2.p), ), ) def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: self.port_1.h_outflow = connected_h["port_2"] self.port_2.h_outflow = connected_h["port_1"] class AmesimPnl0001(ThermodynamicVolumeComponent): """AMESim PNL0001 C-R pneumatic pipe with compressibility and friction.""" MODEL_TYPE = "amesim_pnl0001" MODEL_VERSION = "0.4.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = ( AMESIM_GAS_INDEX_PARAMETER, ParameterDefinition( "diam", 0.01, label="管径", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, description="管路的有效内径，用于计算流通面积、储气容积和摩擦压降。", ), ParameterDefinition( "le", 1.0, label="管长", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, description="管路的有效长度，用于计算储气容积、换热面积和摩擦压降。", ), ParameterDefinition( "rr", 1.0e-5, label="相对粗糙度", quantity="dimensionless", unit="", minimum=0.0, maximum=0.1, description="管壁绝对粗糙度与管径之比，用于计算 Darcy 摩擦因子。", ), ParameterDefinition( "k", 1.35, label="多方指数", quantity="dimensionless", unit="", minimum=0.0, minimum_exclusive=True, maximum=2.0, description=( "mode=1 多方过程使用的指数；当前公开求解器保留该 AMESim " "配置，尚未实现多方指数对状态方程的修正。" ), visible_when=(_DYNAMIC_PIPE_POLYTROPIC_MODE,), ), ParameterDefinition( "kth", 0.0, label="换热系数", quantity="heat_transfer_coefficient", unit="W/(m2*K)", minimum=0.0, description="mode=2 带换热过程使用的气体与外部环境对流换热系数。", visible_when=(_DYNAMIC_PIPE_HEAT_EXCHANGE_MODE,), ), ParameterDefinition( "extemp", 293.15, label="外部温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, description="mode=2 带换热过程使用的外部环境绝对温度。", visible_when=(_DYNAMIC_PIPE_HEAT_EXCHANGE_MODE,), ), ParameterDefinition( "mode", 2.0, label="热模型", quantity="dimensionless", unit="", minimum=1.0, maximum=2.0, editor="choice", options=( ParameterOption(1.0, "多方过程"), ParameterOption(2.0, "带换热"), ), description=( "AMESim 原始编码：1 为多方过程，2 为带换热。当前公开求解器在" "多方模式下关闭环境换热，在带换热模式下按换热系数和外部温度" "计算环境换热。" ), ), ParameterDefinition( "p0", 100000.0, label="初始压力", quantity="pressure", unit="Pa", minimum=0.0, minimum_exclusive=True, description="仿真开始时管内气体的绝对压力。", ), ParameterDefinition( "T0", 293.15, label="初始温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, description="仿真开始时管内气体的绝对温度。", ), ) RESULT_VARIABLES = THERMODYNAMIC_VOLUME_RESULT_VARIABLES + ( ResultVariableDefinition( "re", label="Reynolds 数", quantity="dimensionless", unit="", category="derived", order=100, ), ResultVariableDefinition( "cm", label="质量流量参数", quantity="dimensionless", unit="", category="derived", order=110, ), ResultVariableDefinition( "v", label="平均气体速度", quantity="velocity", unit="m/s", category="derived", order=120, ), ResultVariableDefinition( "ff", label="摩擦因子", quantity="dimensionless", unit="", category="derived", order=130, ), ) DISPLAY = ComponentDisplaySpec( label="PNL0001 C-R 动态管路", library_id="amesim", category_id="flow", symbol="amesim_pnl0001", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=30, parameter_groups=_DYNAMIC_PIPE_PARAMETER_GROUPS, ) def __init__( self, name: str, medium: GasMedium, *, diam: float = 0.01, le: float = 1.0, rr: float = 1.0e-5, k: float = 1.35, kth: float = 0.0, extemp: float = 293.15, gi: float = 1.0, mode: float = 2.0, p0: float = 100000.0, T0: float = 293.15, ) -> None: super().__init__(name=name) self.set_parameter_values( { "diam": diam, "le": le, "rr": rr, "k": k, "kth": kth, "extemp": extemp, "gi": gi, "mode": mode, "p0": p0, "T0": T0, } ) self.medium = medium self.diam = float(diam) self.le = float(le) self.rr = float(rr) self.k = float(k) self.kth = float(kth) self.extemp = float(extemp) self.gi = normalize_amesim_gas_index(gi) self.mode = self._integer_parameter("mode", mode) self.p0 = float(p0) self.T0 = float(T0) self.area = pi * self.diam * self.diam / 4.0 self.volume = self.area * self.le self.exchange_area = pi * self.diam * self.le m0 = medium.density(self.p0, self.T0) * self.volume U0 = m0 * medium.specific_internal_energy_at_pressure(self.p0, self.T0) self.state = VolumeState(m=m0, U=U0) initial_h = medium.specific_enthalpy_at_pressure(self.p0, self.T0) self.port_1 = self.register_declared_port("port_1") self.port_1.p = self.p0 self.port_1.h_outflow = initial_h self.port_2 = self.register_declared_port("port_2") self.port_2.p = self.p0 self.port_2.h_outflow = initial_h self._connected_h: dict[str, float] = {} @staticmethod def _integer_parameter(name: str, value: float) -> int: rounded = round(value) if not isclose(value, rounded, rel_tol=0.0, abs_tol=1.0e-12): raise ValueError(f"PNL0001 parameter {name} must be an integer value.") integer = int(rounded) if name == "mode" and integer not in {1, 2}: raise ValueError("PNL0001 parameter mode must be one of 1, 2.") return integer @classmethod def create( cls, *, name: str, medium: GasMedium, parameters: Mapping[str, float], ) -> "AmesimPnl0001": return cls( name=name, medium=medium, diam=parameters["diam"], le=parameters["le"], rr=parameters["rr"], k=parameters["k"], kth=parameters["kth"], extemp=parameters["extemp"], gi=parameters["gi"], mode=parameters["mode"], p0=parameters["p0"], T0=parameters["T0"], ) 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: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) self.port_1.h_outflow = props.h self.port_2.p = props.p self.port_2.h_outflow = props.h return props def refresh_thermodynamic_ports(self) -> ThermodynamicProperties: return self.properties() def thermal_energy_flow_w(self, temperature: float) -> float: if self.mode == 1: return 0.0 return self.kth * self.exchange_area * (self.extemp - temperature) def _dynamic_viscosity(self, temperature_k: float) -> float: return self.medium.dynamic_viscosity(temperature_k) def reynolds_number(self, mass_flow: float, temperature: float) -> float: viscosity = self._dynamic_viscosity(temperature) return 4.0 * abs(mass_flow) / (pi * self.diam * viscosity) def friction_factor(self, reynolds_number: float) -> float: return AmesimPnl00r.friction_factor(self, reynolds_number) def darcy_pressure_drop( self, mass_flow: float, *, density: float, temperature: float, ) -> float: if mass_flow == 0.0: return 0.0 reynolds = self.reynolds_number(mass_flow, temperature) friction = self.friction_factor(reynolds) velocity = mass_flow / (density * self.area) magnitude = ( friction * (self.le / self.diam) * density * velocity * velocity / 2.0 ) return magnitude if mass_flow > 0.0 else -magnitude def _mass_flow_for_pressure_drop( self, pressure_drop: float, *, density: float, temperature: float, ) -> float: if pressure_drop <= 0.0: return 0.0 upper = 1.0e-9 while self.darcy_pressure_drop( upper, density=density, temperature=temperature, ) < pressure_drop: upper *= 10.0 if upper > 1.0e3: return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) if self.darcy_pressure_drop( middle, density=density, temperature=temperature, ) < pressure_drop: lower = middle else: upper = middle return 0.5 * (lower + upper) @lru_cache(maxsize=32768) def _one_way_pn2pipefr_mass_flow( self, *, upstream_pressure: float, downstream_pressure: float, upstream_temperature: float, resistance_length: float, ) -> float: """AMESim pn2pipefr-style compressible friction flow.""" p_up = max(float(upstream_pressure), 1.0) p_down = max(min(float(downstream_pressure), p_up), 0.0) T_up = max(float(upstream_temperature), 1.0) if resistance_length <= 0.0: raise ValueError("Pipe resistance length must be positive.") gamma_s = self.medium.isentropic_density_pressure_factor( p_up, T_up, p_down, ) gamma_s = min(max(gamma_s, 1.0e-9), 1.0 - 1.0e-9) density = max(self.medium.density(p_up, T_up), 1.0e-12) pressure_ratio = max(p_down / p_up, 0.0) critical_ratio = (2.0 * gamma_s / (gamma_s + 1.0)) ** ( 1.0 / (1.0 - gamma_s) ) def mass_flow_parameter(ratio: float) -> float: if ratio <= critical_ratio: value = ( sqrt(2.0 / (1.0 + gamma_s) * density * T_up / p_up) * (2.0 * gamma_s / (gamma_s + 1.0)) ** (gamma_s / (1.0 - gamma_s)) ) effective_ratio = critical_ratio else: expansion = ratio ** (2.0 * gamma_s) - ratio ** (1.0 + gamma_s) value = sqrt( max( 2.0 / (1.0 - gamma_s) * density * T_up / p_up * expansion, 0.0, ) ) effective_ratio = ratio accuracy = 0.9999 reference_expansion = ( accuracy ** (2.0 * gamma_s) - accuracy ** (1.0 + gamma_s) ) reference = sqrt( max( 2.0 / (1.0 - gamma_s) * density * T_up / p_up * reference_expansion, 0.0, ) ) if value > 0.0 and reference > 0.0: smoothing_argument = ( 12.0 * abs(value / reference) * log(effective_ratio) / log(accuracy) ) value *= tanh(max(smoothing_argument, 0.0)) return value def target_flow(mass_flow: float) -> float: reynolds = self.reynolds_number(mass_flow, T_up) friction = self.friction_factor(reynolds) flow_coefficient = sqrt( self.diam / (resistance_length * friction) ) return ( flow_coefficient * self.area * p_up * mass_flow_parameter(pressure_ratio) / sqrt(T_up) ) flow_coefficient = sqrt(self.diam / (resistance_length * 0.02)) magnitude = ( flow_coefficient * self.area * p_up * mass_flow_parameter(pressure_ratio) / sqrt(T_up) ) for _iteration in range(16): next_magnitude = target_flow(magnitude) if abs(next_magnitude - magnitude) <= max( 1.0e-12, abs(magnitude) * 1.0e-9, ): return next_magnitude magnitude = 0.5 * (magnitude + next_magnitude) return magnitude def mass_flow(self, p_1: float, p_2: float, temperature: float) -> float: if isclose(p_1, p_2, rel_tol=0.0, abs_tol=1.0e-8): return 0.0 pressure_difference = p_1 - p_2 upstream_temperature = max(float(temperature), 1.0) resistance_length = getattr(self, "resistance_length", self.le) magnitude = self._one_way_pn2pipefr_mass_flow( upstream_pressure=max(p_1, p_2), downstream_pressure=min(p_1, p_2), upstream_temperature=upstream_temperature, resistance_length=resistance_length, ) return magnitude if pressure_difference > 0.0 else -magnitude def component_result_values(self) -> Mapping[str, float]: props = self.properties() flow = self.mass_flow(self.port_1.p, props.p, props.T) upstream_pressure = max(self.port_1.p, props.p, 1.0) density = max(self.medium.density(upstream_pressure, props.T), 1.0e-12) reynolds = self.reynolds_number(flow, props.T) return { "m": self.state.m, "U": self.state.U, "p": props.p, "T": props.T, "rho": props.rho, "u": props.u, "h": props.h, "re": reynolds, "cm": ( abs(flow) * sqrt(props.T) / max(self.area * upstream_pressure, 1.0e-18) ), "v": flow / (density * self.area), "ff": _reported_friction_factor(self.friction_factor(reynolds)), } def pressure_flow_equation_values(self) -> tuple[float, ...]: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) return ( self.port_2.p - props.p, self.port_1.m_flow - self.mass_flow(self.port_1.p, props.p, props.T), ) def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) return ( EquationResidual( id=f"{self.name}:port_2_pressure_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_2.p", f"{self.name}.state"), role="effort", value=self.port_2.p - props.p, ), EquationResidual( id=f"{self.name}:port_1_pressure_flow_relation", owner="component", owner_id=self.name, relation="constitutive", variables=( f"{self.name}.port_1.p", f"{self.name}.port_2.p", f"{self.name}.port_1.m_flow", ), role="flow", value=self.port_1.m_flow - self.mass_flow(self.port_1.p, props.p, props.T), ), ) def state_derivative_from_ports( self, connected_h: Mapping[str, float], ) -> list[float]: props = self.properties() inlet_h_1 = self.connection_inlet_enthalpy( port_m_flow=self.port_1.m_flow, connected_h=connected_h["port_1"], internal_h=props.h, ) inlet_h_2 = self.connection_inlet_enthalpy( port_m_flow=self.port_2.m_flow, connected_h=connected_h["port_2"], internal_h=props.h, ) derivative = VolumeState( m=self.port_1.m_flow + self.port_2.m_flow, U=( self.port_1.m_flow * inlet_h_1 + self.port_2.m_flow * inlet_h_2 + self.thermal_energy_flow_w(props.T) ), ) return derivative.as_vector() class AmesimPnl0002(AmesimPnl0001): """AMESim PNL0002 R-C-R pneumatic pipe with one center compliance.""" MODEL_TYPE = "amesim_pnl0002" MODEL_VERSION = "0.6.0" PRESSURE_FLOW_DEPENDS_ON_STREAM = True PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = AmesimPnl0001.PARAMETERS RESULT_VARIABLES = AmesimPnl0001.RESULT_VARIABLES DISPLAY = ComponentDisplaySpec( label="PNL0002 R-C-R 动态管路", library_id="amesim", category_id="flow", symbol="amesim_pnl0002", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=40, parameter_groups=_DYNAMIC_PIPE_PARAMETER_GROUPS, ) @classmethod def create( cls, *, name: str, medium: GasMedium, parameters: Mapping[str, float], ) -> "AmesimPnl0002": return cls( name=name, medium=medium, diam=parameters["diam"], le=parameters["le"], rr=parameters["rr"], k=parameters["k"], kth=parameters["kth"], extemp=parameters["extemp"], gi=parameters["gi"], mode=parameters["mode"], p0=parameters["p0"], T0=parameters["T0"], ) @property def resistance_length(self) -> float: return self.le / 2.0 def properties(self) -> ThermodynamicProperties: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) self.port_1.h_outflow = props.h self.port_2.h_outflow = props.h return props def darcy_pressure_drop( self, mass_flow: float, *, density: float, temperature: float, ) -> float: if mass_flow == 0.0: return 0.0 reynolds = self.reynolds_number(mass_flow, temperature) friction = self.friction_factor(reynolds) velocity = mass_flow / (density * self.area) magnitude = ( friction * (self.resistance_length / self.diam) * density * velocity * velocity / 2.0 ) return magnitude if mass_flow > 0.0 else -magnitude def port_mass_flow( self, port_pressure: float, center_pressure: float, center_temperature: float, *, port_name: str | None = None, ) -> float: upstream_temperature = center_temperature if ( port_name is not None and port_name in self._connected_h and port_pressure > center_pressure ): inlet_h = self._connected_h[port_name] upstream_temperature = self.medium.temperature_from_pressure_enthalpy( max(port_pressure, 1.0), inlet_h, ) return self.mass_flow( port_pressure, center_pressure, max(upstream_temperature, 1.0), ) def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None: self._connected_h = dict(connected_h) def update_flow_temperature_references( self, connected_h: Mapping[str, float], ) -> None: self._connected_h = dict(connected_h) def state_derivative_from_ports( self, connected_h: Mapping[str, float], ) -> list[float]: # Junctions allocate their energy-balanced outlet enthalpy per port. # The separate cache is only the temperature input to pn2pipefr. return super().state_derivative_from_ports(connected_h) def component_result_values(self) -> Mapping[str, float]: props = self.properties() flow_1 = self.port_mass_flow( self.port_1.p, props.p, props.T, port_name="port_1", ) flow_2 = self.port_mass_flow( self.port_2.p, props.p, props.T, port_name="port_2", ) resistance_diagnostics: list[tuple[float, float, float, float]] = [] for port_name, port, flow in ( ("port_1", self.port_1, flow_1), ("port_2", self.port_2, flow_2), ): if flow >= 0.0: upstream_pressure = max(port.p, 1.0) upstream_h = self._connected_h.get(port_name, props.h) upstream_temperature = max( self.medium.temperature_from_pressure_enthalpy( upstream_pressure, upstream_h, ), 1.0, ) else: upstream_pressure = max(props.p, 1.0) upstream_temperature = props.T density = max( self.medium.density(upstream_pressure, upstream_temperature), 1.0e-12, ) reynolds = self.reynolds_number(flow, upstream_temperature) resistance_diagnostics.append( ( reynolds, ( abs(flow) * sqrt(upstream_temperature) / max(self.area * upstream_pressure, 1.0e-18) ), abs(flow) / (density * self.area), self.friction_factor(reynolds), ) ) reynolds, cm, velocity, friction = ( sum(values) / len(resistance_diagnostics) for values in zip(*resistance_diagnostics) ) return { "m": self.state.m, "U": self.state.U, "p": props.p, "T": props.T, "rho": props.rho, "u": props.u, "h": props.h, "re": reynolds, "cm": cm, "v": velocity, "ff": _reported_friction_factor(friction), } def pressure_flow_equation_values(self) -> tuple[float, ...]: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) return ( self.port_1.m_flow - self.port_mass_flow( self.port_1.p, props.p, props.T, port_name="port_1", ), self.port_2.m_flow - self.port_mass_flow( self.port_2.p, props.p, props.T, port_name="port_2", ), ) def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume) return ( EquationResidual( id=f"{self.name}:port_1_pressure_flow_relation", owner="component", owner_id=self.name, relation="constitutive", variables=( f"{self.name}.port_1.p", f"{self.name}.state", f"{self.name}.port_1.m_flow", ), role="flow", value=self.port_1.m_flow - self.port_mass_flow( self.port_1.p, props.p, props.T, port_name="port_1", ), ), EquationResidual( id=f"{self.name}:port_2_pressure_flow_relation", owner="component", owner_id=self.name, relation="constitutive", variables=( f"{self.name}.port_2.p", f"{self.name}.state", f"{self.name}.port_2.m_flow", ), role="flow", value=self.port_2.m_flow - self.port_mass_flow( self.port_2.p, props.p, props.T, port_name="port_2", ), ), ) class AmesimPnl0003(DynamicComponent): """AMESim PNL0003 C-R-C pneumatic pipe with two end compliances.""" state_size = 4 MODEL_TYPE = "amesim_pnl0003" MODEL_VERSION = "0.4.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = AmesimPnl0001.PARAMETERS[:-2] + ( ParameterDefinition( "p1_0", 100000.0, label="端口 1 初始压力", quantity="pressure", unit="Pa", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "T1_0", 293.15, label="端口 1 初始温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "p2_0", 100000.0, label="端口 2 初始压力", quantity="pressure", unit="Pa", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "T2_0", 293.15, label="端口 2 初始温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, ), ) RESULT_VARIABLES = ( ResultVariableDefinition("m1", "端口 1 侧质量", "mass", "kg", "state", 10), ResultVariableDefinition("U1", "端口 1 侧内能", "internal_energy", "J", "state", 20), ResultVariableDefinition("p1", "端口 1 侧压力", "pressure", "Pa", "thermodynamic", 30), ResultVariableDefinition("T1", "端口 1 侧温度", "temperature", "K", "thermodynamic", 40), ResultVariableDefinition("rho1", "端口 1 侧密度", "density", "kg/m³", "thermodynamic", 50), ResultVariableDefinition("u1", "端口 1 侧比内能", "specific_internal_energy", "J/kg", "thermodynamic", 60), ResultVariableDefinition("h1", "端口 1 侧比焓", "specific_enthalpy", "J/kg", "thermodynamic", 70), ResultVariableDefinition("m2", "端口 2 侧质量", "mass", "kg", "state", 80), ResultVariableDefinition("U2", "端口 2 侧内能", "internal_energy", "J", "state", 90), ResultVariableDefinition("p2", "端口 2 侧压力", "pressure", "Pa", "thermodynamic", 100), ResultVariableDefinition("T2", "端口 2 侧温度", "temperature", "K", "thermodynamic", 110), ResultVariableDefinition("rho2", "端口 2 侧密度", "density", "kg/m³", "thermodynamic", 120), ResultVariableDefinition("u2", "端口 2 侧比内能", "specific_internal_energy", "J/kg", "thermodynamic", 130), ResultVariableDefinition("h2", "端口 2 侧比焓", "specific_enthalpy", "J/kg", "thermodynamic", 140), ResultVariableDefinition("dmctr", "中心质量流量", "mass_flow", "kg/s", "derived", 150), ResultVariableDefinition("re", "Reynolds 数", "dimensionless", "", "derived", 160), ResultVariableDefinition("cm", "质量流量参数", "dimensionless", "", "derived", 170), ResultVariableDefinition("v", "平均气体速度", "velocity", "m/s", "derived", 180), ResultVariableDefinition("ff", "摩擦因子", "dimensionless", "", "derived", 190), ) DISPLAY = ComponentDisplaySpec( label="PNL0003 C-R-C 动态管路", library_id="amesim", category_id="flow", symbol="amesim_pnl0003", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=50, parameter_groups=_DYNAMIC_PIPE_PARAMETER_GROUPS, ) def __init__( self, name: str, medium: GasMedium, *, diam: float = 0.01, le: float = 1.0, rr: float = 1.0e-5, k: float = 1.35, kth: float = 0.0, extemp: float = 293.15, gi: float = 1.0, mode: float = 2.0, p1_0: float = 100000.0, T1_0: float = 293.15, p2_0: float = 100000.0, T2_0: float = 293.15, ) -> None: super().__init__(name=name) self.set_parameter_values( { "diam": diam, "le": le, "rr": rr, "k": k, "kth": kth, "extemp": extemp, "gi": gi, "mode": mode, "p1_0": p1_0, "T1_0": T1_0, "p2_0": p2_0, "T2_0": T2_0, } ) self.medium = medium self.diam = float(diam) self.le = float(le) self.rr = float(rr) self.k = float(k) self.kth = float(kth) self.extemp = float(extemp) self.gi = normalize_amesim_gas_index(gi) self.mode = AmesimPnl0001._integer_parameter("mode", mode) self.area = pi * self.diam * self.diam / 4.0 self.volume = self.area * self.le self.compliance_volume = self.volume / 2.0 self.exchange_area = pi * self.diam * self.le self.state_1 = self._initial_state(float(p1_0), float(T1_0)) self.state_2 = self._initial_state(float(p2_0), float(T2_0)) h1 = medium.specific_enthalpy_at_pressure(float(p1_0), float(T1_0)) h2 = medium.specific_enthalpy_at_pressure(float(p2_0), float(T2_0)) self.port_1 = self.register_declared_port("port_1") self.port_1.p = float(p1_0) self.port_1.h_outflow = h1 self.port_2 = self.register_declared_port("port_2") self.port_2.p = float(p2_0) self.port_2.h_outflow = h2 @classmethod def create( cls, *, name: str, medium: GasMedium, parameters: Mapping[str, float], ) -> "AmesimPnl0003": return cls(name=name, medium=medium, **dict(parameters)) def _initial_state(self, pressure: float, temperature: float) -> VolumeState: mass = self.medium.density(pressure, temperature) * self.compliance_volume return VolumeState( m=mass, U=mass * self.medium.specific_internal_energy_at_pressure( pressure, 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(self, state: VolumeState) -> ThermodynamicProperties: return self.medium.properties_from_mU(state.m, state.U, self.compliance_volume) def properties_1(self) -> ThermodynamicProperties: props = self._properties(self.state_1) self.port_1.p = props.p self.port_1.h_outflow = props.h return props def properties_2(self) -> ThermodynamicProperties: props = self._properties(self.state_2) self.port_2.p = props.p self.port_2.h_outflow = props.h return props def refresh_thermodynamic_ports(self) -> tuple[ThermodynamicProperties, ThermodynamicProperties]: return self.properties_1(), self.properties_2() def _dynamic_viscosity(self, temperature_k: float) -> float: return self.medium.dynamic_viscosity(temperature_k) def reynolds_number(self, mass_flow: float, temperature: float) -> float: viscosity = self._dynamic_viscosity(temperature) return 4.0 * abs(mass_flow) / (pi * self.diam * viscosity) def friction_factor(self, reynolds_number: float) -> float: return AmesimPnl00r.friction_factor(self, reynolds_number) def darcy_pressure_drop( self, mass_flow: float, *, density: float, temperature: float, ) -> float: if mass_flow == 0.0: return 0.0 reynolds = self.reynolds_number(mass_flow, temperature) friction = self.friction_factor(reynolds) velocity = mass_flow / (density * self.area) magnitude = friction * (self.le / self.diam) * density * velocity * velocity / 2.0 return magnitude if mass_flow > 0.0 else -magnitude @lru_cache(maxsize=32768) def _mass_flow_for_pressure_drop( self, pressure_drop: float, *, density: float, temperature: float, ) -> float: if pressure_drop <= 0.0: return 0.0 upper = 1.0e-9 while self.darcy_pressure_drop(upper, density=density, temperature=temperature) < pressure_drop: upper *= 10.0 if upper > 1.0e3: return 1.0e3 lower = 0.0 for _ in range(48): middle = 0.5 * (lower + upper) if self.darcy_pressure_drop(middle, density=density, temperature=temperature) < pressure_drop: lower = middle else: upper = middle return 0.5 * (lower + upper) def resistance_mass_flow(self) -> float: port_1 = self._properties(self.state_1) port_2 = self._properties(self.state_2) pressure_difference = port_1.p - port_2.p if isclose(port_1.p, port_2.p, rel_tol=0.0, abs_tol=1.0e-8): 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 _heat_flow_each(self, temperature_1: float, temperature_2: float) -> float: if self.mode == 1: return 0.0 return self.kth * self.exchange_area * (self.extemp - 0.5 * (temperature_1 + temperature_2)) / 2.0 def component_result_values(self) -> Mapping[str, float]: port_1 = self.properties_1() port_2 = self.properties_2() center_flow = self.resistance_mass_flow() upstream = port_1 if center_flow >= 0.0 else port_2 reynolds = self.reynolds_number(center_flow, upstream.T) return { "m1": self.state_1.m, "U1": self.state_1.U, "p1": port_1.p, "T1": port_1.T, "rho1": port_1.rho, "u1": port_1.u, "h1": port_1.h, "m2": self.state_2.m, "U2": self.state_2.U, "p2": port_2.p, "T2": port_2.T, "rho2": port_2.rho, "u2": port_2.u, "h2": port_2.h, "dmctr": center_flow, "re": reynolds, "cm": ( abs(center_flow) * sqrt(upstream.T) / max(self.area * max(port_1.p, port_2.p, 1.0), 1.0e-18) ), "v": center_flow / (max(upstream.rho, 1.0e-12) * self.area), "ff": _reported_friction_factor(self.friction_factor(reynolds)), } def pressure_flow_equation_values(self) -> tuple[float, ...]: port_1 = self._properties(self.state_1) port_2 = self._properties(self.state_2) return ( self.port_1.p - port_1.p, self.port_2.p - port_2.p, ) def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]: port_1 = self._properties(self.state_1) port_2 = self._properties(self.state_2) return ( EquationResidual( id=f"{self.name}:port_1_pressure_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_1.p", f"{self.name}.state"), role="effort", value=self.port_1.p - port_1.p, ), EquationResidual( id=f"{self.name}:port_2_pressure_state", owner="component", owner_id=self.name, relation="state", variables=(f"{self.name}.port_2.p", f"{self.name}.state"), role="effort", value=self.port_2.p - port_2.p, ), ) def state_derivative_from_ports(self, connected_h: Mapping[str, float]) -> list[float]: port_1 = self.properties_1() port_2 = self.properties_2() center_flow = self.resistance_mass_flow() heat_flow_each = self._heat_flow_each(port_1.T, port_2.T) port_1_external_h = self.connection_inlet_enthalpy( port_m_flow=self.port_1.m_flow, connected_h=connected_h["port_1"], internal_h=port_1.h, ) port_2_external_h = self.connection_inlet_enthalpy( port_m_flow=self.port_2.m_flow, connected_h=connected_h["port_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, ) d1 = VolumeState( m=self.port_1.m_flow - center_flow, U=self.port_1.m_flow * port_1_external_h - center_flow * port_1_center_h + heat_flow_each, ) d2 = VolumeState( m=self.port_2.m_flow + center_flow, U=self.port_2.m_flow * port_2_external_h + center_flow * port_2_center_h + heat_flow_each, ) return [*d1.as_vector(), *d2.as_vector()]