重构质量模型并补充回热器分段建模

brayton_cycle/__init__.py

@@ -10,6 +10,17 @@ from .components import (
     Turbine,
 )
 from .cycles import BraytonCycle
+from .mass_models import (
+    apply_component_mass,
+    calculate_component_mass,
+    condenser_radiator_mass,
+    has_default_mass_model,
+    heater_reactor_mass,
+    initialize_recuperator_segments,
+    recuperator_pche_mass,
+    RecuperatorPCHEMassModel,
+    turbine_tac_mass,
+)
 from .optimization import (
     optimize_rc_fixed_param,
     optimize_rc_param,
@@ -34,14 +45,23 @@ __all__ = [
     "Condenser",
     "Heater",
     "Recuperator",
+    "RecuperatorPCHEMassModel",
     "Turbine",
+    "apply_component_mass",
+    "calculate_component_mass",
+    "condenser_radiator_mass",
     "evaluate_rc_efficiency",
+    "has_default_mass_model",
+    "heater_reactor_mass",
+    "initialize_recuperator_segments",
     "local_rc_component_performance_sensitivity",
     "local_rc_design_sensitivity",
     "optimize_rc_fixed_param",
     "optimize_rc_param",
     "plot_optimization_landscape",
     "plot_sweep_optimization_results",
+    "recuperator_pche_mass",
     "scan_rc_efficiency",
     "sweep_and_optimize_rc",
+    "turbine_tac_mass",
 ]

brayton_cycle/components.py

@@ -1,13 +1,9 @@
 # -*- coding: utf-8 -*-
 """Component models used by Brayton cycle simulations."""
 
-import math
-
-import ctREFPROP.ctREFPROP as ct
-
 
 class ComponentMassMixin:
-    """Shared mass-calculation interface for cycle components."""
+    """Shared mass-result storage interface for cycle components."""
 
     component_type = "component"
 
@@ -15,43 +11,16 @@ class ComponentMassMixin:
         self.mass = None
         self.mass_variables = None
 
-    def calculate_mass(self, mass_model=None, **kwargs):
-        """Calculate and store component mass.
+    def set_mass_result(self, result):
+        """Store a mass-model result on this component.
 
-        The future mass model can be provided as either a callable accepting
-        ``(component, **kwargs)`` or an object exposing ``calculate_mass`` or
-        ``calculate_<component_type>_mass``.
+        Mass formulas live in ``brayton_cycle.mass_models``. Components only
+        keep the result after their thermodynamic calculator has populated
+        ``variables``.
         """
         if self.variables is None:
-            raise ValueError("Run component calculator before mass calculation")
+            raise ValueError("Run component calculator before storing mass")
 
-        if mass_model is None:
-            result = self._calculate_mass(**kwargs)
-        else:
-            result = self._run_external_mass_model(mass_model, **kwargs)
-
-        return self._store_mass_result(result)
-
-    def mass_calculator(self, mass_model=None, **kwargs):
-        return self.calculate_mass(mass_model=mass_model, **kwargs)
-
-    def _calculate_mass(self, **kwargs):
-        raise NotImplementedError(
-            f"{self.__class__.__name__} mass model is not implemented yet. "
-            "Pass a mass_model or override _calculate_mass()."
-        )
-
-    def _run_external_mass_model(self, mass_model, **kwargs):
-        method_name = f"calculate_{self.component_type}_mass"
-        if hasattr(mass_model, method_name):
-            return getattr(mass_model, method_name)(self, **kwargs)
-        if hasattr(mass_model, "calculate_mass"):
-            return mass_model.calculate_mass(self, **kwargs)
-        if callable(mass_model):
-            return mass_model(self, **kwargs)
-        raise TypeError("mass_model must be callable or expose a supported method")
-
-    def _store_mass_result(self, result):
         if isinstance(result, dict):
             if "mass" in result:
                 mass = result["mass"]
@@ -70,8 +39,17 @@ class ComponentMassMixin:
         self.variables["mass_variables"] = self.mass_variables
         return self.mass
 
-class Compressor():
+    def clear_mass_result(self):
+        self.mass = None
+        self.mass_variables = None
+        if self.variables is not None:
+            self.variables.pop("mass", None)
+            self.variables.pop("mass_variables", None)
+
+class Compressor(ComponentMassMixin):
     """压缩机类"""
+    component_type = "compressor"
+
     def __init__(self, name, eff):
         """
         初始化参数
@@ -87,6 +65,7 @@ class Compressor():
         self.name = name
         self.eff = eff
         self.variables = None
+        self._init_mass_interface()
         
     def calculator(self, p_in, T_in, p_out, property_calculator):
         # 先计算熵值
@@ -165,40 +144,6 @@ class Turbine(ComponentMassMixin):
             'pi': p_in / p_out
             }
 
-    def _calculate_mass(self, Pe=None, A=0.5, **kwargs):
-        """Estimate turbine/TAC mass with the empirical TAC formula.
-
-        Formula from the provided reference:
-            M = A * sqrt(0.5*pi*(30.522*ln(Pe) - 5.7178))
-
-        If Pe is not supplied, the current turbine specific work ``Wt`` is used
-        directly as Pe without unit conversion.
-        """
-        if not 0.4 <= A <= 0.8:
-            raise ValueError("A should be within the recommended range 0.4-0.8")
-
-        if Pe is None:
-            Pe = abs(self.variables['Wt'])
-
-        if Pe <= 0:
-            raise ValueError("Pe must be positive")
-
-        fit_term = 30.522 * math.log(Pe) - 5.7178
-        if fit_term <= 0:
-            raise ValueError(
-                "Pe is outside the valid logarithmic domain for this "
-                "empirical mass formula"
-            )
-
-        mass = A * math.sqrt(0.5 * math.pi * fit_term)
-        return {
-            'mass': mass,
-            'A': A,
-            'Pe': Pe,
-            'fit_term': fit_term,
-            'formula': 'A*sqrt(0.5*pi*(30.522*ln(Pe)-5.7178))',
-        }
-
 class Recuperator(ComponentMassMixin):
     """换热器类"""
     component_type = "recuperator"
