1197 lines
40 KiB
Python
1197 lines
40 KiB
Python
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()]
|