实现test_mql PNL0003与PNL00R管路

2026-07-17 08:54:38 +00:00
1 parent c4c1b35dee
commit 568be2e632
@@ -8,7 +8,7 @@ from PythonModels.components.amesim_pneumatic import (
    AmesimPneumaticGas,
    diameter_mm_to_area_m2,
 )
-from PythonModels.core.base import DynamicComponent
+from PythonModels.core.base import AlgebraicComponent, DynamicComponent
 from PythonModels.core.medium import ThermodynamicProperties
 from PythonModels.core.ports import PortState
 from PythonModels.core.state import VolumeState
@@ -23,7 +23,89 @@ class AmesimPnl0001Diagnostics:
    pressure_drop_pa: float


-class AmesimPnl0001Pipe(DynamicComponent):
+class _DarcyPipeResistanceMixin:
+    diameter: float
+    length: float
+    relative_roughness: float
+    area: float
+
+    def _mass_flow_for_pressure_drop(
+        self,
+        pressure_drop_pa: float,
+        *,
+        density: float,
+        temperature: float,
+    ) -> float:
+        if pressure_drop_pa <= 0.0:
+            return 0.0
+        upper = 1.0e-9
+        while self._darcy_pressure_drop(
+            upper,
+            density=density,
+            temperature=temperature,
+        ) < pressure_drop_pa:
+            upper *= 10.0
+            if upper > 1.0e3:
+                raise ValueError("unable to bracket pneumatic pipe resistance flow")
+        lower = 0.0
+        for _ in range(80):
+            middle = 0.5 * (lower + upper)
+            if self._darcy_pressure_drop(
+                middle,
+                density=density,
+                temperature=temperature,
+            ) < pressure_drop_pa:
+                lower = middle
+            else:
+                upper = middle
+        return 0.5 * (lower + upper)
+
+    def _darcy_pressure_drop(
+        self,
+        mass_flow_kg_s: float,
+        *,
+        density: float,
+        temperature: float,
+    ) -> float:
+        if mass_flow_kg_s == 0.0:
+            return 0.0
+        reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
+        friction_factor = self._friction_factor(reynolds)
+        velocity = mass_flow_kg_s / (density * self.area)
+        magnitude = (
+            friction_factor
+            * (self.length / self.diameter)
+            * density
+            * velocity
+            * velocity
+            / 2.0
+        )
+        return magnitude if mass_flow_kg_s > 0.0 else -magnitude
+
+    def _reynolds_number(self, mass_flow_kg_s: float, temperature: float) -> float:
+        viscosity = helium_dynamic_viscosity(temperature)
+        return 4.0 * abs(mass_flow_kg_s) / (pi * self.diameter * viscosity)
+
+    def _friction_factor(self, reynolds_number: float) -> float:
+        if reynolds_number <= 0.0:
+            return 64_000_000.0
+        laminar = 64.0 / reynolds_number
+        if reynolds_number <= 2_300.0:
+            return laminar
+        turbulent = 1.0 / (
+            -1.8
+            * log10(
+                (self.relative_roughness / 3.7) ** 1.11
+                + 6.9 / reynolds_number
+            )
+        ) ** 2
+        if reynolds_number >= 4_000.0:
+            return turbulent
+        fraction = (reynolds_number - 2_300.0) / 1_700.0
+        return laminar + fraction * (turbulent - laminar)
+
+
+class AmesimPnl0001Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
    """Physical first-pass implementation of AMESim ``PNL0001`` (C-R).

    Port 2 owns the lumped gas storage. Port 1 is connected through a Darcy
@@ -193,80 +275,289 @@ class AmesimPnl0001Pipe(DynamicComponent):
            U=port_1_m_flow * inlet_h_1 + port_2_m_flow * inlet_h_2 + heat_flow,
        )

-    def _mass_flow_for_pressure_drop(
-        self,
-        pressure_drop_pa: float,
-        *,
-        density: float,
-        temperature: float,
-    ) -> float:
-        if pressure_drop_pa <= 0.0:
-            return 0.0
-        upper = 1.0e-9
-        while self._darcy_pressure_drop(
-            upper,
-            density=density,
-            temperature=temperature,
-        ) < pressure_drop_pa:
-            upper *= 10.0
-            if upper > 1.0e3:
-                raise ValueError("unable to bracket PNL0001 resistance flow")
-        lower = 0.0
-        for _ in range(80):
-            middle = 0.5 * (lower + upper)
-            if self._darcy_pressure_drop(
-                middle,
-                density=density,
-                temperature=temperature,
-            ) < pressure_drop_pa:
-                lower = middle
-            else:
-                upper = middle
-        return 0.5 * (lower + upper)

-    def _darcy_pressure_drop(
+class AmesimPnl0003Pipe(_DarcyPipeResistanceMixin, DynamicComponent):
+    """First-pass AMESim ``PNL0003`` (C-R-C) pipe.
+
+    The two pipe-end compliances are represented as equal half-volume gas
+    stores connected by the same auditable Darcy resistance used for PNL0001.
+    Center flow is positive from port 1 storage to port 2 storage.
+    """
+
+    state_size = 4
+
+    def __init__(
        self,
-        mass_flow_kg_s: float,
+        name: str,
        *,
-        density: float,
-        temperature: float,
-    ) -> float:
-        if mass_flow_kg_s == 0.0:
+        diameter_mm: float,
+        length_m: float,
+        relative_roughness: float,
+        polytropic_constant: float = 1.35,
+        heat_transfer_coefficient: float = 0.0,
+        external_temperature_k: float = 293.15,
+        gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
+        p1_0: float = 101_325.0,
+        T1_0: float = 293.15,
+        p2_0: float = 101_325.0,
+        T2_0: float = 293.15,
+    ) -> None:
+        if diameter_mm <= 0.0:
+            raise ValueError("diameter_mm must be positive")
+        if length_m <= 0.0:
+            raise ValueError("length_m must be positive")
+        if relative_roughness < 0.0:
+            raise ValueError("relative_roughness must be non-negative")
+        if polytropic_constant <= 0.0:
+            raise ValueError("polytropic_constant must be positive")
+        if heat_transfer_coefficient < 0.0:
+            raise ValueError("heat_transfer_coefficient must be non-negative")
+        if external_temperature_k <= 0.0:
+            raise ValueError("external_temperature_k must be positive")
+
+        super().__init__(name=name)
+        self.diameter = diameter_mm * 1.0e-3
+        self.length = length_m
+        self.relative_roughness = relative_roughness
+        self.polytropic_constant = polytropic_constant
+        self.heat_transfer_coefficient = heat_transfer_coefficient
+        self.external_temperature = external_temperature_k
+        self.gas = gas
+        self.area = diameter_mm_to_area_m2(diameter_mm)
+        self.volume = self.area * self.length
+        self.compliance_volume = self.volume / 2.0
+        self.heat_transfer_area = pi * self.diameter * self.length
+
+        self.state_1 = self._initial_state(p1_0, T1_0)
+        self.state_2 = self._initial_state(p2_0, T2_0)
+        self.port_1 = PortState()
+        self.port_2 = PortState()
+
+    def _initial_state(self, pressure: float, temperature: float) -> VolumeState:
+        rho = self.gas.density(pressure, temperature)
+        mass = rho * self.compliance_volume
+        return VolumeState(
+            m=mass,
+            U=mass * self.gas.specific_internal_energy(temperature),
+        )
+
+    def get_state_vector(self) -> list[float]:
+        return [*self.state_1.as_vector(), *self.state_2.as_vector()]
+
+    def set_state_vector(self, values: list[float]) -> None:
+        if len(values) != 4:
+            raise ValueError("PNL0003 state vector requires four values")
+        self.state_1 = VolumeState.from_vector(values[:2])
+        self.state_2 = VolumeState.from_vector(values[2:])
+
+    def properties_1(self) -> ThermodynamicProperties:
+        properties = self._properties(self.state_1)
+        self.port_1.p = properties.p
+        self.port_1.h_outflow = properties.h
+        return properties
+
+    def properties_2(self) -> ThermodynamicProperties:
+        properties = self._properties(self.state_2)
+        self.port_2.p = properties.p
+        self.port_2.h_outflow = properties.h
+        return properties
+
+    def _properties(self, state: VolumeState) -> ThermodynamicProperties:
+        if state.m <= 0.0:
+            raise ValueError("pipe mass must stay positive")
+        temperature = self.gas.temperature_from_internal_energy(state.U / state.m)
+        density = state.m / self.compliance_volume
+        pressure = self.gas.pressure(density, temperature)
+        return ThermodynamicProperties(
+            p=pressure,
+            T=temperature,
+            rho=density,
+            u=state.U / state.m,
+            h=self.gas.specific_enthalpy(temperature),
+        )
+
+    def gas_mass_g(self) -> float:
+        return (self.state_1.m + self.state_2.m) * 1.0e3
+
+    def resistance_mass_flow(self) -> float:
+        """Return center mass flow from port 1 storage to port 2 storage."""
+        port_1 = self.properties_1()
+        port_2 = self.properties_2()
+        pressure_difference = port_1.p - port_2.p
+        if pressure_difference == 0.0:
            return 0.0
+        upstream = port_1 if pressure_difference > 0.0 else port_2
+        magnitude = self._mass_flow_for_pressure_drop(
+            abs(pressure_difference),
+            density=upstream.rho,
+            temperature=upstream.T,
+        )
+        return magnitude if pressure_difference > 0.0 else -magnitude
+
+    def diagnostics(
+        self,
+        *,
+        mass_flow_kg_s: float,
+        temperature_k: float | None = None,
+    ) -> AmesimPnl0001Diagnostics:
+        port_1 = self.properties_1()
+        port_2 = self.properties_2()
+        temperature = temperature_k or (port_1.T if mass_flow_kg_s >= 0.0 else port_2.T)
+        density = port_1.rho if mass_flow_kg_s >= 0.0 else port_2.rho
        reynolds = self._reynolds_number(mass_flow_kg_s, temperature)
        friction_factor = self._friction_factor(reynolds)
        velocity = mass_flow_kg_s / (density * self.area)
-        magnitude = (
-            friction_factor
-            * (self.length / self.diameter)
-            * density
-            * velocity
-            * velocity
+        pressure_drop = self._darcy_pressure_drop(
+            mass_flow_kg_s,
+            density=density,
+            temperature=temperature,
+        )
+        return AmesimPnl0001Diagnostics(
+            mass_flow_kg_s=mass_flow_kg_s,
+            reynolds_number=reynolds,
+            gas_velocity_m_s=velocity,
+            friction_factor=friction_factor,
+            pressure_drop_pa=pressure_drop,
+        )
+
+    def derivatives_from_connections(
+        self,
+        *,
+        port_1_m_flow: float,
+        connected_h_1: float,
+        port_2_m_flow: float,
+        connected_h_2: float,
+    ) -> tuple[VolumeState, VolumeState]:
+        port_1 = self.properties_1()
+        port_2 = self.properties_2()
+        center_flow = self.resistance_mass_flow()
+        heat_flow_each = (
+            self.heat_transfer_coefficient
+            * self.heat_transfer_area
+            * (self.external_temperature - 0.5 * (port_1.T + port_2.T))
            / 2.0
        )
-        return magnitude if mass_flow_kg_s > 0.0 else -magnitude
+        port_1_external_h = self.connection_inlet_enthalpy(
+            port_m_flow=port_1_m_flow,
+            connected_h=connected_h_1,
+            internal_h=port_1.h,
+        )
+        port_2_external_h = self.connection_inlet_enthalpy(
+            port_m_flow=port_2_m_flow,
+            connected_h=connected_h_2,
+            internal_h=port_2.h,
+        )
+        port_1_center_h = self.connection_inlet_enthalpy(
+            port_m_flow=-center_flow,
+            connected_h=port_2.h,
+            internal_h=port_1.h,
+        )
+        port_2_center_h = self.connection_inlet_enthalpy(
+            port_m_flow=center_flow,
+            connected_h=port_1.h,
+            internal_h=port_2.h,
+        )
+        return (
+            VolumeState(
+                m=port_1_m_flow - center_flow,
+                U=(
+                    port_1_m_flow * port_1_external_h
+                    - center_flow * port_1_center_h
+                    + heat_flow_each
+                ),
+            ),
+            VolumeState(
+                m=port_2_m_flow + center_flow,
+                U=(
+                    port_2_m_flow * port_2_external_h
+                    + center_flow * port_2_center_h
+                    + heat_flow_each
+                ),
+            ),
+        )

-    def _reynolds_number(self, mass_flow_kg_s: float, temperature: float) -> float:
-        viscosity = helium_dynamic_viscosity(temperature)
-        return 4.0 * abs(mass_flow_kg_s) / (pi * self.diameter * viscosity)

-    def _friction_factor(self, reynolds_number: float) -> float:
-        if reynolds_number <= 0.0:
-            return 64_000_000.0
-        laminar = 64.0 / reynolds_number
-        if reynolds_number <= 2_300.0:
-            return laminar
-        turbulent = 1.0 / (
-            -1.8
-            * log10(
-                (self.relative_roughness / 3.7) ** 1.11
-                + 6.9 / reynolds_number
-            )
-        ) ** 2
-        if reynolds_number >= 4_000.0:
-            return turbulent
-        fraction = (reynolds_number - 2_300.0) / 1_700.0
-        return laminar + fraction * (turbulent - laminar)
+class AmesimPnl00rPipe(_DarcyPipeResistanceMixin, AlgebraicComponent):
+    """First-pass AMESim ``PNL00R`` (R) pipe resistance."""
+
+    def __init__(
+        self,
+        name: str,
+        *,
+        diameter_mm: float,
+        length_m: float,
+        relative_roughness: float,
+        gas: AmesimPneumaticGas = HELIUM_PNEUMATIC_GAS,
+    ) -> None:
+        if diameter_mm <= 0.0:
+            raise ValueError("diameter_mm must be positive")
+        if length_m <= 0.0:
+            raise ValueError("length_m must be positive")
+        if relative_roughness < 0.0:
+            raise ValueError("relative_roughness must be non-negative")
+
+        super().__init__(name=name)
+        self.diameter = diameter_mm * 1.0e-3
+        self.length = length_m
+        self.relative_roughness = relative_roughness
+        self.gas = gas
+        self.area = diameter_mm_to_area_m2(diameter_mm)
+        self.port_1 = PortState()
+        self.port_2 = PortState()
+
+    def mass_flow(
+        self,
+        *,
+        port_1_pressure_pa: float,
+        port_1_temperature_k: float,
+        port_2_pressure_pa: float,
+        port_2_temperature_k: float,
+    ) -> float:
+        """Return mass flow from port 1 to port 2 in kg/s."""
+        if port_1_pressure_pa <= 0.0 or port_2_pressure_pa <= 0.0:
+            raise ValueError("port pressures must be positive")
+        if port_1_temperature_k <= 0.0 or port_2_temperature_k <= 0.0:
+            raise ValueError("port temperatures must be positive")
+        pressure_difference = port_1_pressure_pa - port_2_pressure_pa
+        if pressure_difference == 0.0:
+            return 0.0
+        upstream_pressure = max(port_1_pressure_pa, port_2_pressure_pa)
+        upstream_temperature = (
+            port_1_temperature_k
+            if pressure_difference > 0.0
+            else port_2_temperature_k
+        )
+        density = self.gas.density(upstream_pressure, upstream_temperature)
+        magnitude = self._mass_flow_for_pressure_drop(
+            abs(pressure_difference),
+            density=density,
+            temperature=upstream_temperature,
+        )
+        return magnitude if pressure_difference > 0.0 else -magnitude
+
+    def diagnostics(
+        self,
+        *,
+        mass_flow_kg_s: float,
+        pressure_pa: float,
+        temperature_k: float,
+    ) -> AmesimPnl0001Diagnostics:
+        density = self.gas.density(pressure_pa, temperature_k)
+        reynolds = self._reynolds_number(mass_flow_kg_s, temperature_k)
+        friction_factor = self._friction_factor(reynolds)
+        velocity = mass_flow_kg_s / (density * self.area)
+        pressure_drop = self._darcy_pressure_drop(
+            mass_flow_kg_s,
+            density=density,
+            temperature=temperature_k,
+        )
+        return AmesimPnl0001Diagnostics(
+            mass_flow_kg_s=mass_flow_kg_s,
+            reynolds_number=reynolds,
+            gas_velocity_m_s=velocity,
+            friction_factor=friction_factor,
+            pressure_drop_pa=pressure_drop,
+        )


 def helium_dynamic_viscosity(temperature_k: float) -> float: