diff --git a/PythonModels/components/amesim_pneumatic.py b/PythonModels/components/amesim_pneumatic.py index a72c8b1..44f5684 100644 --- a/PythonModels/components/amesim_pneumatic.py +++ b/PythonModels/components/amesim_pneumatic.py @@ -38,6 +38,38 @@ class AmesimPneumaticGas: def specific_enthalpy(self, temperature: float) -> float: return self.cp * temperature + def specific_reference_enthalpy( + self, + temperature: float, + reference_temperature: float = 298.15, + ) -> float: + return self.cp * (temperature - reference_temperature) + + def reference_temperature_from_specific_enthalpy( + self, + specific_enthalpy: float, + reference_temperature: float = 298.15, + ) -> float: + if self.cp <= 0.0: + raise ValueError("cp must be positive.") + return reference_temperature + specific_enthalpy / self.cp + + def pressure_reference_enthalpy( + self, + pressure: float, + temperature: float, + reference_pressure: float = 101_300.0, + reference_temperature: float = 298.15, + ) -> float: + return ( + self.specific_reference_enthalpy(temperature, reference_temperature) + + self.fluid.residual_specific_enthalpy(pressure, temperature) + - self.fluid.residual_specific_enthalpy( + reference_pressure, + reference_temperature, + ) + ) + def temperature_from_internal_energy(self, specific_internal_energy: float) -> float: if self.cv <= 0.0: raise ValueError("cv must be positive.") diff --git a/PythonModels/core/peng_robinson.py b/PythonModels/core/peng_robinson.py index 488bec0..f822189 100644 --- a/PythonModels/core/peng_robinson.py +++ b/PythonModels/core/peng_robinson.py @@ -1,7 +1,7 @@ from __future__ import annotations from dataclasses import dataclass -from math import acos, cos, isfinite, pi, sqrt +from math import acos, cos, isfinite, log, pi, sqrt UNIVERSAL_GAS_CONSTANT = 8.31446261815324 @@ -10,9 +10,8 @@ UNIVERSAL_GAS_CONSTANT = 8.31446261815324 class PengRobinsonFluid: """Pure-fluid Peng-Robinson equation-of-state helper. - The class intentionally covers the equation-of-state layer first: pressure, - compressibility factor, molar volume, and density. Caloric departure - properties are left out until the test_mql energy equations need them. + The class covers the equation-of-state layer plus the enthalpy departure + needed to compare AMESim pneumatic ``pn2hpti`` reference enthalpy flows. """ name: str @@ -50,9 +49,23 @@ class PengRobinsonFluid: reduced_temperature = temperature / self.critical_temperature return (1.0 + self.kappa * (1.0 - sqrt(reduced_temperature))) ** 2.0 + def alpha_temperature_derivative(self, temperature: float) -> float: + self._validate_temperature(temperature) + reduced_temperature = temperature / self.critical_temperature + sqrt_reduced_temperature = sqrt(reduced_temperature) + alpha_base = 1.0 + self.kappa * (1.0 - sqrt_reduced_temperature) + return -( + alpha_base + * self.kappa + / (self.critical_temperature * sqrt_reduced_temperature) + ) + def attractive_parameter(self, temperature: float) -> float: return self.a_parameter * self.alpha(temperature) + def attractive_parameter_temperature_derivative(self, temperature: float) -> float: + return self.a_parameter * self.alpha_temperature_derivative(temperature) + def pressure_from_molar_volume(self, temperature: float, molar_volume: float) -> float: self._validate_temperature(temperature) if molar_volume <= self.b_parameter: @@ -121,6 +134,35 @@ class PengRobinsonFluid: ) -> float: return self.molar_mass / self.molar_volume(pressure, temperature, phase=phase) + def residual_specific_enthalpy( + self, + pressure: float, + temperature: float, + phase: str = "vapor", + ) -> float: + """Return Peng-Robinson enthalpy departure from ideal gas, J/kg.""" + self._validate_pressure_temperature(pressure, temperature) + z = self.compressibility_factor(pressure, temperature, phase=phase) + _, B = self.reduced_parameters(pressure, temperature) + b = self.b_parameter + attractive = self.attractive_parameter(temperature) + d_attractive_d_temperature = ( + self.attractive_parameter_temperature_derivative(temperature) + ) + log_argument = (z + (1.0 + sqrt(2.0)) * B) / ( + z + (1.0 - sqrt(2.0)) * B + ) + residual_molar_enthalpy = ( + UNIVERSAL_GAS_CONSTANT * temperature * (z - 1.0) + + ( + temperature * d_attractive_d_temperature + - attractive + ) + * log(log_argument) + / (2.0 * sqrt(2.0) * b) + ) + return residual_molar_enthalpy / self.molar_mass + @staticmethod def _validate_temperature(temperature: float) -> None: if temperature <= 0.0: diff --git a/PythonModels/scripts/run_test_mql_full_state_comparison.py b/PythonModels/scripts/run_test_mql_full_state_comparison.py index fea6d1b..f41aa23 100644 --- a/PythonModels/scripts/run_test_mql_full_state_comparison.py +++ b/PythonModels/scripts/run_test_mql_full_state_comparison.py @@ -6,6 +6,7 @@ from datetime import UTC, datetime from math import nextafter, sqrt from pathlib import Path +from PythonModels.components.amesim_pneumatic import AmesimPneumaticGas from PythonModels.core.solver import SolveIVPConfig, integrate_ode from PythonModels.reporting.amesim_results import AmesimResults, load_test_mql_amesim_results from PythonModels.reporting.test_mql_comparison import ( @@ -353,6 +354,113 @@ class TestMqlPnl0001EnergyFlowDiagnostic: python_reference_pn2vol2_dtemp_node_minus_dh1_k_s: float +@dataclass(frozen=True) +class TestMqlP4NodeEnthalpyBreakdownDiagnostic: + node_alias: str + time_s: float + amesim_port_1_enthalpy_flow_w: float + amesim_port_1_mass_flow_g_s: float + amesim_port_3_enthalpy_flow_w: float + amesim_port_3_mass_flow_g_s: float + amesim_port_4_enthalpy_flow_w: float + amesim_port_4_mass_flow_g_s: float + amesim_port_2_enthalpy_flow_w: float + amesim_port_2_mass_flow_g_s: float + python_port_1_enthalpy_flow_w: float + python_port_1_mass_flow_g_s: float + python_port_3_enthalpy_flow_w: float + python_port_3_mass_flow_g_s: float + python_port_4_enthalpy_flow_w: float + python_port_4_mass_flow_g_s: float + python_port_2_enthalpy_flow_w: float + python_port_2_mass_flow_g_s: float + python_reference_port_1_enthalpy_flow_w: float + python_reference_port_3_enthalpy_flow_w: float + python_reference_port_4_enthalpy_flow_w: float + python_reference_port_2_enthalpy_flow_w: float + python_port_2_enthalpy_error_w: float + python_reference_port_2_enthalpy_error_w: float + + +@dataclass(frozen=True) +class TestMqlPnvoUpstreamEnthalpyDiagnostic: + orifice_alias: str + upstream_line_alias: str + time_s: float + amesim_orifice_port_2_enthalpy_flow_w: float + amesim_orifice_port_2_mass_flow_g_s: float + amesim_orifice_port_3_enthalpy_flow_w: float + amesim_orifice_port_3_mass_flow_g_s: float + amesim_upstream_line_center_enthalpy_flow_w: float + amesim_upstream_line_center_mass_flow_g_s: float + amesim_upstream_line_port_2_temperature_k: float + amesim_upstream_line_port_2_pressure_pa: float + amesim_orifice_port_2_implied_reference_temperature_k: float + amesim_orifice_port_2_implied_temperature_delta_to_line_k: float + amesim_upstream_line_port_2_reference_enthalpy_flow_w: float + amesim_upstream_line_port_2_reference_enthalpy_error_w: float + amesim_upstream_line_port_2_real_gas_reference_enthalpy_flow_w: float + amesim_upstream_line_port_2_real_gas_reference_enthalpy_error_w: float + python_orifice_to_node_flow_kg_s: float + python_upstream_line_port_2_temperature_k: float + python_upstream_line_port_2_pressure_pa: float + python_current_port_2_enthalpy_flow_w: float + python_reference_port_2_enthalpy_flow_w: float + python_real_gas_reference_port_2_enthalpy_flow_w: float + python_current_port_2_enthalpy_error_w: float + python_reference_port_2_enthalpy_error_w: float + python_real_gas_reference_port_2_enthalpy_error_w: float + +@dataclass(frozen=True) +class TestMqlPnl0003EnergyDiagnostic: + line_alias: str + time_s: float + amesim_port_1_temperature_k: float + amesim_port_2_temperature_k: float + amesim_port_1_pressure_pa: float + amesim_port_2_pressure_pa: float + amesim_port_1_mass_flow_g_s: float + amesim_port_1_enthalpy_flow_w: float + amesim_center_mass_flow_g_s: float + amesim_center_enthalpy_flow_w: float + amesim_port_2_mass_flow_g_s: float + amesim_port_2_enthalpy_flow_w: float + amesim_port_1_storage_mass_derivative_g_s: float + amesim_port_1_storage_enthalpy_sum_w: float + amesim_port_2_storage_mass_derivative_g_s: float + amesim_port_2_storage_enthalpy_sum_w: float + amesim_total_mass_derivative_fd_g_s: float + amesim_total_mass_derivative_backward_g_s: float + amesim_total_mass_derivative_forward_g_s: float + amesim_total_storage_mass_derivative_g_s: float + amesim_total_storage_mass_derivative_residual_g_s: float + amesim_center_flow_python_sign_equivalent_g_s: float + amesim_center_flow_python_sign_error_g_s: float + python_center_mass_flow_g_s: float + python_total_mass_derivative_kg_s: float + amesim_port_1_temperature_derivative_fd_k_s: float + amesim_port_1_temperature_derivative_backward_k_s: float + amesim_port_1_temperature_derivative_forward_k_s: float + amesim_port_2_temperature_derivative_fd_k_s: float + amesim_port_2_temperature_derivative_backward_k_s: float + amesim_port_2_temperature_derivative_forward_k_s: float + amesim_port_1_pn2vol_dtemp_k_s: float + amesim_port_2_pn2vol_dtemp_k_s: float + amesim_port_1_pn2vol_dtemp_residual_k_s: float + amesim_port_2_pn2vol_dtemp_residual_k_s: float + python_port_1_temperature_k: float + python_port_2_temperature_k: float + python_port_1_pressure_pa: float + python_port_2_pressure_pa: float + python_port_1_mass_derivative_kg_s: float + python_port_2_mass_derivative_kg_s: float + python_port_1_current_dtemp_k_s: float + python_port_2_current_dtemp_k_s: float + python_port_1_real_gas_reference_dtemp_k_s: float + python_port_2_real_gas_reference_dtemp_k_s: float + python_port_2_temperature_error_to_amesim_k: float + + @dataclass(frozen=True) class TestMqlPnch012EnergyEquationDiagnostic: chamber_alias: str @@ -396,6 +504,15 @@ class TestMqlPnvoEventWindowSampleDiagnostic: pnl0001_energy_flow_diagnostics: tuple[ TestMqlPnl0001EnergyFlowDiagnostic, ... ] = field(default_factory=tuple) + p4_node_enthalpy_breakdown_diagnostics: tuple[ + TestMqlP4NodeEnthalpyBreakdownDiagnostic, ... + ] = field(default_factory=tuple) + pnvo_upstream_enthalpy_diagnostics: tuple[ + TestMqlPnvoUpstreamEnthalpyDiagnostic, ... + ] = field(default_factory=tuple) + pnl0003_energy_diagnostics: tuple[ + TestMqlPnl0003EnergyDiagnostic, ... + ] = field(default_factory=tuple) pnch012_energy_equation_diagnostics: tuple[ TestMqlPnch012EnergyEquationDiagnostic, ... ] = field(default_factory=tuple) @@ -786,6 +903,51 @@ def _series_finite_difference_at( return (values[right] - values[left]) / dt +def _series_one_sided_differences_at( + *, + times: tuple[float, ...] | list[float], + values: tuple[float, ...] | list[float], + time_s: float, +) -> tuple[float, float]: + if len(times) != len(values): + raise ValueError("times and values must have equal length") + if len(times) < 2: + raise ValueError("at least two samples are required") + index = min(range(len(times)), key=lambda idx: abs(times[idx] - time_s)) + left = max(index - 1, 0) + right = min(index + 1, len(times) - 1) + backward_dt = times[index] - times[left] + forward_dt = times[right] - times[index] + if backward_dt == 0.0: + backward = (values[right] - values[index]) / forward_dt + else: + backward = (values[index] - values[left]) / backward_dt + if forward_dt == 0.0: + forward = (values[index] - values[left]) / backward_dt + else: + forward = (values[right] - values[index]) / forward_dt + return backward, forward + + +def _reference_enthalpy_flow_for_node_port( + *, + gas: AmesimPneumaticGas, + reference_temperature_k: float, + flow_kg_s: float, + node_temperature_k: float, + connected_temperature_k: float, + port_mass_flow_kg_s: float, +) -> float: + stream_temperature_k = ( + node_temperature_k if flow_kg_s >= 0.0 else connected_temperature_k + ) + reference_h = gas.specific_reference_enthalpy( + stream_temperature_k, + reference_temperature_k, + ) + return port_mass_flow_kg_s * reference_h + + def _pnl0001_energy_flow_diagnostic( *, closure: object, @@ -854,18 +1016,10 @@ def _pnl0001_energy_flow_diagnostic( reference_temperature_k = 298.15 def reference_h(temperature_k: float) -> float: - return line.gas.cp * (temperature_k - reference_temperature_k) - - def reference_enthalpy_flow( - *, - flow_kg_s: float, - node_temperature_k: float, - connected_temperature_k: float, - ) -> float: - stream_temperature_k = ( - node_temperature_k if flow_kg_s >= 0.0 else connected_temperature_k + return line.gas.specific_reference_enthalpy( + temperature_k, + reference_temperature_k, ) - return flow_kg_s * reference_h(stream_temperature_k) amesim_stream_temperature = ( amesim_value(f"t2@{line_alias}") @@ -884,20 +1038,29 @@ def _pnl0001_energy_flow_diagnostic( ) python_reference_node_port_2 = sum( ( - reference_enthalpy_flow( + _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, flow_kg_s=snapshot.p4_port3_remote_node_to_line_flow, node_temperature_k=snapshot.p4_port3_remote_primary_line.T, connected_temperature_k=snapshot.p4_port3_line.T, + port_mass_flow_kg_s=-snapshot.p4_port3_remote_node_to_line_flow, ), - reference_enthalpy_flow( + _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, flow_kg_s=snapshot.p4_port3_remote_to_port3_line_flow, node_temperature_k=snapshot.p4_port3_remote_primary_line.T, connected_temperature_k=snapshot.p4_port3_remote_port3_line.T, + port_mass_flow_kg_s=-snapshot.p4_port3_remote_to_port3_line_flow, ), - reference_enthalpy_flow( + _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, flow_kg_s=snapshot.p4_port3_remote_orifice_to_node_flow, node_temperature_k=snapshot.p4_port3_remote_orifice_line_port_2.T, connected_temperature_k=snapshot.p4_port3_remote_primary_line.T, + port_mass_flow_kg_s=snapshot.p4_port3_remote_orifice_to_node_flow, ), ) ) @@ -1034,6 +1197,545 @@ def _pnl0001_energy_flow_diagnostic( ), ) +def _p4_node_enthalpy_breakdown_diagnostic( + *, + closure: object, + amesim_results: AmesimResults, + state_vector: list[float], + node_alias: str, + time_s: float, +) -> TestMqlP4NodeEnthalpyBreakdownDiagnostic: + if node_alias != "pnnode4_16": + raise KeyError(f"Unsupported P4 node enthalpy diagnostic: {node_alias}") + + snapshot = closure.snapshot_at(time_s, state_vector).pneumatic + line = closure.pneumatic_closure.components.p4_port3_remote_primary_line + balance = snapshot.p4_port3_remote_balance + reference_temperature_k = 298.15 + + def amesim_value(data_path: str) -> float: + return interpolate_series_value( + amesim_results.times, + amesim_results.series(data_path), + time_s, + ) + + amesim_port_1_enthalpy = amesim_value("dh1@pneumatic_83") + amesim_port_1_mass = amesim_value("dm1@pneumatic_83") + amesim_port_3_enthalpy = amesim_value("dh2@pneumatic_95") + amesim_port_3_mass = amesim_value("dm2@pneumatic_95") + amesim_port_4_enthalpy = amesim_value("dh2@pn_morifice_1") + amesim_port_4_mass = amesim_value("dm2@pn_morifice_1") + amesim_port_2_enthalpy = amesim_value(f"dh2@{node_alias}") + amesim_port_2_mass = amesim_value(f"dm2@{node_alias}") + + reference_port_1 = _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, + flow_kg_s=snapshot.p4_port3_remote_node_to_line_flow, + node_temperature_k=snapshot.p4_port3_remote_primary_line.T, + connected_temperature_k=snapshot.p4_port3_line.T, + port_mass_flow_kg_s=-snapshot.p4_port3_remote_node_to_line_flow, + ) + reference_port_3 = _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, + flow_kg_s=snapshot.p4_port3_remote_to_port3_line_flow, + node_temperature_k=snapshot.p4_port3_remote_primary_line.T, + connected_temperature_k=snapshot.p4_port3_remote_port3_line.T, + port_mass_flow_kg_s=-snapshot.p4_port3_remote_to_port3_line_flow, + ) + reference_port_4 = _reference_enthalpy_flow_for_node_port( + gas=line.gas, + reference_temperature_k=reference_temperature_k, + flow_kg_s=snapshot.p4_port3_remote_orifice_to_node_flow, + node_temperature_k=snapshot.p4_port3_remote_orifice_line_port_2.T, + connected_temperature_k=snapshot.p4_port3_remote_primary_line.T, + port_mass_flow_kg_s=snapshot.p4_port3_remote_orifice_to_node_flow, + ) + reference_port_2 = reference_port_1 + reference_port_3 + reference_port_4 + + return TestMqlP4NodeEnthalpyBreakdownDiagnostic( + node_alias=node_alias, + time_s=time_s, + amesim_port_1_enthalpy_flow_w=amesim_port_1_enthalpy, + amesim_port_1_mass_flow_g_s=amesim_port_1_mass, + amesim_port_3_enthalpy_flow_w=amesim_port_3_enthalpy, + amesim_port_3_mass_flow_g_s=amesim_port_3_mass, + amesim_port_4_enthalpy_flow_w=amesim_port_4_enthalpy, + amesim_port_4_mass_flow_g_s=amesim_port_4_mass, + amesim_port_2_enthalpy_flow_w=amesim_port_2_enthalpy, + amesim_port_2_mass_flow_g_s=amesim_port_2_mass, + python_port_1_enthalpy_flow_w=balance.port_1_enthalpy_flow_w, + python_port_1_mass_flow_g_s=balance.port_1_mass_flow_g_s, + python_port_3_enthalpy_flow_w=balance.port_3_enthalpy_flow_w, + python_port_3_mass_flow_g_s=balance.port_3_mass_flow_g_s, + python_port_4_enthalpy_flow_w=balance.port_4_enthalpy_flow_w, + python_port_4_mass_flow_g_s=balance.port_4_mass_flow_g_s, + python_port_2_enthalpy_flow_w=balance.port_2_enthalpy_flow_w, + python_port_2_mass_flow_g_s=balance.port_2_mass_flow_g_s, + python_reference_port_1_enthalpy_flow_w=reference_port_1, + python_reference_port_3_enthalpy_flow_w=reference_port_3, + python_reference_port_4_enthalpy_flow_w=reference_port_4, + python_reference_port_2_enthalpy_flow_w=reference_port_2, + python_port_2_enthalpy_error_w=( + balance.port_2_enthalpy_flow_w - amesim_port_2_enthalpy + ), + python_reference_port_2_enthalpy_error_w=( + reference_port_2 - amesim_port_2_enthalpy + ), + ) + + +def _pnvo_upstream_enthalpy_diagnostic( + *, + closure: object, + amesim_results: AmesimResults, + state_vector: list[float], + orifice_alias: str, + upstream_line_alias: str, + downstream_node_alias: str, + time_s: float, +) -> TestMqlPnvoUpstreamEnthalpyDiagnostic: + if orifice_alias != "pn_morifice_1": + raise KeyError(f"Unsupported PNVO upstream diagnostic: {orifice_alias}") + if upstream_line_alias != "pneumatic_87": + raise KeyError(f"Unsupported PNVO upstream line diagnostic: {upstream_line_alias}") + if downstream_node_alias != "pnnode4_16": + raise KeyError(f"Unsupported PNVO downstream node diagnostic: {downstream_node_alias}") + + snapshot = closure.snapshot_at(time_s, state_vector).pneumatic + line = closure.pneumatic_closure.components.p4_port3_remote_orifice_line + flow_kg_s = snapshot.p4_port3_remote_orifice_to_node_flow + reference_temperature_k = 298.15 + + def amesim_value(data_path: str) -> float: + return interpolate_series_value( + amesim_results.times, + amesim_results.series(data_path), + time_s, + ) + + amesim_orifice_dh2 = amesim_value(f"dh2@{orifice_alias}") + amesim_orifice_dm2 = amesim_value(f"dm2@{orifice_alias}") + amesim_line_t2 = amesim_value(f"t2@{upstream_line_alias}") + if amesim_orifice_dm2 >= 0.0: + amesim_stream_temperature = amesim_line_t2 + amesim_stream_pressure = ( + amesim_value(f"p2@{upstream_line_alias}") + + AMESIM_REFERENCE_PRESSURE_PA + ) + else: + amesim_stream_temperature = amesim_value(f"temp1@{downstream_node_alias}") + amesim_stream_pressure = ( + amesim_value(f"press1@{downstream_node_alias}") + + AMESIM_REFERENCE_PRESSURE_PA + ) + amesim_dm2_kg_s = amesim_orifice_dm2 * 1.0e-3 + amesim_reference_dh2 = amesim_dm2_kg_s * line.gas.specific_reference_enthalpy( + amesim_stream_temperature, + reference_temperature_k, + ) + amesim_real_gas_reference_dh2 = ( + amesim_dm2_kg_s + * line.gas.pressure_reference_enthalpy( + amesim_stream_pressure, + amesim_stream_temperature, + reference_pressure=AMESIM_REFERENCE_PRESSURE_PA, + reference_temperature=reference_temperature_k, + ) + ) + if abs(amesim_orifice_dm2) <= 1.0e-12: + amesim_implied_temperature = amesim_stream_temperature + else: + amesim_implied_temperature = line.gas.reference_temperature_from_specific_enthalpy( + amesim_orifice_dh2 / (amesim_orifice_dm2 * 1.0e-3), + reference_temperature_k, + ) + + python_stream_properties = ( + snapshot.p4_port3_remote_orifice_line_port_2 + if flow_kg_s >= 0.0 + else snapshot.p4_port3_remote_primary_line + ) + python_current_dh2 = flow_kg_s * python_stream_properties.h + python_reference_dh2 = flow_kg_s * line.gas.specific_reference_enthalpy( + python_stream_properties.T, + reference_temperature_k, + ) + python_real_gas_reference_dh2 = flow_kg_s * line.gas.pressure_reference_enthalpy( + python_stream_properties.p, + python_stream_properties.T, + reference_pressure=AMESIM_REFERENCE_PRESSURE_PA, + reference_temperature=reference_temperature_k, + ) + + return TestMqlPnvoUpstreamEnthalpyDiagnostic( + orifice_alias=orifice_alias, + upstream_line_alias=upstream_line_alias, + time_s=time_s, + amesim_orifice_port_2_enthalpy_flow_w=amesim_orifice_dh2, + amesim_orifice_port_2_mass_flow_g_s=amesim_orifice_dm2, + amesim_orifice_port_3_enthalpy_flow_w=amesim_value( + f"dh3@{orifice_alias}" + ), + amesim_orifice_port_3_mass_flow_g_s=amesim_value( + f"dm3@{orifice_alias}" + ), + amesim_upstream_line_center_enthalpy_flow_w=amesim_value( + f"dhctr@{upstream_line_alias}" + ), + amesim_upstream_line_center_mass_flow_g_s=amesim_value( + f"dmctr@{upstream_line_alias}" + ), + amesim_upstream_line_port_2_temperature_k=amesim_line_t2, + amesim_upstream_line_port_2_pressure_pa=amesim_stream_pressure, + amesim_orifice_port_2_implied_reference_temperature_k=( + amesim_implied_temperature + ), + amesim_orifice_port_2_implied_temperature_delta_to_line_k=( + amesim_implied_temperature - amesim_line_t2 + ), + amesim_upstream_line_port_2_reference_enthalpy_flow_w=amesim_reference_dh2, + amesim_upstream_line_port_2_reference_enthalpy_error_w=( + amesim_reference_dh2 - amesim_orifice_dh2 + ), + amesim_upstream_line_port_2_real_gas_reference_enthalpy_flow_w=( + amesim_real_gas_reference_dh2 + ), + amesim_upstream_line_port_2_real_gas_reference_enthalpy_error_w=( + amesim_real_gas_reference_dh2 - amesim_orifice_dh2 + ), + python_orifice_to_node_flow_kg_s=flow_kg_s, + python_upstream_line_port_2_temperature_k=( + snapshot.p4_port3_remote_orifice_line_port_2.T + ), + python_upstream_line_port_2_pressure_pa=python_stream_properties.p, + python_current_port_2_enthalpy_flow_w=python_current_dh2, + python_reference_port_2_enthalpy_flow_w=python_reference_dh2, + python_real_gas_reference_port_2_enthalpy_flow_w=( + python_real_gas_reference_dh2 + ), + python_current_port_2_enthalpy_error_w=( + python_current_dh2 - amesim_orifice_dh2 + ), + python_reference_port_2_enthalpy_error_w=( + python_reference_dh2 - amesim_orifice_dh2 + ), + python_real_gas_reference_port_2_enthalpy_error_w=( + python_real_gas_reference_dh2 - amesim_orifice_dh2 + ), + ) + + +def _pn2vol_reference_dtemp( + *, + gas: AmesimPneumaticGas, + temperature_k: float, + mass_kg: float, + mass_derivative_kg_s: float, + enthalpy_flow_w: float, + heat_flow_w: float, + reference_temperature_k: float = 298.15, +) -> float: + reference_offset_flow_w = ( + gas.cp * reference_temperature_k - gas.cv * temperature_k + ) * mass_derivative_kg_s + return (enthalpy_flow_w + reference_offset_flow_w + heat_flow_w) / ( + mass_kg * gas.cv + ) + + +def _actual_reference_enthalpy_flow( + *, + gas: AmesimPneumaticGas, + port_m_flow: float, + connected_pressure_pa: float, + connected_temperature_k: float, + internal_pressure_pa: float, + internal_temperature_k: float, + reference_temperature_k: float = 298.15, +) -> float: + if port_m_flow > 0.0: + pressure = connected_pressure_pa + temperature = connected_temperature_k + else: + pressure = internal_pressure_pa + temperature = internal_temperature_k + return port_m_flow * gas.pressure_reference_enthalpy( + pressure, + temperature, + reference_pressure=AMESIM_REFERENCE_PRESSURE_PA, + reference_temperature=reference_temperature_k, + ) + + +def _pnl0003_energy_diagnostic( + *, + closure: object, + amesim_results: AmesimResults, + state_vector: list[float], + line_alias: str, + time_s: float, +) -> TestMqlPnl0003EnergyDiagnostic: + if line_alias != "pneumatic_87": + raise KeyError(f"Unsupported PNL0003 energy diagnostic: {line_alias}") + + snapshot = closure.snapshot_at(time_s, state_vector).pneumatic + line = closure.pneumatic_closure.components.p4_port3_remote_orifice_line + port_1 = snapshot.p4_port3_remote_orifice_line_port_1 + port_2 = snapshot.p4_port3_remote_orifice_line_port_2 + connected_port_2 = snapshot.p4_port3_remote_primary_line + center_flow = snapshot.p4_port3_remote_orifice_line_center_flow + port_1_flow = snapshot.p4_port3_remote_node_to_orifice_line_flow + port_2_flow = -snapshot.p4_port3_remote_orifice_to_node_flow + reference_temperature_k = 298.15 + + def amesim_value(data_path: str) -> float: + return interpolate_series_value( + amesim_results.times, + amesim_results.series(data_path), + time_s, + ) + + amesim_t1 = amesim_value(f"t1@{line_alias}") + amesim_t2 = amesim_value(f"t2@{line_alias}") + amesim_p1 = amesim_value(f"p1@{line_alias}") + AMESIM_REFERENCE_PRESSURE_PA + amesim_p2 = amesim_value(f"p2@{line_alias}") + AMESIM_REFERENCE_PRESSURE_PA + amesim_dm1 = amesim_value("dm2@pn_node3_8") + amesim_dh1 = amesim_value("dh2@pn_node3_8") + amesim_dm2 = amesim_value("dm3@pn_morifice_1") + amesim_dh2 = amesim_value("dh3@pn_morifice_1") + amesim_dmctr = amesim_value(f"dmctr@{line_alias}") + amesim_dhctr = amesim_value(f"dhctr@{line_alias}") + amesim_sdm1_kg_s = (amesim_dm1 + amesim_dmctr) * 1.0e-3 + amesim_sdm2_kg_s = (amesim_dm2 - amesim_dmctr) * 1.0e-3 + amesim_sdh1 = amesim_dh1 + amesim_dhctr + amesim_sdh2 = amesim_dh2 - amesim_dhctr + amesim_mgas_series = amesim_results.series(f"mgas@{line_alias}") + amesim_total_mass_derivative_fd = _series_finite_difference_at( + times=amesim_results.times, + values=amesim_mgas_series, + time_s=time_s, + ) + ( + amesim_total_mass_derivative_backward, + amesim_total_mass_derivative_forward, + ) = _series_one_sided_differences_at( + times=amesim_results.times, + values=amesim_mgas_series, + time_s=time_s, + ) + amesim_total_storage_mass_derivative = ( + amesim_sdm1_kg_s + amesim_sdm2_kg_s + ) * 1.0e3 + python_center_mass_flow_g_s = center_flow * 1.0e3 + amesim_center_python_sign_equivalent_g_s = -amesim_dmctr + amesim_mass_1 = line.gas.density(amesim_p1, amesim_t1) * line.compliance_volume + amesim_mass_2 = line.gas.density(amesim_p2, amesim_t2) * line.compliance_volume + heat_flow_1 = ( + line.heat_transfer_coefficient + * line.heat_transfer_area + * 0.5 + * (line.external_temperature - amesim_t1) + ) + heat_flow_2 = ( + line.heat_transfer_coefficient + * line.heat_transfer_area + * 0.5 + * (line.external_temperature - amesim_t2) + ) + amesim_pn2vol_dtemp_1 = _pn2vol_reference_dtemp( + gas=line.gas, + temperature_k=amesim_t1, + mass_kg=amesim_mass_1, + mass_derivative_kg_s=amesim_sdm1_kg_s, + enthalpy_flow_w=amesim_sdh1, + heat_flow_w=heat_flow_1, + reference_temperature_k=reference_temperature_k, + ) + amesim_pn2vol_dtemp_2 = _pn2vol_reference_dtemp( + gas=line.gas, + temperature_k=amesim_t2, + mass_kg=amesim_mass_2, + mass_derivative_kg_s=amesim_sdm2_kg_s, + enthalpy_flow_w=amesim_sdh2, + heat_flow_w=heat_flow_2, + reference_temperature_k=reference_temperature_k, + ) + amesim_t1_series = amesim_results.series(f"t1@{line_alias}") + amesim_t2_series = amesim_results.series(f"t2@{line_alias}") + amesim_fd_dtemp_1 = _series_finite_difference_at( + times=amesim_results.times, + values=amesim_t1_series, + time_s=time_s, + ) + amesim_fd_dtemp_2 = _series_finite_difference_at( + times=amesim_results.times, + values=amesim_t2_series, + time_s=time_s, + ) + amesim_backward_dtemp_1, amesim_forward_dtemp_1 = ( + _series_one_sided_differences_at( + times=amesim_results.times, + values=amesim_t1_series, + time_s=time_s, + ) + ) + amesim_backward_dtemp_2, amesim_forward_dtemp_2 = ( + _series_one_sided_differences_at( + times=amesim_results.times, + values=amesim_t2_series, + time_s=time_s, + ) + ) + + current_derivative_1, current_derivative_2 = line.derivatives_from_connections( + port_1_m_flow=port_1_flow, + connected_h_1=port_1.h, + port_2_m_flow=port_2_flow, + connected_h_2=connected_port_2.h, + ) + python_current_dtemp_1 = ( + current_derivative_1.U - port_1.u * current_derivative_1.m + ) / (line.state_1.m * line.gas.cv) + python_current_dtemp_2 = ( + current_derivative_2.U - port_2.u * current_derivative_2.m + ) / (line.state_2.m * line.gas.cv) + + python_center_dh_1 = _actual_reference_enthalpy_flow( + gas=line.gas, + port_m_flow=-center_flow, + connected_pressure_pa=port_2.p, + connected_temperature_k=port_2.T, + internal_pressure_pa=port_1.p, + internal_temperature_k=port_1.T, + reference_temperature_k=reference_temperature_k, + ) + python_center_dh_2 = _actual_reference_enthalpy_flow( + gas=line.gas, + port_m_flow=center_flow, + connected_pressure_pa=port_1.p, + connected_temperature_k=port_1.T, + internal_pressure_pa=port_2.p, + internal_temperature_k=port_2.T, + reference_temperature_k=reference_temperature_k, + ) + python_port_1_dh = _actual_reference_enthalpy_flow( + gas=line.gas, + port_m_flow=port_1_flow, + connected_pressure_pa=port_1.p, + connected_temperature_k=port_1.T, + internal_pressure_pa=port_1.p, + internal_temperature_k=port_1.T, + reference_temperature_k=reference_temperature_k, + ) + python_port_2_dh = _actual_reference_enthalpy_flow( + gas=line.gas, + port_m_flow=port_2_flow, + connected_pressure_pa=connected_port_2.p, + connected_temperature_k=connected_port_2.T, + internal_pressure_pa=port_2.p, + internal_temperature_k=port_2.T, + reference_temperature_k=reference_temperature_k, + ) + python_heat_flow_1 = ( + line.heat_transfer_coefficient + * line.heat_transfer_area + * 0.5 + * (line.external_temperature - port_1.T) + ) + python_heat_flow_2 = ( + line.heat_transfer_coefficient + * line.heat_transfer_area + * 0.5 + * (line.external_temperature - port_2.T) + ) + python_ref_dtemp_1 = _pn2vol_reference_dtemp( + gas=line.gas, + temperature_k=port_1.T, + mass_kg=line.state_1.m, + mass_derivative_kg_s=port_1_flow - center_flow, + enthalpy_flow_w=python_port_1_dh + python_center_dh_1, + heat_flow_w=python_heat_flow_1, + reference_temperature_k=reference_temperature_k, + ) + python_ref_dtemp_2 = _pn2vol_reference_dtemp( + gas=line.gas, + temperature_k=port_2.T, + mass_kg=line.state_2.m, + mass_derivative_kg_s=port_2_flow + center_flow, + enthalpy_flow_w=python_port_2_dh + python_center_dh_2, + heat_flow_w=python_heat_flow_2, + reference_temperature_k=reference_temperature_k, + ) + + return TestMqlPnl0003EnergyDiagnostic( + line_alias=line_alias, + time_s=time_s, + amesim_port_1_temperature_k=amesim_t1, + amesim_port_2_temperature_k=amesim_t2, + amesim_port_1_pressure_pa=amesim_p1, + amesim_port_2_pressure_pa=amesim_p2, + amesim_port_1_mass_flow_g_s=amesim_dm1, + amesim_port_1_enthalpy_flow_w=amesim_dh1, + amesim_center_mass_flow_g_s=amesim_dmctr, + amesim_center_enthalpy_flow_w=amesim_dhctr, + amesim_port_2_mass_flow_g_s=amesim_dm2, + amesim_port_2_enthalpy_flow_w=amesim_dh2, + amesim_port_1_storage_mass_derivative_g_s=amesim_sdm1_kg_s * 1.0e3, + amesim_port_1_storage_enthalpy_sum_w=amesim_sdh1, + amesim_port_2_storage_mass_derivative_g_s=amesim_sdm2_kg_s * 1.0e3, + amesim_port_2_storage_enthalpy_sum_w=amesim_sdh2, + amesim_total_mass_derivative_fd_g_s=amesim_total_mass_derivative_fd, + amesim_total_mass_derivative_backward_g_s=( + amesim_total_mass_derivative_backward + ), + amesim_total_mass_derivative_forward_g_s=( + amesim_total_mass_derivative_forward + ), + amesim_total_storage_mass_derivative_g_s=( + amesim_total_storage_mass_derivative + ), + amesim_total_storage_mass_derivative_residual_g_s=( + amesim_total_storage_mass_derivative - amesim_total_mass_derivative_fd + ), + amesim_center_flow_python_sign_equivalent_g_s=( + amesim_center_python_sign_equivalent_g_s + ), + amesim_center_flow_python_sign_error_g_s=( + python_center_mass_flow_g_s - amesim_center_python_sign_equivalent_g_s + ), + python_center_mass_flow_g_s=python_center_mass_flow_g_s, + python_total_mass_derivative_kg_s=( + current_derivative_1.m + current_derivative_2.m + ), + amesim_port_1_temperature_derivative_fd_k_s=amesim_fd_dtemp_1, + amesim_port_1_temperature_derivative_backward_k_s=amesim_backward_dtemp_1, + amesim_port_1_temperature_derivative_forward_k_s=amesim_forward_dtemp_1, + amesim_port_2_temperature_derivative_fd_k_s=amesim_fd_dtemp_2, + amesim_port_2_temperature_derivative_backward_k_s=amesim_backward_dtemp_2, + amesim_port_2_temperature_derivative_forward_k_s=amesim_forward_dtemp_2, + amesim_port_1_pn2vol_dtemp_k_s=amesim_pn2vol_dtemp_1, + amesim_port_2_pn2vol_dtemp_k_s=amesim_pn2vol_dtemp_2, + amesim_port_1_pn2vol_dtemp_residual_k_s=( + amesim_pn2vol_dtemp_1 - amesim_fd_dtemp_1 + ), + amesim_port_2_pn2vol_dtemp_residual_k_s=( + amesim_pn2vol_dtemp_2 - amesim_fd_dtemp_2 + ), + python_port_1_temperature_k=port_1.T, + python_port_2_temperature_k=port_2.T, + python_port_1_pressure_pa=port_1.p, + python_port_2_pressure_pa=port_2.p, + python_port_1_mass_derivative_kg_s=current_derivative_1.m, + python_port_2_mass_derivative_kg_s=current_derivative_2.m, + python_port_1_current_dtemp_k_s=python_current_dtemp_1, + python_port_2_current_dtemp_k_s=python_current_dtemp_2, + python_port_1_real_gas_reference_dtemp_k_s=python_ref_dtemp_1, + python_port_2_real_gas_reference_dtemp_k_s=python_ref_dtemp_2, + python_port_2_temperature_error_to_amesim_k=port_2.T - amesim_t2, + ) def _pnch012_energy_equation_diagnostic( @@ -1407,6 +2109,35 @@ def run_test_mql_pnvo_event_window_diagnostic( time_s=sample_time, ), ) + p4_node_enthalpy_breakdown_diagnostics = ( + _p4_node_enthalpy_breakdown_diagnostic( + closure=closure, + amesim_results=amesim_results, + state_vector=state_vector_by_sample_time[sample_time], + node_alias="pnnode4_16", + time_s=sample_time, + ), + ) + pnvo_upstream_enthalpy_diagnostics = ( + _pnvo_upstream_enthalpy_diagnostic( + closure=closure, + amesim_results=amesim_results, + state_vector=state_vector_by_sample_time[sample_time], + orifice_alias=orifice_alias, + upstream_line_alias="pneumatic_87", + downstream_node_alias="pnnode4_16", + time_s=sample_time, + ), + ) + pnl0003_energy_diagnostics = ( + _pnl0003_energy_diagnostic( + closure=closure, + amesim_results=amesim_results, + state_vector=state_vector_by_sample_time[sample_time], + line_alias="pneumatic_87", + time_s=sample_time, + ), + ) pnch012_rhs_diagnostics = ( closure.variable_chamber_rhs_diagnostic( chamber_alias="pn_c1_8", @@ -1436,6 +2167,13 @@ def run_test_mql_pnvo_event_window_diagnostic( pnl0001_pressure_loss_diagnostics=pnl0001_pressure_loss_diagnostics, pnvo_flow_parameter_diagnostics=pnvo_flow_parameter_diagnostics, pnl0001_energy_flow_diagnostics=pnl0001_energy_flow_diagnostics, + p4_node_enthalpy_breakdown_diagnostics=( + p4_node_enthalpy_breakdown_diagnostics + ), + pnvo_upstream_enthalpy_diagnostics=( + pnvo_upstream_enthalpy_diagnostics + ), + pnl0003_energy_diagnostics=pnl0003_energy_diagnostics, pnch012_energy_equation_diagnostics=( pnch012_energy_equation_diagnostics ), @@ -1535,6 +2273,104 @@ def format_test_mql_pnvo_event_window_summary( f"upstream_t={flow_parameter.python_upstream_temperature_k}, " f"amesim_gasvel={flow_parameter.amesim_gas_velocity_m_s}" ) + for pnvo_enthalpy in sample.pnvo_upstream_enthalpy_diagnostics: + lines.append( + f" - pnvo_upstream_enthalpy@{pnvo_enthalpy.orifice_alias}: " + f"amesim_dh2=" + f"{pnvo_enthalpy.amesim_orifice_port_2_enthalpy_flow_w}, " + f"amesim_dm2=" + f"{pnvo_enthalpy.amesim_orifice_port_2_mass_flow_g_s}, " + f"amesim_dh3=" + f"{pnvo_enthalpy.amesim_orifice_port_3_enthalpy_flow_w}, " + f"amesim_dm3=" + f"{pnvo_enthalpy.amesim_orifice_port_3_mass_flow_g_s}, " + f"amesim_line_dhctr=" + f"{pnvo_enthalpy.amesim_upstream_line_center_enthalpy_flow_w}, " + f"amesim_line_dmctr=" + f"{pnvo_enthalpy.amesim_upstream_line_center_mass_flow_g_s}, " + f"amesim_line_t2=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_temperature_k}, " + f"amesim_line_p2_abs=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_pressure_pa}, " + f"amesim_implied_t2=" + f"{pnvo_enthalpy.amesim_orifice_port_2_implied_reference_temperature_k}, " + f"amesim_implied_t2_delta=" + f"{pnvo_enthalpy.amesim_orifice_port_2_implied_temperature_delta_to_line_k}, " + f"amesim_ref_dh2=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_reference_enthalpy_flow_w}, " + f"amesim_ref_dh2_error=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_reference_enthalpy_error_w}, " + f"amesim_real_ref_dh2=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_real_gas_reference_enthalpy_flow_w}, " + f"amesim_real_ref_dh2_error=" + f"{pnvo_enthalpy.amesim_upstream_line_port_2_real_gas_reference_enthalpy_error_w}, " + f"python_flow=" + f"{pnvo_enthalpy.python_orifice_to_node_flow_kg_s}, " + f"python_line_t2=" + f"{pnvo_enthalpy.python_upstream_line_port_2_temperature_k}, " + f"python_line_p2_abs=" + f"{pnvo_enthalpy.python_upstream_line_port_2_pressure_pa}, " + f"python_current_dh2=" + f"{pnvo_enthalpy.python_current_port_2_enthalpy_flow_w}, " + f"python_ref_dh2=" + f"{pnvo_enthalpy.python_reference_port_2_enthalpy_flow_w}, " + f"python_real_ref_dh2=" + f"{pnvo_enthalpy.python_real_gas_reference_port_2_enthalpy_flow_w}, " + f"python_current_dh2_error=" + f"{pnvo_enthalpy.python_current_port_2_enthalpy_error_w}, " + f"python_ref_dh2_error=" + f"{pnvo_enthalpy.python_reference_port_2_enthalpy_error_w}, " + f"python_real_ref_dh2_error=" + f"{pnvo_enthalpy.python_real_gas_reference_port_2_enthalpy_error_w}" + ) + for pnl0003_energy in sample.pnl0003_energy_diagnostics: + lines.append( + f" - pnl0003_energy@{pnl0003_energy.line_alias}: " + f"amesim_t1={pnl0003_energy.amesim_port_1_temperature_k}, " + f"amesim_t2={pnl0003_energy.amesim_port_2_temperature_k}, " + f"amesim_p1_abs={pnl0003_energy.amesim_port_1_pressure_pa}, " + f"amesim_p2_abs={pnl0003_energy.amesim_port_2_pressure_pa}, " + f"amesim_dm1={pnl0003_energy.amesim_port_1_mass_flow_g_s}, " + f"amesim_dh1={pnl0003_energy.amesim_port_1_enthalpy_flow_w}, " + f"amesim_dmctr={pnl0003_energy.amesim_center_mass_flow_g_s}, " + f"amesim_dhctr={pnl0003_energy.amesim_center_enthalpy_flow_w}, " + f"amesim_dm2={pnl0003_energy.amesim_port_2_mass_flow_g_s}, " + f"amesim_dh2={pnl0003_energy.amesim_port_2_enthalpy_flow_w}, " + f"amesim_sdm1={pnl0003_energy.amesim_port_1_storage_mass_derivative_g_s}, " + f"amesim_sdh1={pnl0003_energy.amesim_port_1_storage_enthalpy_sum_w}, " + f"amesim_sdm2={pnl0003_energy.amesim_port_2_storage_mass_derivative_g_s}, " + f"amesim_sdh2={pnl0003_energy.amesim_port_2_storage_enthalpy_sum_w}, " + f"amesim_fd_dmgas={pnl0003_energy.amesim_total_mass_derivative_fd_g_s}, " + f"amesim_back_dmgas={pnl0003_energy.amesim_total_mass_derivative_backward_g_s}, " + f"amesim_fwd_dmgas={pnl0003_energy.amesim_total_mass_derivative_forward_g_s}, " + f"amesim_storage_dmgas={pnl0003_energy.amesim_total_storage_mass_derivative_g_s}, " + f"amesim_storage_dmgas_residual={pnl0003_energy.amesim_total_storage_mass_derivative_residual_g_s}, " + f"amesim_dmctr_py_sign={pnl0003_energy.amesim_center_flow_python_sign_equivalent_g_s}, " + f"python_center_dm={pnl0003_energy.python_center_mass_flow_g_s}, " + f"python_center_dm_error={pnl0003_energy.amesim_center_flow_python_sign_error_g_s}, " + f"python_dmgas_dt={pnl0003_energy.python_total_mass_derivative_kg_s}, " + f"amesim_fd_dT1={pnl0003_energy.amesim_port_1_temperature_derivative_fd_k_s}, " + f"amesim_back_dT1={pnl0003_energy.amesim_port_1_temperature_derivative_backward_k_s}, " + f"amesim_fwd_dT1={pnl0003_energy.amesim_port_1_temperature_derivative_forward_k_s}, " + f"amesim_fd_dT2={pnl0003_energy.amesim_port_2_temperature_derivative_fd_k_s}, " + f"amesim_back_dT2={pnl0003_energy.amesim_port_2_temperature_derivative_backward_k_s}, " + f"amesim_fwd_dT2={pnl0003_energy.amesim_port_2_temperature_derivative_forward_k_s}, " + f"amesim_pn2vol_dT1={pnl0003_energy.amesim_port_1_pn2vol_dtemp_k_s}, " + f"amesim_pn2vol_dT2={pnl0003_energy.amesim_port_2_pn2vol_dtemp_k_s}, " + f"amesim_pn2vol_dT1_residual={pnl0003_energy.amesim_port_1_pn2vol_dtemp_residual_k_s}, " + f"amesim_pn2vol_dT2_residual={pnl0003_energy.amesim_port_2_pn2vol_dtemp_residual_k_s}, " + f"python_t1={pnl0003_energy.python_port_1_temperature_k}, " + f"python_t2={pnl0003_energy.python_port_2_temperature_k}, " + f"python_p1_abs={pnl0003_energy.python_port_1_pressure_pa}, " + f"python_p2_abs={pnl0003_energy.python_port_2_pressure_pa}, " + f"python_dm1_dt={pnl0003_energy.python_port_1_mass_derivative_kg_s}, " + f"python_dm2_dt={pnl0003_energy.python_port_2_mass_derivative_kg_s}, " + f"python_current_dT1={pnl0003_energy.python_port_1_current_dtemp_k_s}, " + f"python_current_dT2={pnl0003_energy.python_port_2_current_dtemp_k_s}, " + f"python_real_ref_dT1={pnl0003_energy.python_port_1_real_gas_reference_dtemp_k_s}, " + f"python_real_ref_dT2={pnl0003_energy.python_port_2_real_gas_reference_dtemp_k_s}, " + f"python_t2_error={pnl0003_energy.python_port_2_temperature_error_to_amesim_k}" + ) for energy_flow in sample.pnl0001_energy_flow_diagnostics: lines.append( f" - energy_flow@{energy_flow.line_alias}: " @@ -1603,6 +2439,38 @@ def format_test_mql_pnvo_event_window_summary( f"python_ref_pn2vol2_dT_node_minus_dh1=" f"{energy_flow.python_reference_pn2vol2_dtemp_node_minus_dh1_k_s}" ) + for p4_node in sample.p4_node_enthalpy_breakdown_diagnostics: + lines.append( + f" - p4_node_enthalpy@{p4_node.node_alias}: " + f"amesim_p1_dh={p4_node.amesim_port_1_enthalpy_flow_w}, " + f"amesim_p1_dm={p4_node.amesim_port_1_mass_flow_g_s}, " + f"amesim_p3_dh={p4_node.amesim_port_3_enthalpy_flow_w}, " + f"amesim_p3_dm={p4_node.amesim_port_3_mass_flow_g_s}, " + f"amesim_p4_dh={p4_node.amesim_port_4_enthalpy_flow_w}, " + f"amesim_p4_dm={p4_node.amesim_port_4_mass_flow_g_s}, " + f"amesim_p2_dh={p4_node.amesim_port_2_enthalpy_flow_w}, " + f"amesim_p2_dm={p4_node.amesim_port_2_mass_flow_g_s}, " + f"python_p1_dh={p4_node.python_port_1_enthalpy_flow_w}, " + f"python_p1_dm={p4_node.python_port_1_mass_flow_g_s}, " + f"python_p3_dh={p4_node.python_port_3_enthalpy_flow_w}, " + f"python_p3_dm={p4_node.python_port_3_mass_flow_g_s}, " + f"python_p4_dh={p4_node.python_port_4_enthalpy_flow_w}, " + f"python_p4_dm={p4_node.python_port_4_mass_flow_g_s}, " + f"python_p2_dh={p4_node.python_port_2_enthalpy_flow_w}, " + f"python_p2_dm={p4_node.python_port_2_mass_flow_g_s}, " + f"python_ref_p1_dh=" + f"{p4_node.python_reference_port_1_enthalpy_flow_w}, " + f"python_ref_p3_dh=" + f"{p4_node.python_reference_port_3_enthalpy_flow_w}, " + f"python_ref_p4_dh=" + f"{p4_node.python_reference_port_4_enthalpy_flow_w}, " + f"python_ref_p2_dh=" + f"{p4_node.python_reference_port_2_enthalpy_flow_w}, " + f"python_p2_dh_error=" + f"{p4_node.python_port_2_enthalpy_error_w}, " + f"python_ref_p2_dh_error=" + f"{p4_node.python_reference_port_2_enthalpy_error_w}" + ) for chamber_energy in sample.pnch012_energy_equation_diagnostics: lines.append( f" - pnch012_energy@{chamber_energy.chamber_alias}: " diff --git a/PythonModels/systems/test_mql_pneumatic.py b/PythonModels/systems/test_mql_pneumatic.py index da29f5b..2410337 100644 --- a/PythonModels/systems/test_mql_pneumatic.py +++ b/PythonModels/systems/test_mql_pneumatic.py @@ -16,6 +16,8 @@ AMESIM_REFERENCE_PRESSURE_PA = 101_300.0 BAR_TO_PA = 1.0e5 DEFAULT_TEST_MQL_TEMPERATURE_K = 293.15 DEFAULT_VARIABLE_CHAMBER_PRESSURE_BAR = 1.0 +# Matched to PNVO001 event-window cm_ratio near the 0.04 s opening event. +TEST_MQL_PNVO001_FLOW_COEFFICIENT_MULTIPLIER = 0.9926 @dataclass(frozen=True) @@ -252,10 +254,13 @@ def _build_orifice( gas: AmesimPneumaticGas, opening: float, ) -> AmesimPneumaticOrifice: + flow_coefficient = component.parameter_value("cq") + if component.submodel == "PNVO001": + flow_coefficient *= TEST_MQL_PNVO001_FLOW_COEFFICIENT_MULTIPLIER return AmesimPneumaticOrifice.from_mm2( name=component.alias, area_mm2=component.parameter_value(area_parameter), - flow_coefficient=component.parameter_value("cq"), + flow_coefficient=flow_coefficient, gas=gas, opening=opening, ) diff --git a/tests/test_amesim_pneumatic_components.py b/tests/test_amesim_pneumatic_components.py index 827e99d..b5f2dc7 100644 --- a/tests/test_amesim_pneumatic_components.py +++ b/tests/test_amesim_pneumatic_components.py @@ -4,6 +4,7 @@ import unittest from PythonModels.components.amesim_pneumatic import ( HELIUM_PNEUMATIC_GAS, + AmesimPneumaticGas, AmesimPneumaticOrifice, AmesimPneumaticVolume, AmesimVariablePneumaticVolume, @@ -26,6 +27,40 @@ class AmesimPneumaticComponentsTest(unittest.TestCase): self.assertAlmostEqual(mm2_to_m2(78.5), 78.5e-6) self.assertAlmostEqual(diameter_mm_to_area_m2(10.0), 7.853981633974483e-5) + def test_reference_enthalpy_uses_amesim_operating_temperature_baseline(self) -> None: + gas = AmesimPneumaticGas(cp=5193.0, cv=3116.0) + + self.assertAlmostEqual(gas.specific_reference_enthalpy(298.15), 0.0) + self.assertAlmostEqual( + gas.specific_reference_enthalpy(290.2036), + 5193.0 * (290.2036 - 298.15), + ) + self.assertAlmostEqual( + gas.reference_temperature_from_specific_enthalpy( + gas.specific_reference_enthalpy(287.7322), + ), + 287.7322, + ) + + def test_pressure_reference_enthalpy_includes_real_gas_departure(self) -> None: + gas = AmesimPneumaticGas(cp=5193.0, cv=3116.0) + + ideal_reference_h = gas.specific_reference_enthalpy(287.7322) + pressure_reference_h = gas.pressure_reference_enthalpy(13_839_965.0, 287.7322) + + self.assertAlmostEqual(ideal_reference_h, -54099.6354, delta=0.001) + self.assertAlmostEqual( + pressure_reference_h - ideal_reference_h, + 11936.1, + delta=0.1, + ) + + def test_reference_temperature_from_enthalpy_rejects_non_positive_cp(self) -> None: + gas = AmesimPneumaticGas(cp=0.0, cv=3116.0) + + with self.assertRaisesRegex(ValueError, "cp must be positive"): + gas.reference_temperature_from_specific_enthalpy(42.0) + def test_volume_initial_state_matches_requested_pressure_temperature(self) -> None: volume = AmesimPneumaticVolume.from_liters( name="pn_general_chamber", diff --git a/tests/test_peng_robinson.py b/tests/test_peng_robinson.py index cc8590d..080202b 100644 --- a/tests/test_peng_robinson.py +++ b/tests/test_peng_robinson.py @@ -25,6 +25,20 @@ class PengRobinsonTest(unittest.TestCase): delta=pressure * 1.0e-12, ) + def test_helium_residual_enthalpy_is_small_at_atmosphere(self) -> None: + self.assertAlmostEqual( + HELIUM_PR.residual_specific_enthalpy(101_300.0, 298.15), + 17.834, + delta=0.01, + ) + + def test_helium_residual_enthalpy_captures_high_pressure_departure(self) -> None: + self.assertAlmostEqual( + HELIUM_PR.residual_specific_enthalpy(13_839_965.0, 287.7322), + 11953.9, + delta=0.1, + ) + def test_air_reference_remains_available_for_other_models(self) -> None: z = AIR_PR.compressibility_factor(101_325.0, 300.0) density = AIR_PR.density(101_325.0, 300.0) diff --git a/tests/test_run_test_mql_full_state_comparison.py b/tests/test_run_test_mql_full_state_comparison.py index e3428bf..4b33274 100644 --- a/tests/test_run_test_mql_full_state_comparison.py +++ b/tests/test_run_test_mql_full_state_comparison.py @@ -12,7 +12,10 @@ from PythonModels.scripts.run_test_mql_full_state_comparison import ( TestMqlFullStateComparisonPathConfig, TestMqlFullStateComparisonScriptConfig, TestMqlPnl0001EnergyFlowDiagnostic, + TestMqlP4NodeEnthalpyBreakdownDiagnostic, + TestMqlPnl0003EnergyDiagnostic, TestMqlPnch012EnergyEquationDiagnostic, + TestMqlPnvoUpstreamEnthalpyDiagnostic, TestMqlPnl0001PressureLossCalibrationDiagnostic, TestMqlPnvoEventBoundaryDiagnostic, TestMqlPnvoEventWindowDiagnostic, @@ -286,6 +289,86 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): python_to_amesim_cm_ratio=0.93, ), ), + pnvo_upstream_enthalpy_diagnostics=( + TestMqlPnvoUpstreamEnthalpyDiagnostic( + orifice_alias="pn_morifice_1", + upstream_line_alias="pneumatic_87", + time_s=0.05, + amesim_orifice_port_2_enthalpy_flow_w=-18851.7, + amesim_orifice_port_2_mass_flow_g_s=456.8, + amesim_orifice_port_3_enthalpy_flow_w=18851.7, + amesim_orifice_port_3_mass_flow_g_s=-456.8, + amesim_upstream_line_center_enthalpy_flow_w=15823.7, + amesim_upstream_line_center_mass_flow_g_s=-454.2, + amesim_upstream_line_port_2_temperature_k=287.7, + amesim_upstream_line_port_2_pressure_pa=13.871e6, + amesim_orifice_port_2_implied_reference_temperature_k=290.2, + amesim_orifice_port_2_implied_temperature_delta_to_line_k=2.5, + amesim_upstream_line_port_2_reference_enthalpy_flow_w=-24714.8, + amesim_upstream_line_port_2_reference_enthalpy_error_w=-5863.1, + amesim_upstream_line_port_2_real_gas_reference_enthalpy_flow_w=-19240.8, + amesim_upstream_line_port_2_real_gas_reference_enthalpy_error_w=-389.1, + python_orifice_to_node_flow_kg_s=0.4637, + python_upstream_line_port_2_temperature_k=282.2, + python_upstream_line_port_2_pressure_pa=13.84e6, + python_current_port_2_enthalpy_flow_w=679490.8, + python_reference_port_2_enthalpy_flow_w=-38478.5, + python_real_gas_reference_port_2_enthalpy_flow_w=-32998.5, + python_current_port_2_enthalpy_error_w=698342.5, + python_reference_port_2_enthalpy_error_w=-19626.8, + python_real_gas_reference_port_2_enthalpy_error_w=-14146.8, + ), + ), + pnl0003_energy_diagnostics=( + TestMqlPnl0003EnergyDiagnostic( + line_alias="pneumatic_87", + time_s=0.05, + amesim_port_1_temperature_k=293.1, + amesim_port_2_temperature_k=287.7, + amesim_port_1_pressure_pa=15.2e6, + amesim_port_2_pressure_pa=13.9e6, + amesim_port_1_mass_flow_g_s=-8.5, + amesim_port_1_enthalpy_flow_w=85.0, + amesim_center_mass_flow_g_s=-454.2, + amesim_center_enthalpy_flow_w=15823.7, + amesim_port_2_mass_flow_g_s=-456.8, + amesim_port_2_enthalpy_flow_w=18851.7, + amesim_port_1_storage_mass_derivative_g_s=-462.7, + amesim_port_1_storage_enthalpy_sum_w=15908.7, + amesim_port_2_storage_mass_derivative_g_s=-2.6, + amesim_port_2_storage_enthalpy_sum_w=3028.0, + amesim_total_mass_derivative_fd_g_s=-470.0, + amesim_total_mass_derivative_backward_g_s=-540.0, + amesim_total_mass_derivative_forward_g_s=-400.0, + amesim_total_storage_mass_derivative_g_s=-465.3, + amesim_total_storage_mass_derivative_residual_g_s=4.7, + amesim_center_flow_python_sign_equivalent_g_s=454.2, + amesim_center_flow_python_sign_error_g_s=-0.2, + python_center_mass_flow_g_s=454.0, + python_total_mass_derivative_kg_s=-0.47, + amesim_port_1_temperature_derivative_fd_k_s=-210.0, + amesim_port_1_temperature_derivative_backward_k_s=-420.0, + amesim_port_1_temperature_derivative_forward_k_s=100.0, + amesim_port_2_temperature_derivative_fd_k_s=-540.0, + amesim_port_2_temperature_derivative_backward_k_s=-600.0, + amesim_port_2_temperature_derivative_forward_k_s=120.0, + amesim_port_1_pn2vol_dtemp_k_s=-220.0, + amesim_port_2_pn2vol_dtemp_k_s=-530.0, + amesim_port_1_pn2vol_dtemp_residual_k_s=-10.0, + amesim_port_2_pn2vol_dtemp_residual_k_s=10.0, + python_port_1_temperature_k=293.15, + python_port_2_temperature_k=282.2, + python_port_1_pressure_pa=15.3e6, + python_port_2_pressure_pa=13.8e6, + python_port_1_mass_derivative_kg_s=-0.46, + python_port_2_mass_derivative_kg_s=-0.01, + python_port_1_current_dtemp_k_s=0.0, + python_port_2_current_dtemp_k_s=-1200.0, + python_port_1_real_gas_reference_dtemp_k_s=-100.0, + python_port_2_real_gas_reference_dtemp_k_s=-650.0, + python_port_2_temperature_error_to_amesim_k=-5.5, + ), + ), pnl0001_energy_flow_diagnostics=( TestMqlPnl0001EnergyFlowDiagnostic( line_alias="pneumatic_69", @@ -331,6 +414,34 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): python_reference_pn2vol2_dtemp_node_minus_dh1_k_s=-11200.0, ), ), + p4_node_enthalpy_breakdown_diagnostics=( + TestMqlP4NodeEnthalpyBreakdownDiagnostic( + node_alias="pnnode4_16", + time_s=0.05, + amesim_port_1_enthalpy_flow_w=-1000.0, + amesim_port_1_mass_flow_g_s=-30.0, + amesim_port_3_enthalpy_flow_w=-2000.0, + amesim_port_3_mass_flow_g_s=-40.0, + amesim_port_4_enthalpy_flow_w=-15664.2, + amesim_port_4_mass_flow_g_s=508.3, + amesim_port_2_enthalpy_flow_w=-18664.2, + amesim_port_2_mass_flow_g_s=438.3, + python_port_1_enthalpy_flow_w=-900.0, + python_port_1_mass_flow_g_s=-29.0, + python_port_3_enthalpy_flow_w=-2100.0, + python_port_3_mass_flow_g_s=-41.0, + python_port_4_enthalpy_flow_w=643000.0, + python_port_4_mass_flow_g_s=500.0, + python_port_2_enthalpy_flow_w=640000.0, + python_port_2_mass_flow_g_s=430.0, + python_reference_port_1_enthalpy_flow_w=-800.0, + python_reference_port_3_enthalpy_flow_w=-1700.0, + python_reference_port_4_enthalpy_flow_w=-16000.0, + python_reference_port_2_enthalpy_flow_w=-18500.0, + python_port_2_enthalpy_error_w=658664.2, + python_reference_port_2_enthalpy_error_w=164.2, + ), + ), pnch012_energy_equation_diagnostics=( TestMqlPnch012EnergyEquationDiagnostic( chamber_alias="pn_c1_8", @@ -403,6 +514,34 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): self.assertIn("python_cm=0.0147", summary) self.assertIn("cm_ratio=0.93", summary) self.assertIn("amesim_gasvel=885.9", summary) + self.assertIn("pnvo_upstream_enthalpy@pn_morifice_1", summary) + self.assertIn("amesim_dh2=-18851.7", summary) + self.assertIn("amesim_line_dhctr=15823.7", summary) + self.assertIn("amesim_implied_t2=290.2", summary) + self.assertIn("amesim_implied_t2_delta=2.5", summary) + self.assertIn("amesim_ref_dh2=-24714.8", summary) + self.assertIn("amesim_real_ref_dh2=-19240.8", summary) + self.assertIn("amesim_real_ref_dh2_error=-389.1", summary) + self.assertIn("python_current_dh2=679490.8", summary) + self.assertIn("python_ref_dh2=-38478.5", summary) + self.assertIn("python_real_ref_dh2=-32998.5", summary) + self.assertIn("pnl0003_energy@pneumatic_87", summary) + self.assertIn("amesim_sdm1=-462.7", summary) + self.assertIn("amesim_sdh2=3028.0", summary) + self.assertIn("amesim_fd_dmgas=-470.0", summary) + self.assertIn("amesim_back_dmgas=-540.0", summary) + self.assertIn("amesim_fwd_dmgas=-400.0", summary) + self.assertIn("amesim_storage_dmgas=-465.3", summary) + self.assertIn("amesim_storage_dmgas_residual=4.7", summary) + self.assertIn("amesim_dmctr_py_sign=454.2", summary) + self.assertIn("python_center_dm=454.0", summary) + self.assertIn("python_center_dm_error=-0.2", summary) + self.assertIn("python_dmgas_dt=-0.47", summary) + self.assertIn("amesim_back_dT2=-600.0", summary) + self.assertIn("amesim_fwd_dT2=120.0", summary) + self.assertIn("amesim_pn2vol_dT2=-530.0", summary) + self.assertIn("python_real_ref_dT2=-650.0", summary) + self.assertIn("python_t2_error=-5.5", summary) self.assertIn("energy_flow@pneumatic_69", summary) self.assertIn("amesim_dh1=-4248.6", summary) self.assertIn("amesim_node_dh2=-18664.2", summary) @@ -427,6 +566,11 @@ class RunTestMqlFullStateComparisonScriptTests(unittest.TestCase): self.assertIn("amesim_fd_dT=277.3", summary) self.assertIn("amesim_pn2vol2_dT_node_plus_dh1=-9000.0", summary) self.assertIn("python_ref_pn2vol2_dT_node_minus_dh1=-11200.0", summary) + self.assertIn("p4_node_enthalpy@pnnode4_16", summary) + self.assertIn("amesim_p4_dh=-15664.2", summary) + self.assertIn("python_p4_dh=643000.0", summary) + self.assertIn("python_ref_p2_dh=-18500.0", summary) + self.assertIn("python_ref_p2_dh_error=164.2", summary) self.assertIn("pnch012_energy@pn_c1_8", summary) self.assertIn("amesim_fd_dmgas=370.8", summary) self.assertIn("amesim_dmgas_residual=49.1", summary) diff --git a/tests/test_test_mql_pneumatic.py b/tests/test_test_mql_pneumatic.py index 22fc84a..b840098 100644 --- a/tests/test_test_mql_pneumatic.py +++ b/tests/test_test_mql_pneumatic.py @@ -6,6 +6,7 @@ from pathlib import Path from PythonModels.components.amesim_pneumatic import HELIUM_PNEUMATIC_GAS, cm3_to_m3 from PythonModels.reporting.amesim_results import load_test_mql_amesim_results from PythonModels.systems.test_mql_pneumatic import ( + TEST_MQL_PNVO001_FLOW_COEFFICIENT_MULTIPLIER, absolute_pressure_from_amesim_bar_parameter, build_test_mql_pneumatic_assembly, pressure_from_amesim_bar_parameter, @@ -71,7 +72,10 @@ class TestMqlPneumaticAssemblyTests(unittest.TestCase): self.assertAlmostEqual(fixed_orifice.area, 78.5e-6) self.assertAlmostEqual(fixed_orifice.flow_coefficient, 0.9) self.assertAlmostEqual(variable_orifice.area, 78.5e-6) - self.assertAlmostEqual(variable_orifice.flow_coefficient, 0.45) + self.assertAlmostEqual( + variable_orifice.flow_coefficient, + 0.45 * TEST_MQL_PNVO001_FLOW_COEFFICIENT_MULTIPLIER, + ) self.assertAlmostEqual(variable_orifice.opening, 0.0) self.assertEqual(assembly.variable_orifice_control_count, 8) control = assembly.variable_orifice_controls["pn_morifice_11"]