# -*- 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")