diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py index 6cdbf65..8600e2f 100644 --- a/brayton_cycle/__init__.py +++ b/brayton_cycle/__init__.py @@ -10,6 +10,12 @@ from .components import ( Turbine, ) from .cycles import BraytonCycle +from .mass_models import ( + calculate_condenser_mass, + calculate_heater_mass, + calculate_turbine_mass, + calculate_water_outlet_temperature, +) from .optimization import ( optimize_rc_fixed_param, optimize_rc_param, @@ -35,6 +41,10 @@ __all__ = [ "Heater", "Recuperator", "Turbine", + "calculate_condenser_mass", + "calculate_heater_mass", + "calculate_turbine_mass", + "calculate_water_outlet_temperature", "evaluate_rc_efficiency", "local_rc_component_performance_sensitivity", "local_rc_design_sensitivity", diff --git a/brayton_cycle/components.py b/brayton_cycle/components.py index 55ffb9c..bfe945d 100644 --- a/brayton_cycle/components.py +++ b/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." - ) - 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', - } + 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, + ) class Condenser(ComponentMassMixin): """冷凝器类""" @@ -477,112 +419,20 @@ 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( - 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 + 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, + ) class Concentrator(ComponentMassMixin): """汇流组件""" diff --git a/brayton_cycle/mass_models.py b/brayton_cycle/mass_models.py new file mode 100644 index 0000000..6ea0784 --- /dev/null +++ b/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