@@ -389,52 +334,6 @@ class Heater(ComponentMassMixin):
             'Q_in': Q_input,
             }
 
-    def _calculate_mass(
-        self,
-        P_heat_mwt=None,
-        mass_flow_rate=None,
-        shielding_mass_ton=2.8,
-        include_shielding=True,
-        **kwargs,
-    ):
-        """Estimate reactor and shielding mass for the heater module.
-
-        Reactor empirical formula from the provided reference:
-            M_reactor = 0.2195 * P_heat + 0.09836
-
-        P_heat is the reactor thermal power in MWt, and masses are in tons.
-        The shielding mass is added as a constant 2.8 ton by default.
-
-        If P_heat_mwt is not supplied, it is estimated from Q_in and
-        mass_flow_rate, assuming Q_in is kJ/kg and mass_flow_rate is kg/s:
-            P_heat_mwt = abs(Q_in) * mass_flow_rate / 1000
-        """
-        if P_heat_mwt is None:
-            if mass_flow_rate is None:
-                raise ValueError(
-                    "Heater mass calculation requires P_heat_mwt, or "
-                    "mass_flow_rate to estimate P_heat_mwt from Q_in."
-                )
-            P_heat_mwt = abs(self.variables['Q_in']) * mass_flow_rate / 1000
-
-        if P_heat_mwt <= 0:
-            raise ValueError("P_heat_mwt must be positive")
-        if shielding_mass_ton < 0:
-            raise ValueError("shielding_mass_ton must be non-negative")
-
-        reactor_mass_ton = 0.2195 * P_heat_mwt + 0.09836
-        shielding_mass = shielding_mass_ton if include_shielding else 0.0
-        total_mass_ton = reactor_mass_ton + shielding_mass
-
-        return {
-            'mass': total_mass_ton,
-            'reactor_mass_ton': reactor_mass_ton,
-            'shielding_mass_ton': shielding_mass,
-            'P_heat_mwt': P_heat_mwt,
-            'include_shielding': include_shielding,
-            'formula': '0.2195*P_heat+0.09836+shielding_mass',
-        }
-
 class Condenser(ComponentMassMixin):
     """冷凝器类"""
     component_type = "condenser"
@@ -463,127 +362,6 @@ class Condenser(ComponentMassMixin):
             'Q_out': Q_output,
             }
 
-    def _calculate_mass(
-        self,
-        Qc_kw=None,
-        co2_mass_flow_rate=None,
-        coolant_outlet_T=None,
-        water_inlet_T=None,
-        water_mass_flow_rate=None,
-        water_pressure_kpa=101.325,
-        refprop_path="C:/Program Files (x86)/REFPROP 10.0+/REFPROP",
-        emissivity=0.92,
-        surface_temperature=210.15,
-        area_density=6.75,
-        **kwargs,
-    ):
-        """Estimate radiator mass for the condenser module.
-
-        Radiator heat rejection model:
-            Qc = phi * sigma * A_rad * (T**4 - T0**4)
-            M_rad = kappa * A_rad
-
-        The NaK coolant in the reference is represented here by water. If
-        coolant_outlet_T is not supplied, water outlet temperature is evaluated
-        with REFPROP from water_inlet_T, water_mass_flow_rate, and Qc_kw.
-
-        Args:
-            Qc_kw: Total heat rejection in kW.
-            co2_mass_flow_rate: Optional kg/s. If Qc_kw is not supplied, Qc_kw
-                is estimated as Q_out * co2_mass_flow_rate, assuming Q_out is
-                kJ/kg.
-            coolant_outlet_T: Radiator coolant outlet temperature T in K.
-            water_inlet_T: Water inlet temperature in K for REFPROP calculation.
-            water_mass_flow_rate: Water mass flow rate in kg/s.
-            water_pressure_kpa: Water pressure for REFPROP calculation.
-            refprop_path: REFPROP root path.
-            emissivity: Radiator surface emissivity phi.
-            surface_temperature: Ambient/surface temperature T0 in K. Default is
-                -63 degC for Mars, 210.15 K.
-            area_density: Radiator face density kappa in kg/m2.
-        """
-        if Qc_kw is None:
-            if co2_mass_flow_rate is None:
-                raise ValueError(
-                    "Condenser mass calculation requires Qc_kw, or "
-                    "co2_mass_flow_rate to estimate Qc_kw from Q_out."
-                )
-            Qc_kw = abs(self.variables['Q_out']) * co2_mass_flow_rate
-
-        if Qc_kw <= 0:
-            raise ValueError("Qc_kw must be positive")
-        if emissivity <= 0:
-            raise ValueError("emissivity must be positive")
-        if area_density <= 0:
-            raise ValueError("area_density must be positive")
-
-        if coolant_outlet_T is None:
-            coolant_outlet_T = self._calculate_water_outlet_temperature(
-                Qc_kw=Qc_kw,
-                water_inlet_T=water_inlet_T,
-                water_mass_flow_rate=water_mass_flow_rate,
-                water_pressure_kpa=water_pressure_kpa,
-                refprop_path=refprop_path,
-            )
-
-        temperature_term = coolant_outlet_T**4 - surface_temperature**4
-        if temperature_term <= 0:
-            raise ValueError(
-                "coolant_outlet_T must be higher than surface_temperature for "
-                "radiative heat rejection"
-            )
-
-        stefan_boltzmann = 5.670374419e-8
-        Qc_w = Qc_kw * 1000
-        area_m2 = Qc_w / (emissivity * stefan_boltzmann * temperature_term)
-        mass_kg = area_density * area_m2
-
-        return {
-            'mass': mass_kg,
-            'radiator_area_m2': area_m2,
-            'Qc_kw': Qc_kw,
-            'coolant_outlet_T': coolant_outlet_T,
-            'surface_temperature': surface_temperature,
-            'emissivity': emissivity,
-            'area_density': area_density,
-            'mass_unit': 'kg',
-            'formula': 'M_rad=kappa*Qc/(phi*sigma*(T^4-T0^4))',
-        }
-
-    def _calculate_water_outlet_temperature(
-        self,
-        Qc_kw,
-        water_inlet_T,
-        water_mass_flow_rate,
-        water_pressure_kpa,
-        refprop_path,
-    ):
-        if water_inlet_T is None or water_mass_flow_rate is None:
-            raise ValueError(
-                "Provide coolant_outlet_T directly, or provide water_inlet_T "
-                "and water_mass_flow_rate for REFPROP water calculation."
-            )
-        if water_mass_flow_rate <= 0:
-            raise ValueError("water_mass_flow_rate must be positive")
-
-        water = ct.REFPROPFunctionLibrary(refprop_path)
-        water.SETUPdll(1, 'WATER.FLD', 'HMX.BNC', 'DEF')
-        water.SETUPdll(2, 'SI', 'SI', 'DEF')
-        z = [1.0]
-        mw = water.WMOLdll(z)
-
-        inlet = water.TPFLSHdll(water_inlet_T, water_pressure_kpa, z)
-        if inlet.ierr > 0:
-            raise ValueError(f"REFPROP water inlet calculation error: {inlet.ierr}")
-
-        h_in_mass = inlet.h / mw
-        h_out_mass = h_in_mass + Qc_kw / water_mass_flow_rate
-        outlet = water.PHFLSHdll(water_pressure_kpa, h_out_mass * mw, z)
-        if outlet.ierr > 0:
-            raise ValueError(f"REFPROP water outlet calculation error: {outlet.ierr}")
-
-        return outlet.T
-
 class Concentrator(ComponentMassMixin):
     """汇流组件"""
     component_type = "concentrator"

