From 1a4b050d5336fdf06f712959ec90b6592b7a5c24 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Wed, 24 Jun 2026 11:23:58 +0800 Subject: [PATCH] =?UTF-8?q?=E9=87=8D=E6=9E=84=E8=B4=A8=E9=87=8F=E6=A8=A1?= =?UTF-8?q?=E5=9E=8B=E5=B9=B6=E8=A1=A5=E5=85=85=E5=9B=9E=E7=83=AD=E5=99=A8?= =?UTF-8?q?=E5=88=86=E6=AE=B5=E5=BB=BA=E6=A8=A1?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- brayton_cycle/__init__.py | 20 + brayton_cycle/components.py | 258 +---------- brayton_cycle/cycles.py | 36 +- brayton_cycle/mass_models.py | 861 +++++++++++++++++++++++++++++++++++ brayton_cycle/properties.py | 13 +- 5 files changed, 940 insertions(+), 248 deletions(-) create mode 100644 brayton_cycle/mass_models.py diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py index 6cdbf65..f5bbf78 100644 --- a/brayton_cycle/__init__.py +++ b/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", ] diff --git a/brayton_cycle/components.py b/brayton_cycle/components.py index 55ffb9c..2a70006 100644 --- a/brayton_cycle/components.py +++ b/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__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" diff --git a/brayton_cycle/cycles.py b/brayton_cycle/cycles.py index 3451c30..da558c2 100644 --- a/brayton_cycle/cycles.py +++ b/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.""" diff --git a/brayton_cycle/mass_models.py b/brayton_cycle/mass_models.py new file mode 100644 index 0000000..0e18ca8 --- /dev/null +++ b/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") diff --git a/brayton_cycle/properties.py b/brayton_cycle/properties.py index 231487b..d6a8999 100644 --- a/brayton_cycle/properties.py +++ b/brayton_cycle/properties.py @@ -60,6 +60,15 @@ class CO2PropertyCalculator(): # 补充提取的物性，这里由于后续还要使用，不进行参数变换 properties['D'] = result.D - properties['cp'] = result.Cp - properties['cv'] = result.Cv, + properties['density_kg_m3'] = result.D * self.mw + properties['cp'] = result.Cp + 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