From facb914417582494c1eddc1d2c46b7735359aba2 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Wed, 24 Jun 2026 21:55:52 +0800 Subject: [PATCH] =?UTF-8?q?=E6=8D=A2=E7=83=AD=E5=99=A8=E8=B4=A8=E9=87=8F?= =?UTF-8?q?=E6=A8=A1=E5=9E=8B=E5=BB=BA=E7=AB=8B-=E5=85=B3=E8=81=94?= =?UTF-8?q?=E5=BC=8F=E5=87=BD=E6=95=B0=E5=AE=9E=E7=8E=B0=E3=80=81PCHE?= =?UTF-8?q?=E7=BB=93=E6=9E=84=E5=8F=82=E6=95=B0=E5=88=9D=E5=A7=8B=E5=8C=96?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- brayton_cycle/__init__.py | 34 +++++ brayton_cycle/mass_models.py | 288 +++++++++++++++++++++++++++++++++++ brayton_cycle/properties.py | 32 +++- 3 files changed, 351 insertions(+), 3 deletions(-) diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py index 8600e2f..f1eea90 100644 --- a/brayton_cycle/__init__.py +++ b/brayton_cycle/__init__.py @@ -11,10 +11,27 @@ from .components import ( ) 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, @@ -39,19 +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", ] diff --git a/brayton_cycle/mass_models.py b/brayton_cycle/mass_models.py index 6ea0784..54fdffc 100644 --- a/brayton_cycle/mass_models.py +++ b/brayton_cycle/mass_models.py @@ -1,9 +1,297 @@ # -*- 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. diff --git a/brayton_cycle/properties.py b/brayton_cycle/properties.py index 231487b..8f4de59 100644 --- a/brayton_cycle/properties.py +++ b/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 @@ -60,6 +60,32 @@ class CO2PropertyCalculator(): # 补充提取的物性，这里由于后续还要使用，不进行参数变换 properties['D'] = result.D - properties['cp'] = result.Cp - properties['cv'] = result.Cv, + 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