from __future__ import annotations from collections.abc import Mapping from math import isclose, log10, pi, sqrt from app.simulation.core.base import AlgebraicComponent, DynamicComponent, ThermodynamicVolumeComponent from app.simulation.core.catalog import ComponentDisplaySpec, PortDisplaySpec from app.simulation.core.equations import EquationResidual from app.simulation.core.metadata import ( ParameterDefinition, ResultVariableDefinition, THERMODYNAMIC_VOLUME_RESULT_VARIABLES, ) from app.simulation.core.medium import IdealGasMedium, ThermodynamicProperties from app.simulation.core.ports import PortDefinition from app.simulation.core.state import VolumeState 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.1.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = ( ParameterDefinition( "diam", 0.01, label="管径", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "le", 1.0, label="管长", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "rr", 1.0e-5, label="相对粗糙度", quantity="dimensionless", unit="", minimum=0.0, maximum=0.1, ), ParameterDefinition( "gi", 1.0, label="气体类型索引", quantity="dimensionless", unit="", minimum=1.0, maximum=99.0, ), ) 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="pipe", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=20, ) def __init__( self, name: str, medium: IdealGasMedium, *, 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 = self._integer_parameter("gi", 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: IdealGasMedium, 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) if port.h_outflow > 0.0: return max(port.h_outflow / self.medium.cp_ref, 1.0) return self.medium.T_ref @staticmethod def _dynamic_viscosity(temperature_k: float) -> float: if temperature_k <= 0.0: raise ValueError("temperature_k must be positive") reference_temperature = 293.15 reference_viscosity = 1.82e-5 sutherland_constant = 110.4 return ( reference_viscosity * (temperature_k / reference_temperature) ** 1.5 * (reference_temperature + sutherland_constant) / (temperature_k + sutherland_constant) ) 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 turbulent = 1.0 / ( -1.8 * log10((self.rr / 3.7) ** 1.11 + 6.9 / reynolds_number) ) ** 2 if reynolds_number >= 4000.0: return turbulent fraction = (reynolds_number - 2300.0) / 1700.0 return laminar + fraction * (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 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: raise ValueError("unable to bracket PNL00R 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: lower = middle else: upper = middle return 0.5 * (lower + upper) def mass_flow(self, p_1: float, p_2: float) -> float: if p_1 == p_2: 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) / max(self.area * upstream_pressure, 1.0e-18) return { "re": reynolds, "cm": cm, "v": velocity, "ff": self.friction_factor(reynolds), } 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.1.0" PORTS = ( PortDefinition.pneumatic("port_1", nominal_role="bidirectional"), PortDefinition.pneumatic("port_2", nominal_role="bidirectional"), ) PARAMETERS = ( ParameterDefinition( "diam", 0.01, label="管径", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "le", 1.0, label="管长", quantity="length", unit="m", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "rr", 1.0e-5, label="相对粗糙度", quantity="dimensionless", unit="", minimum=0.0, maximum=0.1, ), ParameterDefinition( "k", 1.35, label="多方指数", quantity="dimensionless", unit="", minimum=0.0, minimum_exclusive=True, maximum=2.0, ), ParameterDefinition( "kth", 0.0, label="换热系数", quantity="heat_transfer_coefficient", unit="W/(m2*K)", minimum=0.0, ), ParameterDefinition( "extemp", 293.15, label="外部温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "gi", 1.0, label="气体类型索引", quantity="dimensionless", unit="", minimum=1.0, maximum=99.0, ), ParameterDefinition( "mode", 2.0, label="热模型", quantity="dimensionless", unit="", minimum=1.0, maximum=2.0, ), ParameterDefinition( "p0", 100000.0, label="初始压力", quantity="pressure", unit="Pa", minimum=0.0, minimum_exclusive=True, ), ParameterDefinition( "T0", 293.15, label="初始温度", quantity="temperature", unit="K", minimum=0.0, minimum_exclusive=True, ), ) 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="pipe", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=30, ) def __init__( self, name: str, medium: IdealGasMedium, *, 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 = self._integer_parameter("gi", 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 = self.p0 * self.volume / (medium.R_gas * self.T0) U0 = m0 * medium.specific_internal_energy(self.T0) self.state = VolumeState(m=m0, U=U0) initial_h = medium.specific_enthalpy(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 @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.") return int(rounded) @classmethod def create( cls, *, name: str, medium: IdealGasMedium, 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) @staticmethod def _dynamic_viscosity(temperature_k: float) -> float: return AmesimPnl00r._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: raise ValueError("unable to bracket PNL0001 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: lower = middle else: upper = middle return 0.5 * (lower + upper) def mass_flow(self, p_1: float, p_2: float, temperature: float) -> float: if p_1 == p_2: return 0.0 pressure_difference = p_1 - p_2 upstream_pressure = max(p_1, p_2, 1.0) density = max(self.medium.density(upstream_pressure, temperature), 1.0e-12) magnitude = self._mass_flow_for_pressure_drop( abs(pressure_difference), density=density, temperature=temperature, ) 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) / max(self.area * upstream_pressure, 1.0e-18), "v": flow / (density * self.area), "ff": self.friction_factor(reynolds), } 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.1.0" 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="pipe", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=40, ) @classmethod def create( cls, *, name: str, medium: IdealGasMedium, 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, ) -> float: return self.mass_flow(port_pressure, center_pressure, center_temperature) 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) flow_2 = self.port_mass_flow(self.port_2.p, props.p, props.T) diagnostic_flow = flow_1 if abs(flow_1) >= abs(flow_2) else flow_2 upstream_pressure = max(self.port_1.p, self.port_2.p, props.p, 1.0) density = max(self.medium.density(upstream_pressure, props.T), 1.0e-12) reynolds = self.reynolds_number(diagnostic_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(diagnostic_flow) / max(self.area * upstream_pressure, 1.0e-18), "v": diagnostic_flow / (density * self.area), "ff": self.friction_factor(reynolds), } 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), ), 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), ), ) class AmesimPnl0003(DynamicComponent): """AMESim PNL0003 C-R-C pneumatic pipe with two end compliances.""" state_size = 4 MODEL_TYPE = "amesim_pnl0003" MODEL_VERSION = "0.1.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="pipe", ports=( PortDisplaySpec("port_1", "left", order=10), PortDisplaySpec("port_2", "right", order=20), ), order=50, ) def __init__( self, name: str, medium: IdealGasMedium, *, 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 = AmesimPnl0001._integer_parameter("gi", 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(float(T1_0)) h2 = medium.specific_enthalpy(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: IdealGasMedium, parameters: Mapping[str, float], ) -> "AmesimPnl0003": return cls(name=name, medium=medium, **dict(parameters)) def _initial_state(self, pressure: float, temperature: float) -> VolumeState: mass = pressure * self.compliance_volume / (self.medium.R_gas * temperature) return VolumeState(m=mass, U=mass * self.medium.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(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() @staticmethod def _dynamic_viscosity(temperature_k: float) -> float: return AmesimPnl00r._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: raise ValueError("unable to bracket PNL0003 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: 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 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 _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) / 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": self.friction_factor(reynolds), } 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()]