diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py index 12d4b5b..8e72662 100644 --- a/app/simulation/components/amesim/flow/pipes.py +++ b/app/simulation/components/amesim/flow/pipes.py @@ -1667,6 +1667,16 @@ class AmesimPnl0003(DynamicComponent): viscosity = self._dynamic_viscosity(temperature) return 4.0 * abs(mass_flow) / (pi * self.diam * viscosity) + def reported_reynolds_number( + self, + mass_flow: float, + temperature: float, + ) -> float: + """Return AMESim's derived re without changing pipe dynamics.""" + + viscosity = self.medium.diagnostic_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) @@ -1734,6 +1744,7 @@ class AmesimPnl0003(DynamicComponent): 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) + reported_reynolds = self.reported_reynolds_number(center_flow, upstream.T) friction = self.friction_factor(reynolds) return { "m1": self.state_1.m, @@ -1751,7 +1762,7 @@ class AmesimPnl0003(DynamicComponent): "u2": port_2.u, "h2": port_2.h, "dmctr": center_flow, - "re": reynolds, + "re": reported_reynolds, "cm": _bare_pipe_mass_flow_parameter( center_flow, area=self.area, diff --git a/app/simulation/components/amesim/media/mediums.py b/app/simulation/components/amesim/media/mediums.py index 8618e13..0256f4f 100644 --- a/app/simulation/components/amesim/media/mediums.py +++ b/app/simulation/components/amesim/media/mediums.py @@ -2,7 +2,7 @@ from __future__ import annotations from collections.abc import Callable, Sequence from dataclasses import dataclass -from math import isfinite +from math import exp, isfinite, log from typing import ClassVar from app.simulation.core.errors import RecoverableTrialStateError @@ -55,6 +55,12 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): fluid: ClassVar[PengRobinsonFluid] = HELIUM_PR nasa_cp_over_R: ClassVar[float] = 2.5 nasa_enthalpy_constant_K: ClassVar[float] = -745.375 + nasa_viscosity_coefficients: ClassVar[tuple[float, float, float, float]] = ( + 0.7501594, + 35.76324, + -2212.129, + 0.9212635, + ) name: str = "AMESimHeliumPengRobinson" R_gas: float = HELIUM_PR.specific_gas_constant @@ -73,6 +79,19 @@ class AmesimHeliumPengRobinsonMedium(IdealGasMedium): del T return self.cv + def diagnostic_dynamic_viscosity(self, T: float) -> float: + """Return the AMESim NASA-table viscosity used by pipe diagnostics. + + pn2pipefr reports Reynolds number with sagum viscosity. Keep this + separate from dynamic_viscosity so matching that diagnostic cannot + alter the already-validated pipe flow or friction dynamics. + """ + + if T <= 0.0: + raise ValueError("Temperature must be positive.") + a, b, c, d = self.nasa_viscosity_coefficients + return 1.0e-7 * exp(a * log(T) + b / T + c / (T * T) + d) + @profile_property("density") @cache_property_calculation("density") def density(self, p: float, T: float) -> float: diff --git a/app/simulation/core/medium.py b/app/simulation/core/medium.py index e394fe8..fd6da1f 100644 --- a/app/simulation/core/medium.py +++ b/app/simulation/core/medium.py @@ -81,6 +81,8 @@ class GasMedium(Protocol): def dynamic_viscosity(self, T: float) -> float: ... + def diagnostic_dynamic_viscosity(self, T: float) -> float: ... + def specific_internal_energy(self, T: float) -> float: ... def specific_internal_energy_at_pressure(self, p: float, T: float) -> float: ... @@ -182,6 +184,16 @@ class IdealGasMedium: / (T + self.sutherland_constant) ) + def diagnostic_dynamic_viscosity(self, T: float) -> float: + """Return the viscosity convention used by derived diagnostics. + + Most media use the same transport property for dynamics and reported + diagnostics. Reference-library media may override this without + changing a calibrated constitutive flow relation. + """ + + return self.dynamic_viscosity(T) + @profile_property("specific_internal_energy") def specific_internal_energy(self, T: float) -> float: delta_T = T - self.T_ref diff --git a/tests/test_amesim_helium_medium.py b/tests/test_amesim_helium_medium.py index 0778e36..d4054e8 100644 --- a/tests/test_amesim_helium_medium.py +++ b/tests/test_amesim_helium_medium.py @@ -72,8 +72,33 @@ class AmesimHeliumPengRobinsonMediumTests(unittest.TestCase): temperature, ) self.assertAlmostEqual(medium.dynamic_viscosity(293.15), 1.96e-5) + self.assertAlmostEqual( + medium.diagnostic_dynamic_viscosity(293.15), + 1.9616132203938165e-5, + ) self.assertEqual(medium.sutherland_constant, 79.4) + def test_pipe_diagnostic_viscosity_reproduces_amesim_nasa_table(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + + self.assertEqual( + medium.nasa_viscosity_coefficients, + (0.7501594, 35.76324, -2212.129, 0.9212635), + ) + self.assertAlmostEqual( + medium.diagnostic_dynamic_viscosity(233.92077861710192), + 1.683135986478913e-5, + delta=1.0e-19, + ) + self.assertNotAlmostEqual( + medium.diagnostic_dynamic_viscosity(233.92077861710192), + medium.dynamic_viscosity(233.92077861710192), + delta=1.0e-8, + ) + + with self.assertRaises(ValueError): + medium.diagnostic_dynamic_viscosity(0.0) + def test_pressure_transport_enthalpy_includes_peng_robinson_departure( self, ) -> None: diff --git a/tests/test_amesim_pnl0002_pnl0003_component.py b/tests/test_amesim_pnl0002_pnl0003_component.py index e470f1b..8e4ef1e 100644 --- a/tests/test_amesim_pnl0002_pnl0003_component.py +++ b/tests/test_amesim_pnl0002_pnl0003_component.py @@ -1,9 +1,12 @@ from __future__ import annotations import unittest -from math import log10, sqrt +from math import log10, pi, sqrt from app.simulation.components.amesim.flow.pipes import AmesimPnl0002, AmesimPnl0003 +from app.simulation.components.amesim.media.mediums import ( + AmesimHeliumPengRobinsonMedium, +) from app.simulation.core.medium import IdealGasMedium from app.simulation.registry import COMPONENT_MODEL_REGISTRY @@ -462,6 +465,33 @@ class AmesimPnl0003ComponentTests(unittest.TestCase): self.assertIn("dmctr", values) self.assertIn("re", values) + def test_reported_reynolds_uses_amesim_helium_nasa_viscosity_only(self) -> None: + medium = AmesimHeliumPengRobinsonMedium() + pipe = AmesimPnl0003("pnl_3", medium, diam=0.02) + mass_flow = 0.04555349965141288 + upstream_temperature = 233.92077861710192 + + # Frozen AMESim dense row: t=0.81 s, branch 2, D=0.02 m. + expected_amesim_reynolds = 172298.96343252173 + reported = pipe.reported_reynolds_number( + mass_flow, + upstream_temperature, + ) + dynamic = pipe.reynolds_number(mass_flow, upstream_temperature) + + self.assertAlmostEqual(reported, expected_amesim_reynolds, delta=1.0e-9) + self.assertNotAlmostEqual(reported, dynamic, delta=1.0) + self.assertAlmostEqual( + dynamic, + 4.0 + * abs(mass_flow) + / ( + pi + * pipe.diam + * medium.dynamic_viscosity(upstream_temperature) + ), + ) + def test_result_cm_excludes_center_resistance_flow_coefficient(self) -> None: pipe = AmesimPnl0003( "pnl_3",