架构调整，质量计算模型分块包装

brayton_cycle/mass_models.py

@@ -0,0 +1,196 @@
+# -*- coding: utf-8 -*-
+"""Mass-estimation models for Brayton cycle components."""
+
+import math
+
+
+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