换热器质量模型建立-关联式函数实现、PCHE结构参数初始化

brayton_cycle/__init__.py

@@ -10,6 +10,29 @@ from .components import (
     Turbine,
 )
 from .cycles import BraytonCycle
+from .mass_models import (
+    PCHEGeometry,
+    PCHEMaterial,
+    PCHEDesignOptions,
+    PCHE_STATE_UNIT_CONVENTION,
+    calculate_condenser_mass,
+    calculate_heater_mass,
+    calculate_turbine_mass,
+    calculate_water_outlet_temperature,
+    flow_velocity,
+    friction_factor,
+    gnielinski_f_star,
+    heat_transfer_coefficient,
+    nusselt_number,
+    overall_heat_transfer_coefficient,
+    pche_channel_flow_area,
+    pche_channel_wetted_perimeter,
+    pche_equivalent_diameter,
+    prandtl_number,
+    pressure_drop_pa,
+    reynolds_number,
+    validate_pche_state,
+)
 from .optimization import (
     optimize_rc_fixed_param,
     optimize_rc_param,
@@ -33,15 +56,36 @@ __all__ = [
     "Concentrator",
     "Condenser",
     "Heater",
+    "PCHEGeometry",
+    "PCHEMaterial",
+    "PCHEDesignOptions",
+    "PCHE_STATE_UNIT_CONVENTION",
     "Recuperator",
     "Turbine",
+    "calculate_condenser_mass",
+    "calculate_heater_mass",
+    "calculate_turbine_mass",
+    "calculate_water_outlet_temperature",
+    "flow_velocity",
+    "friction_factor",
+    "gnielinski_f_star",
+    "heat_transfer_coefficient",
     "evaluate_rc_efficiency",
     "local_rc_component_performance_sensitivity",
     "local_rc_design_sensitivity",
+    "nusselt_number",
     "optimize_rc_fixed_param",
     "optimize_rc_param",
+    "overall_heat_transfer_coefficient",
+    "pche_channel_flow_area",
+    "pche_channel_wetted_perimeter",
+    "pche_equivalent_diameter",
     "plot_optimization_landscape",
     "plot_sweep_optimization_results",
+    "prandtl_number",
+    "pressure_drop_pa",
+    "reynolds_number",
     "scan_rc_efficiency",
     "sweep_and_optimize_rc",
+    "validate_pche_state",
 ]

brayton_cycle/components.py

@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
 """Component models used by Brayton cycle simulations."""
 
-import math
-
-import ctREFPROP.ctREFPROP as ct
+from .mass_models import (
+    calculate_condenser_mass,
+    calculate_heater_mass,
+    calculate_turbine_mass,
+)
 
 
 class ComponentMassMixin:
@@ -166,38 +168,7 @@ class Turbine(ComponentMassMixin):
             }
 
     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))',
-        }
+        return calculate_turbine_mass(self, Pe=Pe, A=A, **kwargs)
 
 class Recuperator(ComponentMassMixin):
     """换热器类"""
@@ -397,43 +368,14 @@ class Heater(ComponentMassMixin):
         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."
+        return calculate_heater_mass(
+            self,
+            P_heat_mwt=P_heat_mwt,
+            mass_flow_rate=mass_flow_rate,
+            shielding_mass_ton=shielding_mass_ton,
+            include_shielding=include_shielding,
+            **kwargs,
         )
-            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):
     """冷凝器类"""