brayton_cycle/cycles.py

@@ -9,6 +9,7 @@ from .components import (
     Recuperator,
     Turbine,
 )
+from .mass_models import apply_component_mass, has_default_mass_model
 from .properties import CO2PropertyCalculator
 
 
@@ -57,13 +58,31 @@ class BraytonCycle:
             else:
                 yield value
 
-    def calculate_component_masses(self, mass_model=None, **kwargs):
+    def calculate_component_masses(self, mass_model=None, strict=False, **kwargs):
         """Calculate mass for each initialized component."""
         component_masses = {}
         for component in self.iter_components():
-            if not hasattr(component, "calculate_mass"):
+            if not hasattr(component, "set_mass_result"):
                 continue
-            mass = component.calculate_mass(mass_model=mass_model, **kwargs)
+            if mass_model is None and not has_default_mass_model(component):
+                if strict:
+                    raise NotImplementedError(
+                        f"No default mass model for component type "
+                        f"{component.component_type!r}."
+                    )
+                continue
+
+            try:
+                mass = apply_component_mass(
+                    component,
+                    mass_model=mass_model,
+                    **kwargs,
+                )
+            except NotImplementedError:
+                if strict:
+                    raise
+                continue
+
             key = component.variables.get("name", component.name)
             component_masses[key] = {
                 "component_type": component.component_type,
@@ -74,18 +93,23 @@ class BraytonCycle:
         self.component_masses = component_masses
         return self.component_masses
 
-    def cycle_mass_calculator(self, mass_model=None, **kwargs):
+    def cycle_mass_calculator(self, mass_model=None, strict=False, **kwargs):
         """Calculate and return total cycle component mass."""
         component_masses = self.calculate_component_masses(
             mass_model=mass_model,
+            strict=strict,
             **kwargs,
         )
         self.total_mass = sum(item["mass"] for item in component_masses.values())
         return self.total_mass
 
-    def calculate_total_mass(self, mass_model=None, **kwargs):
+    def calculate_total_mass(self, mass_model=None, strict=False, **kwargs):
         """Compatibility alias for cycle_mass_calculator."""
-        return self.cycle_mass_calculator(mass_model=mass_model, **kwargs)
+        return self.cycle_mass_calculator(
+            mass_model=mass_model,
+            strict=strict,
+            **kwargs,
+        )
 
     def SC(self, T_low, T_high, p_low, p_high, param=None):
         """Simple Brayton cycle."""

brayton_cycle/mass_models.py

@@ -0,0 +1,861 @@
+# -*- coding: utf-8 -*-
+"""Mass assessment models for Brayton-cycle components."""
+
+import math
+
+
+def _variables_from(component_or_variables):
+    if isinstance(component_or_variables, dict):
+        return component_or_variables
+    if hasattr(component_or_variables, "variables"):
+        if component_or_variables.variables is None:
+            raise ValueError("Run component calculator before mass calculation")
+        return component_or_variables.variables
+    raise TypeError("Expected component variables dict or component object")
+
+
+def turbine_tac_mass(variables, Pe=None, A=0.5, **kwargs):
+    """Estimate turbine/TAC mass with the empirical TAC formula.
+
+    Formula from the provided reference:
+        M = A * sqrt(0.5*pi*(30.522*ln(Pe) - 5.7178))
+
+    If Pe is not supplied, the turbine specific work ``Wt`` is used directly
+    as Pe without unit conversion.
+    """
+    variables = _variables_from(variables)
+
+    if not 0.4 <= A <= 0.8:
+        raise ValueError("A should be within the recommended range 0.4-0.8")
+
+    if Pe is None:
+        Pe = abs(variables["Wt"])
+
+    if Pe <= 0:
+        raise ValueError("Pe must be positive")
+
+    fit_term = 30.522 * math.log(Pe) - 5.7178
+    if fit_term <= 0:
+        raise ValueError(
+            "Pe is outside the valid logarithmic domain for this empirical "
+            "mass formula"
+        )
+
+    mass = A * math.sqrt(0.5 * math.pi * fit_term)
+    return {
+        "mass": mass,
+        "A": A,
+        "Pe": Pe,
+        "fit_term": fit_term,
+        "model": "turbine_tac_mass",
+        "formula": "A*sqrt(0.5*pi*(30.522*ln(Pe)-5.7178))",
+    }
+
+
+def heater_reactor_mass(
+    variables,
+    P_heat_mwt=None,
+    mass_flow_rate=None,
+    shielding_mass_ton=2.8,
+    include_shielding=True,
+    **kwargs,
+):
+    """Estimate reactor and shielding mass for the heater module.
+
+    Reactor empirical formula:
+        M_reactor = 0.2195 * P_heat + 0.09836
+
+    P_heat is the reactor thermal power in MWt, and masses are in tons. If
+    P_heat_mwt is not supplied, it is estimated from Q_in and mass_flow_rate,
+    assuming Q_in is kJ/kg and mass_flow_rate is kg/s.
+    """
+    variables = _variables_from(variables)
+
+    if P_heat_mwt is None:
+        if mass_flow_rate is None:
+            raise ValueError(
+                "Heater mass calculation requires P_heat_mwt, or "
+                "mass_flow_rate to estimate P_heat_mwt from Q_in."
+            )
+        P_heat_mwt = abs(variables["Q_in"]) * mass_flow_rate / 1000
+
+    if P_heat_mwt <= 0:
+        raise ValueError("P_heat_mwt must be positive")
+    if shielding_mass_ton < 0:
+        raise ValueError("shielding_mass_ton must be non-negative")
+
+    reactor_mass_ton = 0.2195 * P_heat_mwt + 0.09836
+    shielding_mass = shielding_mass_ton if include_shielding else 0.0
+    total_mass_ton = reactor_mass_ton + shielding_mass
+
+    return {
+        "mass": total_mass_ton,
+        "reactor_mass_ton": reactor_mass_ton,
+        "shielding_mass_ton": shielding_mass,
+        "P_heat_mwt": P_heat_mwt,
+        "include_shielding": include_shielding,
+        "mass_unit": "ton",
+        "model": "heater_reactor_mass",
+        "formula": "0.2195*P_heat+0.09836+shielding_mass",
+    }
+
+
+def condenser_radiator_mass(
+    variables,
+    Qc_kw=None,
+    co2_mass_flow_rate=None,
+    coolant_outlet_T=None,
+    water_inlet_T=None,
+    water_mass_flow_rate=None,
+    water_pressure_kpa=101.325,
+    refprop_path="C:/Program Files (x86)/REFPROP 10.0+/REFPROP",
+    emissivity=0.92,
+    surface_temperature=210.15,
+    area_density=6.75,
+    **kwargs,
+):
+    """Estimate radiator mass for the condenser module.
+
+    Radiator heat rejection model:
+        Qc = phi * sigma * A_rad * (T**4 - T0**4)
+        M_rad = kappa * A_rad
+
+    The NaK coolant in the reference is represented here by water. If
+    coolant_outlet_T is not supplied, water outlet temperature is evaluated
+    with REFPROP from water_inlet_T, water_mass_flow_rate, and Qc_kw.
+    """
+    variables = _variables_from(variables)
+
+    if Qc_kw is None:
+        if co2_mass_flow_rate is None:
+            raise ValueError(
+                "Condenser mass calculation requires Qc_kw, or "
+                "co2_mass_flow_rate to estimate Qc_kw from Q_out."
+            )
+        Qc_kw = abs(variables["Q_out"]) * co2_mass_flow_rate
+
+    if Qc_kw <= 0:
+        raise ValueError("Qc_kw must be positive")
+    if emissivity <= 0:
+        raise ValueError("emissivity must be positive")
+    if area_density <= 0:
+        raise ValueError("area_density must be positive")
+
+    if coolant_outlet_T is None:
+        coolant_outlet_T = calculate_water_outlet_temperature(
+            Qc_kw=Qc_kw,
+            water_inlet_T=water_inlet_T,
+            water_mass_flow_rate=water_mass_flow_rate,
+            water_pressure_kpa=water_pressure_kpa,
+            refprop_path=refprop_path,
+        )
+
+    temperature_term = coolant_outlet_T**4 - surface_temperature**4
+    if temperature_term <= 0:
+        raise ValueError(
+            "coolant_outlet_T must be higher than surface_temperature for "
+            "radiative heat rejection"
+        )
+
+    stefan_boltzmann = 5.670374419e-8
+    Qc_w = Qc_kw * 1000
+    area_m2 = Qc_w / (emissivity * stefan_boltzmann * temperature_term)
+    mass_kg = area_density * area_m2
+
+    return {
+        "mass": mass_kg,
+        "radiator_area_m2": area_m2,
+        "Qc_kw": Qc_kw,
+        "coolant_outlet_T": coolant_outlet_T,
+        "surface_temperature": surface_temperature,
+        "emissivity": emissivity,
+        "area_density": area_density,
+        "mass_unit": "kg",
+        "model": "condenser_radiator_mass",
+        "formula": "M_rad=kappa*Qc/(phi*sigma*(T^4-T0^4))",
+    }
+
+
+class RecuperatorPCHEMassModel:
+    """PCHE recuperator mass-model interface and segment initializer."""
+
+    def initialize_segments(
+        self,
+        variables,
+        channel_count,
+        heat_transfer_length_m,
+        mass_flow_rate=None,
+        channel_diameter_m=0.002,
+        property_calculator=None,
+        hot_inlet_state=None,
+        hot_outlet_state=None,
+        cold_inlet_state=None,
+        cold_outlet_state=None,
+    ):
+        """Initialize equal-length PCHE recuperator segments."""
+        variables = _variables_from(variables)
+        hot_inlet_state = self._state_from(
+            variables,
+            "hot_inlet_state",
+            hot_inlet_state,
+        )
+        hot_outlet_state = self._state_from(
+            variables,
+            "hot_outlet_state",
+            hot_outlet_state,
+        )
+        cold_inlet_state = self._state_from(
+            variables,
+            "cold_inlet_state",
+            cold_inlet_state,
+        )
+        cold_outlet_state = self._state_from(
+            variables,
+            "cold_outlet_state",
+            cold_outlet_state,
+        )
+
+        if not isinstance(channel_count, int):
+            raise TypeError("channel_count must be an integer")
+        if channel_count <= 0:
+            raise ValueError("channel_count must be positive")
+        if heat_transfer_length_m <= 0:
+            raise ValueError("heat_transfer_length_m must be positive")
+        if mass_flow_rate is not None and mass_flow_rate <= 0:
+            raise ValueError("mass_flow_rate must be positive")
+
+        segment_length_m = heat_transfer_length_m / channel_count
+        segment_mass_flow_rate = (
+            None if mass_flow_rate is None else mass_flow_rate / channel_count
+        )
+        channel_geometry = self._semicircular_channel_geometry(channel_diameter_m)
+        segments = []
+        for index in range(channel_count):
+            start = index / channel_count
+            end = (index + 1) / channel_count
+            hot_segment_in = self._interpolate_state(
+                hot_inlet_state,
+                hot_outlet_state,
+                start,
+            )
+            hot_segment_out = self._interpolate_state(
+                hot_inlet_state,
+                hot_outlet_state,
+                end,
+            )
+            cold_segment_in = self._interpolate_state(
+                cold_inlet_state,
+                cold_outlet_state,
+                1 - end,
+            )
+            cold_segment_out = self._interpolate_state(
+                cold_inlet_state,
+                cold_outlet_state,
+                1 - start,
+            )
+            hot_side = self._side_flow_result(
+                property_calculator,
+                inlet_state=hot_segment_in,
+                outlet_state=hot_segment_out,
+                mass_flow_rate=segment_mass_flow_rate,
+                flow_area_m2=channel_geometry["flow_area_m2"],
+                length_m=segment_length_m,
+                hydraulic_diameter_m=channel_geometry["hydraulic_diameter_m"],
+            )
+            cold_side = self._side_flow_result(
+                property_calculator,
+                inlet_state=cold_segment_in,
+                outlet_state=cold_segment_out,
+                mass_flow_rate=segment_mass_flow_rate,
+                flow_area_m2=channel_geometry["flow_area_m2"],
+                length_m=segment_length_m,
+                hydraulic_diameter_m=channel_geometry["hydraulic_diameter_m"],
+            )
+            segments.append(
+                {
+                    "index": index,
+                    "hot_inlet_state": hot_segment_in,
+                    "hot_outlet_state": hot_side["outlet_state"],
+                    "cold_inlet_state": cold_segment_in,
+                    "cold_outlet_state": cold_side["outlet_state"],
+                    "hot_mean_temperature": hot_side["mean_temperature"],
+                    "hot_mean_pressure": hot_side["mean_pressure"],
+                    "hot_flow_velocity_m_s": hot_side["flow_velocity_m_s"],
+                    "hot_heat_transfer_coefficient": (
+                        hot_side["heat_transfer_coefficient"]
+                    ),
+                    "cold_mean_temperature": cold_side["mean_temperature"],
+                    "cold_mean_pressure": cold_side["mean_pressure"],
+                    "cold_flow_velocity_m_s": cold_side["flow_velocity_m_s"],
+                    "cold_heat_transfer_coefficient": (
+                        cold_side["heat_transfer_coefficient"]
+                    ),
+                    "geometry": {
+                        "channel_count": channel_count,
+                        "length_m": segment_length_m,
+                        "mass_flow_rate": segment_mass_flow_rate,
+                        **channel_geometry,
+                    },
+                }
+            )
+
+        return {
+            "channel_count": channel_count,
+            "heat_transfer_length_m": heat_transfer_length_m,
+            "mass_flow_rate": mass_flow_rate,
+            **channel_geometry,
+            "segment_count": channel_count,
+            "segment_length_m": segment_length_m,
+            "segment_mass_flow_rate": segment_mass_flow_rate,
+            "segments": segments,
+        }
+
+    def calculate_mass(
+        self,
+        variables,
+        hot_inlet_state=None,
+        hot_outlet_state=None,
+        cold_inlet_state=None,
+        cold_outlet_state=None,
+        Q_exchange=None,
+        mass_flow_rate=None,
+        channel_diameter_m=0.002,
+        property_calculator=None,
+        channel_count=None,
+        heat_transfer_length_m=None,
+        number_of_units=None,
+        channel_length_m=None,
+        total_channel_length_m=None,
+        cross_section_area_m2=None,
+        material_density=8360.0,
+        heat_transfer_area_m2=None,
+        pressure_drop_hot_kpa=None,
+        pressure_drop_cold_kpa=None,
+        segment_count=None,
+        **kwargs,
+    ):
+        """Calculate PCHE recuperator mass from initialized design data."""
+        variables = _variables_from(variables)
+        hot_inlet_state = self._state_from(
+            variables,
+            "hot_inlet_state",
+            hot_inlet_state,
+        )
+        hot_outlet_state = self._state_from(
+            variables,
+            "hot_outlet_state",
+            hot_outlet_state,
+        )
+        cold_inlet_state = self._state_from(
+            variables,
+            "cold_inlet_state",
+            cold_inlet_state,
+        )
+        cold_outlet_state = self._state_from(
+            variables,
+            "cold_outlet_state",
+            cold_outlet_state,
+        )
+
+        if channel_count is None:
+            channel_count = number_of_units
+        if heat_transfer_length_m is None:
+            heat_transfer_length_m = channel_length_m
+        if number_of_units is None:
+            number_of_units = channel_count
+        if channel_length_m is None:
+            channel_length_m = heat_transfer_length_m
+
+        if Q_exchange is None:
+            Q_exchange = variables.get("Q_exchange")
+        if material_density <= 0:
+            raise ValueError("material_density must be positive")
+        if cross_section_area_m2 is not None and cross_section_area_m2 <= 0:
+            raise ValueError("cross_section_area_m2 must be positive")
+        if number_of_units is not None and number_of_units <= 0:
+            raise ValueError("number_of_units must be positive")
+        if channel_length_m is not None and channel_length_m <= 0:
+            raise ValueError("channel_length_m must be positive")
+        if total_channel_length_m is not None and total_channel_length_m <= 0:
+            raise ValueError("total_channel_length_m must be positive")
+        if heat_transfer_area_m2 is not None and heat_transfer_area_m2 <= 0:
+            raise ValueError("heat_transfer_area_m2 must be positive")
+        if mass_flow_rate is not None and mass_flow_rate <= 0:
+            raise ValueError("mass_flow_rate must be positive")
+
+        segment_data = None
+        if channel_count is not None and heat_transfer_length_m is not None:
+            segment_data = self.initialize_segments(
+                variables,
+                channel_count=channel_count,
+                heat_transfer_length_m=heat_transfer_length_m,
+                mass_flow_rate=mass_flow_rate,
+                channel_diameter_m=channel_diameter_m,
+                property_calculator=property_calculator,
+                hot_inlet_state=hot_inlet_state,
+                hot_outlet_state=hot_outlet_state,
+                cold_inlet_state=cold_inlet_state,
+                cold_outlet_state=cold_outlet_state,
+            )
+            segment_count = segment_data["segment_count"]
+
+        if total_channel_length_m is None:
+            if number_of_units is not None and channel_length_m is not None:
+                total_channel_length_m = number_of_units * channel_length_m
+
+        if total_channel_length_m is None or cross_section_area_m2 is None:
+            raise NotImplementedError(
+                "Segmented PCHE recuperator mass calculation is not "
+                "implemented yet. Provide final geometric design results "
+                "(number_of_units + channel_length_m + cross_section_area_m2, "
+                "or total_channel_length_m + cross_section_area_m2) to close "
+                "M_heat=N*L*Area*rho."
+            )
+
+        mass_kg = total_channel_length_m * cross_section_area_m2
+        mass_kg *= material_density
+        return {
+            "mass": mass_kg,
+            "mass_unit": "kg",
+            "model": "recuperator_pche_mass",
+            "formula": "M_heat=N*L*Area*rho",
+            "hot_inlet_state": hot_inlet_state,
+            "hot_outlet_state": hot_outlet_state,
+            "cold_inlet_state": cold_inlet_state,
+            "cold_outlet_state": cold_outlet_state,
+            "Q_exchange": Q_exchange,
+            "mass_flow_rate": mass_flow_rate,
+            **self._semicircular_channel_geometry(channel_diameter_m),
+            "channel_count": channel_count,
+            "heat_transfer_length_m": heat_transfer_length_m,
+            "number_of_units": number_of_units,
+            "channel_length_m": channel_length_m,
+            "total_channel_length_m": total_channel_length_m,
+            "cross_section_area_m2": cross_section_area_m2,
+            "material_density": material_density,
+            "heat_transfer_area_m2": heat_transfer_area_m2,
+            "pressure_drop_hot_kpa": pressure_drop_hot_kpa,
+            "pressure_drop_cold_kpa": pressure_drop_cold_kpa,
+            "segment_count": segment_count,
+            "segment_length_m": (
+                None if segment_data is None else segment_data["segment_length_m"]
+            ),
+            "segment_mass_flow_rate": (
+                None
+                if segment_data is None
+                else segment_data["segment_mass_flow_rate"]
+            ),
+            "segments": None if segment_data is None else segment_data["segments"],
+            "extra_parameters": dict(kwargs),
+        }
+
+    def calculate_recuperator_mass(self, component, **kwargs):
+        return self.calculate_mass(component, **kwargs)
+
+    @staticmethod
+    def _density_kg_m3(property_calculator, T, P):
+        if property_calculator is None:
+            return None
+        properties = property_calculator.calculate_properties(T=T, P=P)
+        if "density_kg_m3" in properties:
+            return properties["density_kg_m3"]
+        if "rho" in properties:
+            return properties["rho"]
+        if "density" in properties:
+            return properties["density"]
+        if "D" not in properties:
+            raise ValueError("property_calculator result must include density")
+        density = properties["D"]
+        mw = getattr(property_calculator, "mw", None)
+        if mw is None:
+            return density
+        return density * mw
+
+    @staticmethod
+    def _flow_velocity(mass_flow_rate, density_kg_m3, flow_area_m2):
+        if mass_flow_rate is None or density_kg_m3 is None:
+            return None
+        if density_kg_m3 <= 0:
+            raise ValueError("density_kg_m3 must be positive")
+        return mass_flow_rate / (density_kg_m3 * flow_area_m2)
+
+    def _side_flow_result(
+        self,
+        property_calculator,
+        inlet_state,
+        outlet_state,
+        mass_flow_rate,
+        flow_area_m2,
+        length_m,
+        hydraulic_diameter_m,
+    ):
+        outlet_state = dict(outlet_state)
+        mean_temperature = (inlet_state["T"] + outlet_state["T"]) / 2
+        mean_pressure = (inlet_state["P"] + outlet_state["P"]) / 2
+        density = self._density_kg_m3(
+            property_calculator,
+            T=mean_temperature,
+            P=mean_pressure,
+        )
+        velocity = self._flow_velocity(mass_flow_rate, density, flow_area_m2)
+        pressure_drop_kpa = self._pressure_drop_kpa(
+            property_calculator=property_calculator,
+            T=mean_temperature,
+            P=mean_pressure,
+            density_kg_m3=density,
+            velocity_m_s=velocity,
+            length_m=length_m,
+            hydraulic_diameter_m=hydraulic_diameter_m,
+        )
+
+        if pressure_drop_kpa is not None:
+            outlet_state["P"] = self.calculate_outlet_pressure(
+                inlet_state["P"],
+                pressure_drop_kpa,
+            )
+            mean_pressure = (inlet_state["P"] + outlet_state["P"]) / 2
+            density = self._density_kg_m3(
+                property_calculator,
+                T=mean_temperature,
+                P=mean_pressure,
+            )
+            velocity = self._flow_velocity(mass_flow_rate, density, flow_area_m2)
+
+        convection = self._convection_result(
+            property_calculator,
+            T=mean_temperature,
+            P=mean_pressure,
+            density_kg_m3=density,
+            velocity_m_s=velocity,
+            hydraulic_diameter_m=hydraulic_diameter_m,
+        )
+        return {
+            "outlet_state": outlet_state,
+            "mean_temperature": mean_temperature,
+            "mean_pressure": mean_pressure,
+            "flow_velocity_m_s": velocity,
+            "heat_transfer_coefficient": (
+                convection["heat_transfer_coefficient"]
+            ),
+        }
+
+    @staticmethod
+    def calculate_outlet_pressure(inlet_pressure_kpa, pressure_drop_kpa):
+        outlet_pressure_kpa = inlet_pressure_kpa - pressure_drop_kpa
+        if outlet_pressure_kpa <= 0:
+            raise ValueError("Calculated outlet pressure must be positive")
+        return outlet_pressure_kpa
+
+    def calculate_segment_outlet_pressure(
+        self,
+        inlet_pressure_kpa,
+        property_calculator,
+        T,
+        P,
+        density_kg_m3,
+        velocity_m_s,
+        length_m,
+        hydraulic_diameter_m,
+    ):
+        pressure_drop_kpa = self._pressure_drop_kpa(
+            property_calculator=property_calculator,
+            T=T,
+            P=P,
+            density_kg_m3=density_kg_m3,
+            velocity_m_s=velocity_m_s,
+            length_m=length_m,
+            hydraulic_diameter_m=hydraulic_diameter_m,
+        )
+        if pressure_drop_kpa is None:
+            return None
+        return self.calculate_outlet_pressure(inlet_pressure_kpa, pressure_drop_kpa)
+
+    def _pressure_drop_kpa(
+        self,
+        property_calculator,
+        T,
+        P,
+        density_kg_m3,
+        velocity_m_s,
+        length_m,
+        hydraulic_diameter_m,
+    ):
+        if (
+            property_calculator is None
+            or density_kg_m3 is None
+            or velocity_m_s is None
+        ):
+            return None
+        properties = property_calculator.calculate_properties(T=T, P=P)
+        viscosity = self._pick_positive_property(
+            properties,
+            (
+                "viscosity_pa_s",
+                "viscosity",
+                "dynamic_viscosity",
+                "mu",
+            ),
+        )
+        reynolds = (
+            density_kg_m3
+            * velocity_m_s
+            * hydraulic_diameter_m
+            / viscosity
+        )
+        friction_factor = self._pressure_drop_friction_factor(reynolds)
+        pressure_drop_pa = (
+            friction_factor
+            * length_m
+            / hydraulic_diameter_m
+            * density_kg_m3
+            * velocity_m_s**2
+            / 2
+        )
+        return pressure_drop_pa / 1000
+
+    @staticmethod
+    def _pressure_drop_friction_factor(reynolds):
+        if reynolds <= 0:
+            raise ValueError("Reynolds number must be positive")
+        if reynolds < 2300:
+            return 64 / reynolds
+        return 0.3164 / reynolds**0.25
+
+    def _convection_result(
+        self,
+        property_calculator,
+        T,
+        P,
+        density_kg_m3,
+        velocity_m_s,
+        hydraulic_diameter_m,
+    ):
+        if property_calculator is None or density_kg_m3 is None:
+            return {
+                "heat_transfer_coefficient": None,
+            }
+        if velocity_m_s is None:
+            return {
+                "heat_transfer_coefficient": None,
+            }
+
+        properties = property_calculator.calculate_properties(T=T, P=P)
+        viscosity = self._pick_positive_property(
+            properties,
+            (
+                "viscosity_pa_s",
+                "viscosity",
+                "dynamic_viscosity",
+                "mu",
+            ),
+        )
+        thermal_conductivity = self._pick_positive_property(
+            properties,
+            (
+                "thermal_conductivity_w_m_k",
+                "thermal_conductivity",
+                "conductivity",
+                "lambda",
+                "k",
+            ),
+        )
+        cp = self._specific_heat_j_kg_k(properties, property_calculator)
+
+        reynolds = (
+            density_kg_m3
+            * velocity_m_s
+            * hydraulic_diameter_m
+            / viscosity
+        )
+        prandtl = cp * viscosity / thermal_conductivity
+        nusselt, _ = self._nusselt_and_friction_factor(
+            reynolds,
+            prandtl,
+        )
+        return {
+            "heat_transfer_coefficient": (
+                nusselt * thermal_conductivity / hydraulic_diameter_m
+            ),
+        }
+
+    @staticmethod
+    def _nusselt_and_friction_factor(reynolds, prandtl):
+        if reynolds <= 0:
+            raise ValueError("Reynolds number must be positive")
+        if prandtl <= 0:
+            raise ValueError("Prandtl number must be positive")
+        if reynolds < 2300:
+            return 4.36, 64 / reynolds
+
+        friction_factor = 1 / (1.82 * math.log10(reynolds) - 1.64) ** 2
+        numerator = (friction_factor / 8) * (reynolds - 1000) * prandtl
+        denominator = (
+            1
+            + 12.7
+            * math.sqrt(friction_factor / 8)
+            * (prandtl ** (2 / 3) - 1)
+        )
+        return numerator / denominator, friction_factor
+
+    @staticmethod
+    def _pick_positive_property(properties, keys):
+        for key in keys:
+            if key in properties:
+                value = properties[key]
+                if value <= 0:
+                    raise ValueError(f"{key} must be positive")
+                return value
+        raise ValueError(
+            "property_calculator result must include one of: "
+            + ", ".join(keys)
+        )
+
+    @staticmethod
+    def _specific_heat_j_kg_k(properties, property_calculator):
+        for key in ("cp_j_kg_k", "specific_heat_j_kg_k", "cp_mass"):
+            if key in properties:
+                value = properties[key]
+                if value <= 0:
+                    raise ValueError(f"{key} must be positive")
+                return value
+        if "cp" not in properties:
+            raise ValueError("property_calculator result must include cp")
+
+        cp = properties["cp"]
+        if cp <= 0:
+            raise ValueError("cp must be positive")
+        mw = getattr(property_calculator, "mw", None)
+        if mw is None:
+            return cp
+        return cp / mw * 1000
+
+    @staticmethod
+    def _semicircular_channel_geometry(channel_diameter_m):
+        if channel_diameter_m <= 0:
+            raise ValueError("channel_diameter_m must be positive")
+        flow_area_m2 = math.pi * channel_diameter_m**2 / 8
+        wetted_perimeter_m = math.pi * channel_diameter_m / 2 + channel_diameter_m
+        hydraulic_diameter_m = 4 * flow_area_m2 / wetted_perimeter_m
+        return {
+            "channel_diameter_m": channel_diameter_m,
+            "flow_area_m2": flow_area_m2,
+            "wetted_perimeter_m": wetted_perimeter_m,
+            "hydraulic_diameter_m": hydraulic_diameter_m,
+        }
+
+    @staticmethod
+    def _state_from(variables, key, override):
+        state = override if override is not None else variables.get(key)
+        if state is None:
+            raise ValueError(f"Missing {key} for heat exchanger mass calculation")
+        return dict(state)
+
+    @staticmethod
+    def _interpolate_state(start, end, fraction):
+        return {
+            key: start[key] + (end[key] - start[key]) * fraction
+            for key in ("T", "P")
+        }
+
+
+_DEFAULT_RECUPERATOR_PCHE_MODEL = RecuperatorPCHEMassModel()
+
+
+def initialize_recuperator_segments(*args, **kwargs):
+    return _DEFAULT_RECUPERATOR_PCHE_MODEL.initialize_segments(*args, **kwargs)
+
+
+def recuperator_pche_mass(*args, **kwargs):
+    return _DEFAULT_RECUPERATOR_PCHE_MODEL.calculate_mass(*args, **kwargs)
+
+
+def calculate_water_outlet_temperature(
+    Qc_kw,
+    water_inlet_T,
+    water_mass_flow_rate,
+    water_pressure_kpa,
+    refprop_path,
+):
+    """Calculate substitute water outlet temperature with REFPROP."""
+    if water_inlet_T is None or water_mass_flow_rate is None:
+        raise ValueError(
+            "Provide coolant_outlet_T directly, or provide water_inlet_T "
+            "and water_mass_flow_rate for REFPROP water calculation."
+        )
+    if water_mass_flow_rate <= 0:
+        raise ValueError("water_mass_flow_rate must be positive")
+
+    import ctREFPROP.ctREFPROP as ct
+
+    water = ct.REFPROPFunctionLibrary(refprop_path)
+    water.SETUPdll(1, "WATER.FLD", "HMX.BNC", "DEF")
+    water.SETUPdll(2, "SI", "SI", "DEF")
+    z = [1.0]
+    mw = water.WMOLdll(z)
+
+    inlet = water.TPFLSHdll(water_inlet_T, water_pressure_kpa, z)
+    if inlet.ierr > 0:
+        raise ValueError(f"REFPROP water inlet calculation error: {inlet.ierr}")
+
+    h_in_mass = inlet.h / mw
+    h_out_mass = h_in_mass + Qc_kw / water_mass_flow_rate
+    outlet = water.PHFLSHdll(water_pressure_kpa, h_out_mass * mw, z)
+    if outlet.ierr > 0:
+        raise ValueError(f"REFPROP water outlet calculation error: {outlet.ierr}")
+
+    return outlet.T
+
+
+DEFAULT_MASS_MODELS = {
+    "turbine": turbine_tac_mass,
+    "heater": heater_reactor_mass,
+    "condenser": condenser_radiator_mass,
+    "recuperator": recuperator_pche_mass,
+}
+
+
+def has_default_mass_model(component_or_type):
+    """Return whether a component type has a built-in mass model."""
+    component_type = (
+        component_or_type
+        if isinstance(component_or_type, str)
+        else getattr(component_or_type, "component_type", None)
+    )
+    return component_type in DEFAULT_MASS_MODELS
+
+
+def calculate_component_mass(component, mass_model=None, **kwargs):
+    """Calculate mass for one component without mutating it."""
+    if component.variables is None:
+        raise ValueError("Run component calculator before mass calculation")
+
+    if mass_model is None:
+        try:
+            model_function = DEFAULT_MASS_MODELS[component.component_type]
+        except KeyError as exc:
+            raise NotImplementedError(
+                f"No default mass model for component type "
+                f"{component.component_type!r}."
+            ) from exc
+        return model_function(component.variables, **kwargs)
+
+    return _run_external_mass_model(component, mass_model, **kwargs)
+
+
+def apply_component_mass(component, mass_model=None, **kwargs):
+    """Calculate mass with a mass model and store it on the component."""
+    result = calculate_component_mass(component, mass_model=mass_model, **kwargs)
+    return component.set_mass_result(result)
+
+
+def _run_external_mass_model(component, mass_model, **kwargs):
+    method_name = f"calculate_{component.component_type}_mass"
+    if hasattr(mass_model, method_name):
+        return getattr(mass_model, method_name)(component, **kwargs)
+    if hasattr(mass_model, "calculate_component_mass"):
+        return mass_model.calculate_component_mass(component, **kwargs)
+    if hasattr(mass_model, "calculate_mass"):
+        return mass_model.calculate_mass(component, **kwargs)
+    if callable(mass_model):
+        return mass_model(component, **kwargs)
+    raise TypeError("mass_model must be callable or expose a supported method")

brayton_cycle/properties.py

@@ -60,6 +60,15 @@ class CO2PropertyCalculator():
         
         # 补充提取的物性，这里由于后续还要使用，不进行参数变换
         properties['D'] = result.D
+        properties['density_kg_m3'] = result.D * self.mw
         properties['cp'] = result.Cp
-        properties['cv'] = result.Cv,     
+        properties['cp_j_kg_k'] = result.Cp / self.mw * 1000
+        properties['cv'] = result.Cv
+        properties['cv_j_kg_k'] = result.Cv / self.mw * 1000
+
+        transport = self.rp.TRNPRPdll(result.T, result.D, self.z)
+        if transport.ierr > 0:
+            raise ValueError(f"REFPROP transport calculation error:{transport.ierr}")
+        properties['viscosity_pa_s'] = transport.eta * 1e-6
+        properties['thermal_conductivity_w_m_k'] = transport.tcx
         return properties