@@ -477,113 +419,21 @@ class Condenser(ComponentMassMixin):
         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(
+        return calculate_condenser_mass(
+            self,
             Qc_kw=Qc_kw,
+            co2_mass_flow_rate=co2_mass_flow_rate,
+            coolant_outlet_T=coolant_outlet_T,
             water_inlet_T=water_inlet_T,
             water_mass_flow_rate=water_mass_flow_rate,
             water_pressure_kpa=water_pressure_kpa,
             refprop_path=refprop_path,
+            emissivity=emissivity,
+            surface_temperature=surface_temperature,
+            area_density=area_density,
+            **kwargs,
         )
 
-        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/mass_models.py

@@ -0,0 +1,484 @@
+# -*- coding: utf-8 -*-
+"""Mass-estimation models for Brayton cycle components."""
+
+from dataclasses import dataclass
+import math
+
+
+PCHE_STATE_UNIT_CONVENTION = {
+    "P": "kPa",
+    "T": "K",
+    "h": "kJ/kg",
+    "mass_flow_rate": "kg/s",
+    "density": "kg/m3",
+    "viscosity": "Pa*s",
+    "thermal_conductivity": "W/(m*K)",
+    "specific_heat": "J/(kg*K)",
+    "length": "m",
+    "area": "m2",
+    "pressure_drop": "Pa",
+    "mass": "kg",
+}
+
+
+@dataclass(frozen=True)
+class PCHEGeometry:
+    """Default PCHE channel geometry from Yuan et al."""
+
+    channel_diameter_m: float = 0.002
+    channel_pitch_m: float = 0.0024
+    plate_thickness_m: float = 0.0015
+
+    def __post_init__(self):
+        _require_positive("channel_diameter_m", self.channel_diameter_m)
+        _require_positive("channel_pitch_m", self.channel_pitch_m)
+        _require_positive("plate_thickness_m", self.plate_thickness_m)
+
+    @property
+    def flow_area_m2(self):
+        return pche_channel_flow_area(self)
+
+    @property
+    def wetted_perimeter_m(self):
+        return pche_channel_wetted_perimeter(self)
+
+    @property
+    def equivalent_diameter_m(self):
+        return pche_equivalent_diameter(self)
+
+
+@dataclass(frozen=True)
+class PCHEMaterial:
+    """Default Inconel 617 material data from Yuan et al."""
+
+    density_kg_m3: float = 8360.0
+    thermal_conductivity_w_m_k: float = 21.0
+
+    def __post_init__(self):
+        _require_positive("density_kg_m3", self.density_kg_m3)
+        _require_positive(
+            "thermal_conductivity_w_m_k",
+            self.thermal_conductivity_w_m_k,
+        )
+
+
+@dataclass(frozen=True)
+class PCHEDesignOptions:
+    """Numerical and design limits for the Yuan-style PCHE calculation."""
+
+    num_segments: int = 50
+    allowable_pressure_drop_ratio: float = 0.01
+    outlet_temperature_tolerance_k: float = 1.0e-3
+    max_iterations: int = 100
+    reynolds_transition: float = 2300.0
+
+    def __post_init__(self):
+        if self.num_segments < 1:
+            raise ValueError("num_segments must be at least 1")
+        if self.max_iterations < 1:
+            raise ValueError("max_iterations must be at least 1")
+        _require_positive(
+            "allowable_pressure_drop_ratio",
+            self.allowable_pressure_drop_ratio,
+        )
+        _require_positive(
+            "outlet_temperature_tolerance_k",
+            self.outlet_temperature_tolerance_k,
+        )
+        _require_positive("reynolds_transition", self.reynolds_transition)
+
+
+def _require_positive(name, value):
+    _require_finite(name, value)
+    if value <= 0:
+        raise ValueError(f"{name} must be positive")
+
+
+def _require_finite(name, value):
+    if not math.isfinite(value):
+        raise ValueError(f"{name} must be finite")
+
+
+def validate_pche_state(state, state_name="state", require_transport=False):
+    """Validate the state fields expected by the PCHE mass model."""
+    required = ("P", "T", "h")
+    missing = [key for key in required if key not in state]
+    if missing:
+        raise ValueError(f"{state_name} missing required keys: {', '.join(missing)}")
+
+    _require_positive(f"{state_name}.P", state["P"])
+    _require_positive(f"{state_name}.T", state["T"])
+    _require_finite(f"{state_name}.h", state["h"])
+
+    if require_transport:
+        transport_keys = (
+            "rho",
+            "cp_mass",
+            "mu",
+            "thermal_conductivity",
+        )
+        missing = [key for key in transport_keys if key not in state]
+        if missing:
+            raise ValueError(
+                f"{state_name} missing transport keys: {', '.join(missing)}"
+            )
+        for key in transport_keys:
+            _require_positive(f"{state_name}.{key}", state[key])
+
+
+def pche_channel_flow_area(geometry=None):
+    """Return the semicircular PCHE channel flow area in m2."""
+    geometry = geometry or PCHEGeometry()
+    diameter = geometry.channel_diameter_m
+    return math.pi * diameter**2 / 8.0
+
+
+def pche_channel_wetted_perimeter(geometry=None):
+    """Return the Yuan PCHE channel wetted perimeter in m."""
+    geometry = geometry or PCHEGeometry()
+    diameter = geometry.channel_diameter_m
+    return math.pi * diameter / 2.0 + diameter
+
+
+def pche_equivalent_diameter(geometry=None):
+    """Return the Yuan PCHE equivalent diameter in m."""
+    area = pche_channel_flow_area(geometry)
+    perimeter = pche_channel_wetted_perimeter(geometry)
+    return 4.0 * area / perimeter
+
+
+def flow_velocity(
+    mass_flow_rate_kg_s,
+    density_kg_m3,
+    flow_area_m2,
+    parallel_channels=1,
+):
+    """Return average channel velocity in m/s."""
+    _require_positive("mass_flow_rate_kg_s", mass_flow_rate_kg_s)
+    _require_positive("density_kg_m3", density_kg_m3)
+    _require_positive("flow_area_m2", flow_area_m2)
+    if parallel_channels < 1:
+        raise ValueError("parallel_channels must be at least 1")
+
+    channel_mass_flow = mass_flow_rate_kg_s / parallel_channels
+    return channel_mass_flow / (density_kg_m3 * flow_area_m2)
+
+
+def reynolds_number(
+    density_kg_m3,
+    velocity_m_s,
+    hydraulic_diameter_m,
+    viscosity_pa_s,
+):
+    """Return Reynolds number."""
+    _require_positive("density_kg_m3", density_kg_m3)
+    _require_positive("velocity_m_s", velocity_m_s)
+    _require_positive("hydraulic_diameter_m", hydraulic_diameter_m)
+    _require_positive("viscosity_pa_s", viscosity_pa_s)
+    return density_kg_m3 * velocity_m_s * hydraulic_diameter_m / viscosity_pa_s
+
+
+def prandtl_number(
+    specific_heat_j_kg_k,
+    viscosity_pa_s,
+    thermal_conductivity_w_m_k,
+):
+    """Return Prandtl number."""
+    _require_positive("specific_heat_j_kg_k", specific_heat_j_kg_k)
+    _require_positive("viscosity_pa_s", viscosity_pa_s)
+    _require_positive(
+        "thermal_conductivity_w_m_k",
+        thermal_conductivity_w_m_k,
+    )
+    return specific_heat_j_kg_k * viscosity_pa_s / thermal_conductivity_w_m_k
+
+
+def gnielinski_f_star(reynolds):
+    """Return the turbulent f* term used in Yuan et al.'s Nu correlation."""
+    _require_positive("reynolds", reynolds)
+    denominator = 1.82 * math.log10(reynolds) - 1.64
+    if denominator == 0:
+        raise ValueError("reynolds gives a zero denominator in f* correlation")
+    return 1.0 / denominator**2
+
+
+def nusselt_number(reynolds, prandtl, transition_re=2300.0):
+    """Return Nusselt number using Yuan et al.'s laminar/turbulent formulas."""
+    _require_positive("reynolds", reynolds)
+    _require_positive("prandtl", prandtl)
+    _require_positive("transition_re", transition_re)
+
+    if reynolds < transition_re:
+        return 4.36
+
+    f_star = gnielinski_f_star(reynolds)
+    numerator = (f_star / 8.0) * (reynolds - 1000.0) * prandtl
+    denominator = 1.0 + 12.7 * math.sqrt(f_star / 8.0) * (
+        prandtl ** (2.0 / 3.0) - 1.0
+    )
+    if denominator == 0:
+        raise ValueError("Nusselt correlation denominator is zero")
+    return numerator / denominator
+
+
+def heat_transfer_coefficient(
+    nusselt,
+    thermal_conductivity_w_m_k,
+    hydraulic_diameter_m,
+):
+    """Return convective heat-transfer coefficient in W/(m2*K)."""
+    _require_positive("nusselt", nusselt)
+    _require_positive(
+        "thermal_conductivity_w_m_k",
+        thermal_conductivity_w_m_k,
+    )
+    _require_positive("hydraulic_diameter_m", hydraulic_diameter_m)
+    return nusselt * thermal_conductivity_w_m_k / hydraulic_diameter_m
+
+
+def overall_heat_transfer_coefficient(
+    hot_h_w_m2_k,
+    cold_h_w_m2_k,
+    wall_thickness_m,
+    wall_thermal_conductivity_w_m_k,
+):
+    """Return total heat-transfer coefficient K in W/(m2*K)."""
+    _require_positive("hot_h_w_m2_k", hot_h_w_m2_k)
+    _require_positive("cold_h_w_m2_k", cold_h_w_m2_k)
+    _require_positive("wall_thickness_m", wall_thickness_m)
+    _require_positive(
+        "wall_thermal_conductivity_w_m_k",
+        wall_thermal_conductivity_w_m_k,
+    )
+
+    resistance = (
+        1.0 / hot_h_w_m2_k
+        + 1.0 / cold_h_w_m2_k
+        + wall_thickness_m / wall_thermal_conductivity_w_m_k
+    )
+    return 1.0 / resistance
+
+
+def friction_factor(reynolds, transition_re=2300.0):
+    """Return Darcy friction factor from Yuan et al.'s piecewise relation."""
+    _require_positive("reynolds", reynolds)
+    _require_positive("transition_re", transition_re)
+    if reynolds < transition_re:
+        return 64.0 / reynolds
+    return 0.3164 / reynolds**0.25
+
+
+def pressure_drop_pa(
+    friction_factor_value,
+    length_m,
+    hydraulic_diameter_m,
+    density_kg_m3,
+    velocity_m_s,
+):
+    """Return channel pressure drop in Pa."""
+    _require_positive("friction_factor_value", friction_factor_value)
+    _require_positive("length_m", length_m)
+    _require_positive("hydraulic_diameter_m", hydraulic_diameter_m)
+    _require_positive("density_kg_m3", density_kg_m3)
+    _require_positive("velocity_m_s", velocity_m_s)
+
+    return (
+        friction_factor_value
+        * length_m
+        / hydraulic_diameter_m
+        * density_kg_m3
+        * velocity_m_s**2
+        / 2.0
+    )
+
+
+def calculate_turbine_mass(component, 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(component.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))",
+    }
+
+
+def calculate_heater_mass(
+    component,
+    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(component.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",
+    }
+
+
+def calculate_condenser_mass(
+    component,
+    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.
+    """
+    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(component.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",
+        "formula": "M_rad=kappa*Qc/(phi*sigma*(T^4-T0^4))",
+    }
+
+
+def calculate_water_outlet_temperature(
+    Qc_kw,
+    water_inlet_T,
+    water_mass_flow_rate,
+    water_pressure_kpa,
+    refprop_path,
+):
+    """Calculate water outlet temperature for the condenser mass model."""
+    import ctREFPROP.ctREFPROP as ct
+
+    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

brayton_cycle/properties.py

@@ -15,7 +15,7 @@ class CO2PropertyCalculator():
         self.rp.SETUPdll(2, 'SI', 'SI', 'DEF')
         self.z = [1.0]
         self.mw = self.rp.WMOLdll(self.z)
-    def calculate_properties(self, T=None, P=None, h=None, s=None):
+    def calculate_properties(self, T=None, P=None, h=None, s=None, include_transport=False):
         """计算二氧化碳物性"""
         if T is not None and P is not None:
             # 已知Tp
@@ -62,4 +62,30 @@ class CO2PropertyCalculator():
         properties['D'] = result.D
         properties['cp'] = result.Cp
         properties['cv'] = result.Cv,
+        properties['rho'] = result.D * self.mw
+        properties['cp_mass'] = result.Cp / self.mw * 1000.0
+        properties['cv_mass'] = result.Cv / self.mw * 1000.0
+        if include_transport:
+            self._add_transport_properties(properties, result)
         return properties
+
+    def _add_transport_properties(self, properties, result):
+        """Add CO2 transport properties needed by PCHE correlations."""
+        transport = self.rp.TRNPRPdll(properties['T'], result.D, self.z)
+        if getattr(transport, 'ierr', 0) > 0:
+            raise ValueError(f"REFPROP transport calculation error:{transport.ierr}")
+
+        viscosity_micro_pa_s = getattr(transport, 'eta', None)
+        thermal_conductivity = getattr(transport, 'tcx', None)
+        if viscosity_micro_pa_s is None:
+            viscosity_micro_pa_s = getattr(transport, 'visc', None)
+        if thermal_conductivity is None:
+            thermal_conductivity = getattr(transport, 'tcond', None)
+        if viscosity_micro_pa_s is None or thermal_conductivity is None:
+            raise ValueError("REFPROP transport result missing viscosity or conductivity")
+
+        viscosity_pa_s = viscosity_micro_pa_s * 1.0e-6
+        properties['mu'] = viscosity_pa_s
+        properties['viscosity'] = viscosity_pa_s
+        properties['lambda'] = thermal_conductivity
+        properties['thermal_conductivity'] = thermal_conductivity