From 1d09599e8ccc6d4deaa457d03efab3fade7adb56 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Sat, 20 Jun 2026 14:12:51 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9E=B6=E6=9E=84=E8=B0=83=E6=95=B4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 6 +- brayton_cycle/__init__.py | 23 + brayton_cycle/components.py | 331 ++++++++++++ brayton_cycle/cycles.py | 475 ++++++++++++++++ brayton_cycle/properties.py | 65 +++ brayton_cycle_test.py | 1019 ++--------------------------------- 6 files changed, 948 insertions(+), 971 deletions(-) create mode 100644 brayton_cycle/__init__.py create mode 100644 brayton_cycle/components.py create mode 100644 brayton_cycle/cycles.py create mode 100644 brayton_cycle/properties.py diff --git a/.gitignore b/.gitignore index 0cafc1c..a1506c7 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,5 @@ -.venv/ \ No newline at end of file +.venv/ +.pip-cache/ +__pycache__/ +*.pyc +*.pyo diff --git a/brayton_cycle/__init__.py b/brayton_cycle/__init__.py new file mode 100644 index 0000000..2e26aab --- /dev/null +++ b/brayton_cycle/__init__.py @@ -0,0 +1,23 @@ +"""Tools for CO2 Brayton cycle simulation and optimization.""" + +from .components import ( + Compressor, + Concentrator, + Condenser, + Heater, + Recuperator, + Turbine, +) +from .cycles import BraytonCycle +from .properties import CO2PropertyCalculator + +__all__ = [ + "BraytonCycle", + "CO2PropertyCalculator", + "Compressor", + "Concentrator", + "Condenser", + "Heater", + "Recuperator", + "Turbine", +] diff --git a/brayton_cycle/components.py b/brayton_cycle/components.py new file mode 100644 index 0000000..a89ddbb --- /dev/null +++ b/brayton_cycle/components.py @@ -0,0 +1,331 @@ +# -*- coding: utf-8 -*- +"""Component models used by Brayton cycle simulations.""" + +class Compressor(): + """压缩机类""" + def __init__(self, name, eff): + """ + 初始化参数 + + name: 名称 + eff: 等熵效率 + Wc: 压缩功 + inlet_state: 入口参数 + outlet_state: 出口参数 + outlet_state_is: 等熵状态下出口参数 + """ + + self.name = name + self.eff = eff + self.variables = None + + def calculator(self, p_in, T_in, p_out, property_calculator): + # 先计算熵值 + inlet_state = property_calculator.calculate_properties(T=T_in, P=p_in) + mw = property_calculator.mw + s = inlet_state['s'] + h_in = inlet_state['h'] + outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s) + h_out_is = outlet_state_is['h'] + h_out = h_in + (h_out_is - h_in) / self.eff + outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out) + Wc = h_out - h_in + T_out = outlet_state['T'] + # 计算结果 + self.variables = { + 'name': self.name, + 'inlet_state':{ + 'P': p_in, # 压强(kPa) + 'T': T_in, + 'h': h_in/mw, # 比焓(J/mol)->(J/kg) + 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) + }, + 'outlet_state':{ + 'P': p_out, # 压强(kPa) + 'T': T_out, + 'h': h_out/mw, # 比焓(kJ/mol)->(kJ/kg) + 's': s/mw, # 比熵(kJ/mol.K)->(kJ/kg.K) + }, + 'eff': self.eff, + 'Wc': Wc/mw, # 压缩功(kJ/mol)->(kJ/kg) + 'pi': p_out / p_in + } + +class Turbine(): + """透平类""" + def __init__(self, name, eff): + """ + 初始化参数 + name: 名称 + eff: 透平效率 + """ + self.name = name + self.eff = eff + self.variables = None + + def calculator(self, p_in, T_in, p_out, property_calculator): + """涡轮参数计算""" + inlet_state = property_calculator.calculate_properties(P=p_in, T=T_in) + mw = property_calculator.mw + s = inlet_state['s'] + h_in = inlet_state['h'] + outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s) + h_out_is = outlet_state_is['h'] + h_out = h_in - (h_in - h_out_is) * self.eff + outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out) + T_out = outlet_state['T'] + Wt = h_in - h_out + self.variables = { + 'name': self.name, + 'inlet_state':{ + 'P': p_in, # 压强(kPa) + 'T': T_in, + 'h': h_in/mw, # 比焓(J/mol)->(J/kg) + 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) + }, + 'outlet_state':{ + 'P': p_out, # 压强(kPa) + 'T': T_out, + 'h': h_out/mw, # 比焓(J/mol)->(J/kg) + 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) + }, + 'eff': self.eff, + 'Wt': Wt/mw, # 透平做功(J/mol)->(J/kg) + 'pi': p_in / p_out + } + +class Recuperator(): + """换热器类""" + def __init__(self, name, eff, x=0): + """ + 初始化参数 + name: 名称 + eff: 换热效率(基于焓的计算方法) + """ + self.name = name + self.eff = eff + self.Q_ex = None + self.variables = None + self.x = x + + def calculator(self, cold_inlet_state, hot_inlet_state, bypass_info, ploss=0.0, property_calculator=None): + """ + 计算换热器两侧参数，默认逆流 + bypass_info: 是否存在分流，0为不存在，1为存在 + + 下标含义 + ---------- + 1: 换热器冷端入口 + 2: 换热器冷端出口 + 3: 换热器热端入口 + 4: 换热器热端出口 + ass: 迭代中间变量，假设值 + """ + + # 计算两入口参数, 这里单位是kg + if property_calculator is None: + if hasattr(ploss, "calculate_properties"): + property_calculator = ploss + ploss = 0.0 + else: + raise ValueError("property_calculator is required") + + p1 = cold_inlet_state['P'] + T1 = cold_inlet_state['T'] + h1 = cold_inlet_state['h'] + p2 = p1 * (1 - ploss) + p3 = hot_inlet_state['P'] + T3 = hot_inlet_state['T'] + h3 = hot_inlet_state['h'] + p4 = p3 * (1 - ploss) + mw = property_calculator.mw + + # 设定质量流量 + m_cold = (1-self.x) if bypass_info == 1 else 1.0 + m_hot = 1.0 + +# ============================================================================= +# # 采用焓差效能的方式来计算换热器进出口参数 +# # 假设最大温差发生在冷端, 计算冷端出口温度和比焓 +# cold_outlet_state = property_calculator.calculate_properties(T=T3, P=p2) +# h2 = cold_outlet_state['h'] / mw +# Q_ass_cold = m_cold * abs(h2 - h1) +# +# # 假设最大温差发生在热端, 计算热端出口温度和比焓 +# hot_outlet_state = property_calculator.calculate_properties(T=T1, P=p4) +# h4 = hot_outlet_state['h'] / mw +# Q_ass_hot = m_hot * abs(h3 - h4) +# +# # 比较两个可能的Q，取最小值与焓差效能的乘积作为实际换热量 +# self.Q_ex = min(Q_ass_hot, Q_ass_cold) * self.eff +# # 由实际换热量计算出口焓和出口状态 +# h2 = h1 + self.Q_ex / m_cold +# h4 = h3 - self.Q_ex / m_hot +# cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw) +# hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw) +# T2 = cold_outlet_state['T'] +# T4 = hot_outlet_state['T'] +# ============================================================================= + + # 采用温差效能的方式来计算换热器进出口参数 + T_ass_max = abs(T1 - T3) + # 假设最大温差发生在冷端, 计算冷端出口温度和比焓 + T2 = T1 + T_ass_max * self.eff + cold_outlet_state = property_calculator.calculate_properties(T=T2, P=p2) + h2 = cold_outlet_state['h'] / mw + Q_ass_cold = m_cold * abs(h2 - h1) + + # 假设最大温差发生在热端, 计算热端出口温度和比焓 + T4 = T3 - T_ass_max * self.eff + hot_outlet_state = property_calculator.calculate_properties(T=T4, P=p4) + h4 = hot_outlet_state['h'] / mw + Q_ass_hot = m_hot * abs(h3 - h4) + + # 比较两个可能的Q，取最小值与焓差效能的乘积作为实际换热量 + self.Q_ex = min(Q_ass_hot, Q_ass_cold) + # 由实际换热量计算出口焓和出口状态 + h2 = h1 + self.Q_ex / m_cold + h4 = h3 - self.Q_ex / m_hot + cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw) + hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw) + T2 = cold_outlet_state['T'] + T4 = hot_outlet_state['T'] + + # 拼装变量 + self.variables = { + 'name': self.name, + 'cold_inlet_state':{ + 'P': p1, + 'T': T1, + 'h': h1, + 's': cold_inlet_state['s'] + }, + 'cold_outlet_state':{ + 'P': p2, + 'T': cold_outlet_state['T'], + 'h': h2, + 's': cold_outlet_state['s']/mw + }, + 'hot_inlet_state':{ + 'P': p3, + 'T': T3, + 'h': h3, + 's': hot_inlet_state['s'] + }, + 'hot_outlet_state':{ + 'P': p4, + 'T': T4, + 'h': h4, + 's': hot_outlet_state['s']/mw + }, + 'eff': self.eff, + 'Q_exchange': self.Q_ex + } + def check_pinch_point(self, property_calculator, num_segments=20): + """ + 换热器内部夹点校验 + 将换热量均分为 num_segments 段，检查内部每个微元的冷热流体温度 + """ + h_cold_in = self.variables['cold_inlet_state']['h'] + h_hot_in = self.variables['hot_inlet_state']['h'] + p_cold = self.variables['cold_inlet_state']['P'] + p_hot = self.variables['hot_inlet_state']['P'] + + m_cold = (1 - self.x) if self.name == "Low Temperature recuprerator" else 1.0 + m_hot = 1.0 + + dQ = self.Q_ex / num_segments + + # 沿冷流体流动方向步进检查 + for i in range(num_segments + 1): + q_current = i * dQ + + # 当前微元截面的焓值 + h_cold_local = h_cold_in + q_current / m_cold + h_hot_local = (h_hot_in - self.Q_ex / m_hot) + q_current / m_hot + + # 查温度 + T_cold_local = property_calculator.calculate_properties(P=p_cold, h=h_cold_local * property_calculator.mw)['T'] + T_hot_local = property_calculator.calculate_properties(P=p_hot, h=h_hot_local * property_calculator.mw)['T'] + + # 如果热流体温度低于等于冷流体温度 (设定一个 0.1K 的最小逼近温差容差) + if T_hot_local - T_cold_local < 0.1: + return False # 发生温度交叉，物理不可行！ + + return True + +class Heater(): + """加热器类""" + def __init__(self, name): + self.name = name + self.variables = None + def calculator(self, inlet_state, outlet_state): + h_in = inlet_state['h'] + h_out = outlet_state['h'] + Q_input = h_out - h_in + self.variables = { + 'name': self.name, + 'inlet_state':{ + 'P': inlet_state['P'], # 压强(kPa) + 'T': inlet_state['T'], + 'h': inlet_state['h'], + 's': inlet_state['s'], + }, + 'outlet_state':{ + 'P': outlet_state['P'], # 压强(kPa) + 'T': outlet_state['T'], + 'h': outlet_state['h'], + 's': outlet_state['s'], + }, + 'Q_in': Q_input, + } + +class Condenser(): + """冷凝器类""" + def __init__(self, name): + self.name = name + self.variables = None + def calculator(self, inlet_state, outlet_state): + h_in = inlet_state['h'] + h_out = outlet_state['h'] + Q_output = h_in - h_out + self.variables = { + 'name': self.name, + 'inlet_state':{ + 'P': inlet_state['P'], # 压强(kPa) + 'T': inlet_state['T'], + 'h': inlet_state['h'], + 's': inlet_state['s'], + }, + 'outlet_state':{ + 'P': outlet_state['P'], # 压强(kPa) + 'T': outlet_state['T'], + 'h': outlet_state['h'], + 's': outlet_state['s'], + }, + 'Q_out': Q_output, + } + +class Concentrator(): + """汇流组件""" + def __init__(self, name): + self.name = name + self.variables = None + def calculator(self, inlet_state_bypass, inlet_state_mroad, x, property_calculator): + h_in_bypass = inlet_state_bypass['h'] + h_in_mroad = inlet_state_mroad['h'] + p_in = inlet_state_bypass['P'] + mw = property_calculator.mw + h_out = x * h_in_bypass + (1-x) * h_in_mroad + outlet_state = property_calculator.calculate_properties(P=p_in, h=h_out*mw) + self.variables = { + 'name': self.name, + 'inlet_state_bypass': inlet_state_bypass, + 'inlet_state_mroad': inlet_state_mroad, + 'outlet_state':{ + 'P': p_in, + 'T': outlet_state['T'], + 'h': h_out, + 's': outlet_state['s']/mw + } + } diff --git a/brayton_cycle/cycles.py b/brayton_cycle/cycles.py new file mode 100644 index 0000000..42c6b51 --- /dev/null +++ b/brayton_cycle/cycles.py @@ -0,0 +1,475 @@ +# -*- coding: utf-8 -*- +"""Cycle definitions, efficiency calculations, optimization, and plotting.""" + +from scipy.optimize import minimize_scalar +import numpy as np +import matplotlib.pyplot as plt + +from .components import ( + Compressor, + Concentrator, + Condenser, + Heater, + Recuperator, + Turbine, +) +from .properties import CO2PropertyCalculator + + +class BraytonCycle(): + """循环计算""" + def __init__(self, name, refprop_path = "C:/Program Files (x86)/REFPROP 10.0+/REFPROP"): + self.name = name + self.property_calculator = None + self.compressor = None + self.turbine = None + self.recuperator = None + self.heater = None + self.condenser = None + self.concentrator = None + self.refprop_path = refprop_path + self.property_calculator = CO2PropertyCalculator(self.refprop_path) + + def cycle_eff_calculator(self): + if type(self.compressor) == list: + Wc = (self.compressor[0].variables['Wc'] + + self.compressor[1].variables['Wc']) + else: + Wc = self.compressor.variables['Wc'] + Wt = self.turbine.variables['Wt'] + Q_input = self.heater.variables['Q_in'] + cycle_eff = (Wt - Wc) / Q_input + return cycle_eff + + def SC(self, T_low, T_high, p_low, p_high, param = None): + if param == None: + param = { + 'compressor_eff': 0.98, + 'turbine_eff': 0.95 + } + # 定义循环组件 + self.heater = Heater(name = "Main heater") + self.condenser = Condenser(name = "Main condenser") + self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) + self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) + + # 计算循环参数 + self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) + self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) + self.heater.calculator(self.compressor.variables['outlet_state'], self.turbine.variables['inlet_state']) + self.condenser.calculator(self.turbine.variables['outlet_state'], self.compressor.variables['inlet_state']) + + # 计算循环效率 + cycle_eff = self.cycle_eff_calculator() + return cycle_eff + + def SRC(self, T_low, T_high, p_low, p_high, param=None, ploss=0.0): + if param is None or 'recuperator_eff' not in param: + param = { + 'compressor_eff': 0.98, + 'turbine_eff': 0.95, + 'recuperator_eff': 0.85 + } + + # 定义循环组件 + self.heater = Heater(name = "Main heater") + self.condenser = Condenser(name = "Main condenser") + self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) + self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) + self.recuperator = Recuperator(name = "Recuperator", eff = param['recuperator_eff']) + + # 计算循环参数 + self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) + self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) + self.recuperator.calculator(self.compressor.variables['outlet_state'], + self.turbine.variables['outlet_state'], 0, + ploss, self.property_calculator) + self.condenser.calculator(self.recuperator.variables['hot_outlet_state'], + self.compressor.variables['inlet_state']) + self.heater.calculator(self.recuperator.variables['cold_outlet_state'], + self.turbine.variables['inlet_state']) + + # 计算循环效率 + cycle_eff = self.cycle_eff_calculator() + return cycle_eff + + def RC(self, T_low, T_high, p_low, p_high, ploss, param = None): + if (param is None or + 'recompressor_eff' not in param or + 'highT_recuperator_eff' not in param or + 'x' not in param): + param = { + 'compressor_eff': 0.98, + 'turbine_eff': 0.95, + 'recuperator_eff': 0.85, + 'recompressor_eff': 0.98, + 'highT_recuperator_eff': 0.85, + 'x': 0.9 + } + self.x = param['x'] + # 定义循环组件 + self.heater = Heater(name = "Main heater") + self.condenser = Condenser(name = "Main_condenser") + self.turbine = Turbine(name = "Main Turbine", + eff = param['turbine_eff']) + main_compressor = Compressor(name = "Main compressor", + eff = param['compressor_eff']) + recompressor = Compressor(name = "Recompressor", + eff = param['recompressor_eff']) + lrecuperator = Recuperator(name = "Low Temperature recuprerator", + eff = param['recuperator_eff'], + x = self.x) + hrecuperator = Recuperator(name = "High Temperature recuperator", + eff = param['highT_recuperator_eff']) + self.concentrator = Concentrator(name = "Concentrator") + self.recuperator = [] + self.compressor = [] + + + # 计算组件进出口参数 + # 先算压缩机和涡轮 + main_compressor.calculator(p_low, T_low, p_high, self.property_calculator) + self.turbine.calculator(p_high*(1-ploss)**3, + T_high, p_low*(1-ploss)**(-3), + self.property_calculator) + + # 假定低温回热器出口参数并进行迭代 + mw = self.property_calculator.mw + T_hr_inlet = ((self.turbine.variables['outlet_state']['T'] + + main_compressor.variables['outlet_state']['T'])/2) + # T_lr_inlet = self.turbine.variables['outlet_state']['T'] * 1.01 + max_iter = 300 + relax_fac = 0.4 + for i in range(max_iter): + hr_inlet_state_hot_mol = self.property_calculator.calculate_properties(P=p_high*(1-ploss), T=T_hr_inlet) + # 单位转换成mol + hr_inlet_state_hot = { + 'P': hr_inlet_state_hot_mol['P'], + 'T': hr_inlet_state_hot_mol['T'], + 'h': hr_inlet_state_hot_mol['h']/mw, + 's': hr_inlet_state_hot_mol['s']/mw + } + + hrecuperator.calculator(hr_inlet_state_hot, + self.turbine.variables['outlet_state'], 0, + ploss, self.property_calculator) + lrecuperator.calculator(main_compressor.variables['outlet_state'], + hrecuperator.variables['hot_outlet_state'], 1, + ploss, self.property_calculator) + recompressor.calculator(p_in = p_low/(1-ploss), p_out = p_high*(1-ploss), + T_in = lrecuperator.variables['hot_outlet_state']['T'], + property_calculator = self.property_calculator) + self.concentrator.calculator(recompressor.variables['outlet_state'], + lrecuperator.variables['cold_outlet_state'], + self.x, self.property_calculator) + + Tc_outlet = self.concentrator.variables['outlet_state']['T'] + err = abs(Tc_outlet - T_hr_inlet) + if err <= 1e-5: + break + if i == max_iter - 1: + raise ValueError(f"迭代次数超过范围，当前误差{err:.4f}") + + T_hr_inlet = relax_fac * Tc_outlet + (1 - relax_fac) * T_hr_inlet + + hrecuperator.calculator(self.concentrator.variables['outlet_state'], + self.turbine.variables['outlet_state'], 0, + ploss, self.property_calculator) + lrecuperator.calculator(main_compressor.variables['outlet_state'], + hrecuperator.variables['hot_outlet_state'], 1, + ploss, self.property_calculator) + + main_compressor.variables['Wc'] *= (1-self.x) + recompressor.variables['Wc'] *= self.x + self.recuperator.append(lrecuperator) + self.recuperator.append(hrecuperator) + self.compressor.append(main_compressor) + self.compressor.append(recompressor) + self.condenser.calculator(lrecuperator.variables['hot_outlet_state'], + main_compressor.variables['inlet_state']) + self.condenser.variables['Q_out'] *= (1-self.x) + self.heater.calculator(hrecuperator.variables['cold_outlet_state'], + self.turbine.variables['inlet_state']) + # 计算循环效率 + cycle_eff = self.cycle_eff_calculator() + return cycle_eff + + def base_params_single_optimize(self, fixed_var, base_params, target_var_name, bounds): + """ + 组件性能单变量优化器 + : 固定边界条件 + : 默认参数 + : 优化变量名称 + : 变量范围 + """ + + if target_var_name not in base_params: + raise ValueError(f"参数{target_var_name}不在参数字典中") + + def opt_fun(opt_var): + # 复制变量字典 + opt_param = base_params.copy() + # 修改要优化的变量为参数 + opt_param[target_var_name] = opt_var + # 带入循环参数计算 + try: + eff = self.RC( + T_low = fixed_var['T_low'], + T_high= fixed_var['T_high'], + p_low = fixed_var['p_low'], + p_high = fixed_var['p_high'], + ploss = fixed_var['ploss'], + param = opt_param) + return -eff + except Exception as e: + return 0.0 # 遇到物性计算崩溃时返回极差值 + res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') + if res.success: + print(f"✅ 优化完成！") + print(f"👉 最佳 {target_var_name} = {res.x:.4f}") + print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") + else: + print("❌ 优化失败。") + + return res + + def fixed_params_single_optimize(self, fixed_params, target_var_name, params, bounds): + """ + 边界条件单变量优化器 + : 固定边界条件 + : 默认参数 + : 优化变量名称 + : 变量范围 + """ + if target_var_name not in fixed_params: + raise ValueError(f"参数{target_var_name}不在参数字典中") + + def opt_fun(opt_var): + # 复制变量 + opt_params = fixed_params.copy() + # 变量替换 + opt_params[target_var_name] = opt_var + # 带入循环 + try: + eff = self.RC(T_low = opt_params['T_low'], + T_high = opt_params['T_high'], + p_low = opt_params['p_low'], + p_high = opt_params['p_high'], + ploss = opt_params['ploss'], + param = params) + return -eff + except Exception as e: + return 0.0 + res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') + if res.success: + print(f"✅ 优化完成！") + print(f"👉 最佳 {target_var_name} = {res.x:.4f}") + print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") + else: + print("❌ 优化失败。") + return res + + def plot_optimization_landscape(self, fixed_params, params, target_var_name, bounds, res, num_points=50): + """ + 绘制单变量优化地形图 + :param target_var_name: 要扫描和优化的变量名（如 'x'） + :param bounds: 扫描和优化的范围 (min, max) + :param num_points: 扫描的采样点数量，越大曲线越平滑，但计算越慢 + """ + print(f"开始对【{target_var_name}】进行区间扫描，共计算 {num_points} 个点...") + plt.rcParams['font.sans-serif'] = ['SimHei'] # Windows 用黑体 + plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 + # 1. 生成扫描数组 + x_vals = np.linspace(bounds[0], bounds[1], num_points) + eff_vals = [] + valid_x = [] # 记录那些没有报错的 x + if target_var_name not in fixed_params: + # 2. 遍历计算曲线上的点 + for val in x_vals: + current_param = params.copy() + current_param[target_var_name] = val + try: + # 调用你的黑盒物理模型（注意：这里取正效率用于画图） + eff = self.RC( + T_low=fixed_params['T_low'], + T_high=fixed_params['T_high'], + p_low=fixed_params['p_low'], + p_high=fixed_params['p_high'], + ploss=fixed_params['ploss'], + param=current_param + ) + # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 + # 这里确保只有成功的点才画上去 + eff_vals.append(eff * 100) # 乘以 100 转换为百分比 + valid_x.append(val) + except Exception as e: + # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 + pass + else: + for val in x_vals: + current_param = fixed_params.copy() + current_param[target_var_name] = val + try: + # 调用你的黑盒物理模型（注意：这里取正效率用于画图） + eff = self.RC( + T_low=current_param['T_low'], + T_high=current_param['T_high'], + p_low=current_param['p_low'], + p_high=current_param['p_high'], + ploss=current_param['ploss'], + param=params + ) + # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 + # 这里确保只有成功的点才画上去 + eff_vals.append(eff * 100) # 乘以 100 转换为百分比 + valid_x.append(val) + except Exception as e: + # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 + pass + print("扫描完成！正在使用优化器寻找精确最高点...") + + # 4. 开始绘图 + plt.figure(figsize=(8, 6), dpi=120) # 设置画布大小和清晰度 + + # 画出目标函数曲线 + plt.plot(valid_x, eff_vals, linestyle='-', color='#1f77b4', linewidth=2, label='系统热效率曲线') + + # 如果优化成功，用醒目的红星标出最优点 + if res.success: + best_x = res.x + best_eff = -res.fun * 100 + plt.scatter(best_x, best_eff, color='red', marker='*', s=200, zorder=5, label=f'最优点 ({best_x:.4f}, {best_eff:.2f}%)') + + # 画辅助虚线对齐坐标轴 + plt.axvline(x=best_x, color='gray', linestyle='--', alpha=0.6) + plt.axhline(y=best_eff, color='gray', linestyle='--', alpha=0.6) + else: + print("优化结果未输入!") + # 设置图表装饰 + plt.title(f'系统热效率随 {target_var_name} 的变化趋势', fontsize=14) + plt.xlabel(f'优化变量: {target_var_name}', fontsize=12) + plt.ylabel('循环热效率 η (%)', fontsize=12) + plt.grid(True, linestyle=':', alpha=0.7) + plt.legend(fontsize=11) + + # 显示图像 + plt.tight_layout() + plt.show() + + # 调用测试 + # plot_optimization_landscape('x', bounds=(0.6, 0.95), num_points=40) + + + def sweep_and_optimize(self, fixed_params, params, sweep_var, sweep_bounds, opt_var, opt_bounds, num_points=50): + """ + 带内部动态优化的单变量扫描器 (极度通用版) + + :param fixed_params: 固定的边界条件字典 + :param params: 组件性能参数字典 + :param sweep_var: 你要扫描/遍历的变量名 (例如 'T_low') + :param sweep_bounds: 扫描变量的范围 (min, max) + :param opt_var: 在每个扫描点下，你需要动态寻找最优值的变量名 (例如 'x') + :param opt_bounds: 优化变量的搜索范围 (min, max) + :param num_points: 扫描点数 + """ + print(f"\n🚀 开始执行嵌套扫描：") + print(f" - 扫描变量 (X轴): 【{sweep_var}】 范围 {sweep_bounds}") + print(f" - 内部动态优化变量: 【{opt_var}】 范围 {opt_bounds}") + + # 1. 生成扫描节点 + sweep_vals = np.linspace(sweep_bounds[0], sweep_bounds[1], num_points) + + # 记录数据的列表 + valid_sweep_vals = [] + best_effs = [] + best_opt_vals = [] + + # 2. 开始逐点扫描 + for s_val in sweep_vals: + + # 【核心1：每次必须使用干净的字典副本】 + current_fixed = fixed_params.copy() + current_param = params.copy() + + # 判断扫描变量是属于 fixed_params 还是 params，并赋值 + if sweep_var in current_fixed: + current_fixed[sweep_var] = s_val + elif sweep_var in current_param: + current_param[sweep_var] = s_val + else: + raise ValueError(f"找不到扫描变量: {sweep_var}") + + # 3. 定义内部优化目标函数 (闭包) + def inner_objective(guess_val): + # 将优化器猜的值赋给 opt_var + current_param[opt_var] = guess_val + + try: + # 调用黑盒计算 + eff = self.RC( + T_low=current_fixed['T_low'], + T_high=current_fixed['T_high'], + p_low=current_fixed['p_low'], + p_high=current_fixed['p_high'], + ploss=current_fixed['ploss'], + param=current_param + ) + # 【预留口：此处可加入换热器内部夹点校验】 + # 取出两个换热器进行夹点校验 + for rec in self.recuperator: + if not rec.check_pinch_point(self.property_calculator): + return 0.0 # 核心！如果交叉了，直接返回 0 效率，强迫优化器换参数 + return -eff + except Exception: + return 0.0 + + # 4. 调用一维优化器 + res = minimize_scalar(inner_objective, bounds=opt_bounds, method='bounded') + + # 5. 结果校验与存储 + if res.success and -res.fun > 0: + best_eff = -res.fun * 100 + best_opt_val = res.x + + valid_sweep_vals.append(s_val) + best_effs.append(best_eff) + best_opt_vals.append(best_opt_val) + + print(f"✔️ {sweep_var} = {s_val:.2f} | 寻得最优 {opt_var} = {best_opt_val:.4f} | 最高效率 = {best_eff:.2f}%") + else: + print(f"❌ {sweep_var} = {s_val:.2f} | 优化失败或物理无解，已跳过") + + # 6. 调用画图方法 (将画图剥离，保持代码干净) + self._plot_results(sweep_var, valid_sweep_vals, best_effs, opt_var, best_opt_vals) + + return valid_sweep_vals, best_effs, best_opt_vals + + def _plot_results(self, sweep_var, x_data, y_eff_data, opt_var, y_opt_data): + """专门用来画图的内部方法，支持双Y轴""" + if not x_data: + print("没有有效数据可供绘制！") + return + + plt.rcParams['font.sans-serif'] = ['SimHei'] + plt.rcParams['axes.unicode_minus'] = False + + fig, ax1 = plt.subplots(figsize=(9, 6), dpi=120) + + # 画左Y轴：最高效率曲线 + color1 = '#1f77b4' + ax1.set_xlabel(f'扫描变量: {sweep_var}', fontsize=12) + ax1.set_ylabel('最优循环热效率 η (%)', color=color1, fontsize=12) + ax1.plot(x_data, y_eff_data, color=color1, linewidth=2.5, label='系统热效率') + ax1.tick_params(axis='y', labelcolor=color1) + ax1.grid(True, linestyle=':', alpha=0.6) + + # 画右Y轴：对应的最优分流量走势 + ax2 = ax1.twinx() + color2 = '#d62728' + ax2.set_ylabel(f'匹配的最优动态变量: {opt_var}', color=color2, fontsize=12) + ax2.plot(x_data, y_opt_data, color=color2, linestyle='--', linewidth=2, label=f'最优 {opt_var} 值') + ax2.tick_params(axis='y', labelcolor=color2) + + plt.title(f'系统最高效率及对应的最优 {opt_var} 随 {sweep_var} 的变化', fontsize=14) + fig.tight_layout() + plt.show() diff --git a/brayton_cycle/properties.py b/brayton_cycle/properties.py new file mode 100644 index 0000000..231487b --- /dev/null +++ b/brayton_cycle/properties.py @@ -0,0 +1,65 @@ +# -*- coding: utf-8 -*- +"""REFPROP-backed CO2 property helpers.""" + +import ctREFPROP.ctREFPROP as ct + + +class CO2PropertyCalculator(): + """二氧化碳物性计算""" + def __init__(self, refprop_path = None): + """初始化库""" + self.rp = ct.REFPROPFunctionLibrary(refprop_path) + # 设置流体文件 + self.rp.SETUPdll(1, 'CO2.FLD', 'HMX.BNC', 'DEF') + # 设置单位 + 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): + """计算二氧化碳物性""" + if T is not None and P is not None: + # 已知Tp + result = self.rp.TPFLSHdll(T, P, self.z) + properties = { + 'T': T, + 'P': P, + 'h': result.h, + 's': result.s, + } + elif P is not None and h is not None: + # 已知Ph + result = self.rp.PHFLSHdll(P, h, self.z) + properties = { + 'T': result.T, + 'P': P, + 'h': h, + 's': result.s, + } + elif T is not None and h is not None: + # 已知Th + result = self.rp.THFLSHdll(T, h, self.z) + properties = { + 'T': T, + 'P': result.P, + 'h': h, + 's': result.s, + } + elif P is not None and s is not None: + # 已知Ps + result = self.rp.PSFLSHdll(P, s, self.z) + properties = { + 'T': result.T, + 'P': P, + 'h': result.h, + 's': s, + } + else: + raise ValueError("提供的参数不足") + if result.ierr > 0: + raise ValueError(f"REFPROP计算错误:{result.ierr}") + + # 补充提取的物性，这里由于后续还要使用，不进行参数变换 + properties['D'] = result.D + properties['cp'] = result.Cp + properties['cv'] = result.Cv, + return properties diff --git a/brayton_cycle_test.py b/brayton_cycle_test.py index 9de7a8c..b3e0d60 100644 --- a/brayton_cycle_test.py +++ b/brayton_cycle_test.py @@ -1,980 +1,59 @@ # -*- coding: utf-8 -*- -""" -Created on Mon Dec 22 15:35:07 2025 +"""Run an example recompression CO2 Brayton cycle sweep.""" -""" -import ctREFPROP.ctREFPROP as ct -from scipy.optimize import minimize_scalar -import numpy as np -import matplotlib.pyplot as plt +from brayton_cycle import BraytonCycle -class CO2PropertyCalculator(): - """二氧化碳物性计算""" - def __init__(self, refprop_path = None): - """初始化库""" - self.rp = ct.REFPROPFunctionLibrary(refprop_path) - # 设置流体文件 - self.rp.SETUPdll(1, 'CO2.FLD', 'HMX.BNC', 'DEF') - # 设置单位 - 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): - """计算二氧化碳物性""" - if T is not None and P is not None: - # 已知Tp - result = self.rp.TPFLSHdll(T, P, self.z) - properties = { - 'T': T, - 'P': P, - 'h': result.h, - 's': result.s, - } - elif P is not None and h is not None: - # 已知Ph - result = self.rp.PHFLSHdll(P, h, self.z) - properties = { - 'T': result.T, - 'P': P, - 'h': h, - 's': result.s, - } - elif T is not None and h is not None: - # 已知Th - result = self.rp.THFLSHdll(T, h, self.z) - properties = { - 'T': T, - 'P': result.P, - 'h': h, - 's': result.s, - } - elif P is not None and s is not None: - # 已知Ps - result = self.rp.PSFLSHdll(P, s, self.z) - properties = { - 'T': result.T, - 'P': P, - 'h': result.h, - 's': s, - } - else: - raise ValueError("提供的参数不足") - if result.ierr > 0: - raise ValueError(f"REFPROP计算错误:{result.ierr}") - - # 补充提取的物性，这里由于后续还要使用，不进行参数变换 - properties['D'] = result.D - properties['cp'] = result.Cp - properties['cv'] = result.Cv, - return properties - -class Compressor(): - """压缩机类""" - def __init__(self, name, eff): - """ - 初始化参数 - - name: 名称 - eff: 等熵效率 - Wc: 压缩功 - inlet_state: 入口参数 - outlet_state: 出口参数 - outlet_state_is: 等熵状态下出口参数 - """ - - self.name = name - self.eff = eff - self.variables = None - - def calculator(self, p_in, T_in, p_out, property_calculator): - # 先计算熵值 - inlet_state = property_calculator.calculate_properties(T=T_in, P=p_in) - mw = property_calculator.mw - s = inlet_state['s'] - h_in = inlet_state['h'] - outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s) - h_out_is = outlet_state_is['h'] - h_out = h_in + (h_out_is - h_in) / self.eff - outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out) - Wc = h_out - h_in - T_out = outlet_state['T'] - # 计算结果 - self.variables = { - 'name': self.name, - 'inlet_state':{ - 'P': p_in, # 压强(kPa) - 'T': T_in, - 'h': h_in/mw, # 比焓(J/mol)->(J/kg) - 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) - }, - 'outlet_state':{ - 'P': p_out, # 压强(kPa) - 'T': T_out, - 'h': h_out/mw, # 比焓(kJ/mol)->(kJ/kg) - 's': s/mw, # 比熵(kJ/mol.K)->(kJ/kg.K) - }, - 'eff': self.eff, - 'Wc': Wc/mw, # 压缩功(kJ/mol)->(kJ/kg) - 'pi': p_out / p_in - } -class Turbine(): - """透平类""" - def __init__(self, name, eff): - """ - 初始化参数 - name: 名称 - eff: 透平效率 - """ - self.name = name - self.eff = eff - self.variables = None - - def calculator(self, p_in, T_in, p_out, property_calculator): - """涡轮参数计算""" - inlet_state = property_calculator.calculate_properties(P=p_in, T=T_in) - mw = property_calculator.mw - s = inlet_state['s'] - h_in = inlet_state['h'] - outlet_state_is = property_calculator.calculate_properties(P=p_out, s=s) - h_out_is = outlet_state_is['h'] - h_out = h_in - (h_in - h_out_is) * self.eff - outlet_state = property_calculator.calculate_properties(P=p_out, h=h_out) - T_out = outlet_state['T'] - Wt = h_in - h_out - self.variables = { - 'name': self.name, - 'inlet_state':{ - 'P': p_in, # 压强(kPa) - 'T': T_in, - 'h': h_in/mw, # 比焓(J/mol)->(J/kg) - 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) - }, - 'outlet_state':{ - 'P': p_out, # 压强(kPa) - 'T': T_out, - 'h': h_out/mw, # 比焓(J/mol)->(J/kg) - 's': s/mw, # 比熵(J/mol.K)->(J/kg.K) - }, - 'eff': self.eff, - 'Wt': Wt/mw, # 透平做功(J/mol)->(J/kg) - 'pi': p_in / p_out - } - -class Recuperator(): - """换热器类""" - def __init__(self, name, eff, x=0): - """ - 初始化参数 - name: 名称 - eff: 换热效率(基于焓的计算方法) - """ - self.name = name - self.eff = eff - self.Q_ex = None - self.variables = None - self.x = x - - def calculator(self, cold_inlet_state, hot_inlet_state, bypass_info, ploss, property_calculator): - """ - 计算换热器两侧参数，默认逆流 - bypass_info: 是否存在分流，0为不存在，1为存在 - - 下标含义 - ---------- - 1: 换热器冷端入口 - 2: 换热器冷端出口 - 3: 换热器热端入口 - 4: 换热器热端出口 - ass: 迭代中间变量，假设值 - """ - - # 计算两入口参数, 这里单位是kg - p1 = cold_inlet_state['P'] - T1 = cold_inlet_state['T'] - h1 = cold_inlet_state['h'] - p2 = p1 * (1 - ploss) - p3 = hot_inlet_state['P'] - T3 = hot_inlet_state['T'] - h3 = hot_inlet_state['h'] - p4 = p3 * (1 - ploss) - mw = property_calculator.mw - - # 设定质量流量 - m_cold = (1-self.x) if bypass_info == 1 else 1.0 - m_hot = 1.0 - -# ============================================================================= -# # 采用焓差效能的方式来计算换热器进出口参数 -# # 假设最大温差发生在冷端, 计算冷端出口温度和比焓 -# cold_outlet_state = property_calculator.calculate_properties(T=T3, P=p2) -# h2 = cold_outlet_state['h'] / mw -# Q_ass_cold = m_cold * abs(h2 - h1) -# -# # 假设最大温差发生在热端, 计算热端出口温度和比焓 -# hot_outlet_state = property_calculator.calculate_properties(T=T1, P=p4) -# h4 = hot_outlet_state['h'] / mw -# Q_ass_hot = m_hot * abs(h3 - h4) -# -# # 比较两个可能的Q，取最小值与焓差效能的乘积作为实际换热量 -# self.Q_ex = min(Q_ass_hot, Q_ass_cold) * self.eff -# # 由实际换热量计算出口焓和出口状态 -# h2 = h1 + self.Q_ex / m_cold -# h4 = h3 - self.Q_ex / m_hot -# cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw) -# hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw) -# T2 = cold_outlet_state['T'] -# T4 = hot_outlet_state['T'] -# ============================================================================= +def run_rc_sweep(): + brayton = BraytonCycle(name="simple brayton cycle test") - # 采用温差效能的方式来计算换热器进出口参数 - T_ass_max = abs(T1 - T3) - # 假设最大温差发生在冷端, 计算冷端出口温度和比焓 - T2 = T1 + T_ass_max * self.eff - cold_outlet_state = property_calculator.calculate_properties(T=T2, P=p2) - h2 = cold_outlet_state['h'] / mw - Q_ass_cold = m_cold * abs(h2 - h1) - - # 假设最大温差发生在热端, 计算热端出口温度和比焓 - T4 = T3 - T_ass_max * self.eff - hot_outlet_state = property_calculator.calculate_properties(T=T4, P=p4) - h4 = hot_outlet_state['h'] / mw - Q_ass_hot = m_hot * abs(h3 - h4) - - # 比较两个可能的Q，取最小值与焓差效能的乘积作为实际换热量 - self.Q_ex = min(Q_ass_hot, Q_ass_cold) - # 由实际换热量计算出口焓和出口状态 - h2 = h1 + self.Q_ex / m_cold - h4 = h3 - self.Q_ex / m_hot - cold_outlet_state = property_calculator.calculate_properties(P=p2, h=h2*mw) - hot_outlet_state = property_calculator.calculate_properties(P=p4, h=h4*mw) - T2 = cold_outlet_state['T'] - T4 = hot_outlet_state['T'] - - # 拼装变量 - self.variables = { - 'name': self.name, - 'cold_inlet_state':{ - 'P': p1, - 'T': T1, - 'h': h1, - 's': cold_inlet_state['s'] - }, - 'cold_outlet_state':{ - 'P': p2, - 'T': cold_outlet_state['T'], - 'h': h2, - 's': cold_outlet_state['s']/mw - }, - 'hot_inlet_state':{ - 'P': p3, - 'T': T3, - 'h': h3, - 's': hot_inlet_state['s'] - }, - 'hot_outlet_state':{ - 'P': p4, - 'T': T4, - 'h': h4, - 's': hot_outlet_state['s']/mw - }, - 'eff': self.eff, - 'Q_exchange': self.Q_ex - } - def check_pinch_point(self, property_calculator, num_segments=20): - """ - 换热器内部夹点校验 - 将换热量均分为 num_segments 段，检查内部每个微元的冷热流体温度 - """ - h_cold_in = self.variables['cold_inlet_state']['h'] - h_hot_in = self.variables['hot_inlet_state']['h'] - p_cold = self.variables['cold_inlet_state']['P'] - p_hot = self.variables['hot_inlet_state']['P'] - - m_cold = (1 - self.x) if self.name == "Low Temperature recuprerator" else 1.0 - m_hot = 1.0 - - dQ = self.Q_ex / num_segments - - # 沿冷流体流动方向步进检查 - for i in range(num_segments + 1): - q_current = i * dQ - - # 当前微元截面的焓值 - h_cold_local = h_cold_in + q_current / m_cold - h_hot_local = (h_hot_in - self.Q_ex / m_hot) + q_current / m_hot - - # 查温度 - T_cold_local = property_calculator.calculate_properties(P=p_cold, h=h_cold_local * property_calculator.mw)['T'] - T_hot_local = property_calculator.calculate_properties(P=p_hot, h=h_hot_local * property_calculator.mw)['T'] - - # 如果热流体温度低于等于冷流体温度 (设定一个 0.1K 的最小逼近温差容差) - if T_hot_local - T_cold_local < 0.1: - return False # 发生温度交叉，物理不可行！ - - return True # 校验通过 - -class Heater(): - """加热器类""" - def __init__(self, name): - self.name = name - self.variables = None - def calculator(self, inlet_state, outlet_state): - h_in = inlet_state['h'] - h_out = outlet_state['h'] - Q_input = h_out - h_in - self.variables = { - 'name': self.name, - 'inlet_state':{ - 'P': inlet_state['P'], # 压强(kPa) - 'T': inlet_state['T'], - 'h': inlet_state['h'], - 's': inlet_state['s'], - }, - 'outlet_state':{ - 'P': outlet_state['P'], # 压强(kPa) - 'T': outlet_state['T'], - 'h': outlet_state['h'], - 's': outlet_state['s'], - }, - 'Q_in': Q_input, - } - -class Condenser(): - """冷凝器类""" - def __init__(self, name): - self.name = name - self.variables = None - def calculator(self, inlet_state, outlet_state): - h_in = inlet_state['h'] - h_out = outlet_state['h'] - Q_output = h_in - h_out - self.variables = { - 'name': self.name, - 'inlet_state':{ - 'P': inlet_state['P'], # 压强(kPa) - 'T': inlet_state['T'], - 'h': inlet_state['h'], - 's': inlet_state['s'], - }, - 'outlet_state':{ - 'P': outlet_state['P'], # 压强(kPa) - 'T': outlet_state['T'], - 'h': outlet_state['h'], - 's': outlet_state['s'], - }, - 'Q_out': Q_output, - } -class Concentrator(): - """汇流组件""" - def __init__(self, name): - self.name = name - self.variables = None - def calculator(self, inlet_state_bypass, inlet_state_mroad, x, property_calculator): - h_in_bypass = inlet_state_bypass['h'] - h_in_mroad = inlet_state_mroad['h'] - p_in = inlet_state_bypass['P'] - mw = property_calculator.mw - h_out = x * h_in_bypass + (1-x) * h_in_mroad - outlet_state = property_calculator.calculate_properties(P=p_in, h=h_out*mw) - self.variables = { - 'name': self.name, - 'inlet_state_bypass': inlet_state_bypass, - 'inlet_state_mroad': inlet_state_mroad, - 'outlet_state':{ - 'P': p_in, - 'T': outlet_state['T'], - 'h': h_out, - 's': outlet_state['s']/mw - } - } -class BraytonCycle(): - """循环计算""" - def __init__(self, name, refprop_path = "C:/Program Files (x86)/REFPROP 10.0+/REFPROP"): - self.name = name - self.property_calculator = None - self.compressor = None - self.turbine = None - self.recuperator = None - self.heater = None - self.condenser = None - self.concentrator = None - self.refprop_path = refprop_path - self.property_calculator = CO2PropertyCalculator(self.refprop_path) - - def cycle_eff_calculator(self): - if type(self.compressor) == list: - Wc = (self.compressor[0].variables['Wc'] + - self.compressor[1].variables['Wc']) - else: - Wc = self.compressor.variables['Wc'] - Wt = self.turbine.variables['Wt'] - Q_input = self.heater.variables['Q_in'] - cycle_eff = (Wt - Wc) / Q_input - return cycle_eff - - def SC(self, T_low, T_high, p_low, p_high, param = None): - if param == None: - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95 - } - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main condenser") - self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) - self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) - - # 计算循环参数 - self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) - self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) - self.heater.calculator(self.compressor.variables['outlet_state'], self.turbine.variables['inlet_state']) - self.condenser.calculator(self.turbine.variables['outlet_state'], self.compressor.variables['inlet_state']) - - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff - - def SRC(self, T_low, T_high, p_low, p_high, param = None): - if ('recuperator_eff' not in param) or param == None: - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95, - 'recuperator_eff': 0.85 - } - - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main condenser") - self.compressor = Compressor(name = "Main compressor", eff = param['compressor_eff']) - self.turbine = Turbine(name = "Turbine", eff = param['turbine_eff']) - self.recuperator = Recuperator(name = "Recuperator", eff = param['recuperator_eff']) - - # 计算循环参数 - self.compressor.calculator(p_low, T_low, p_high, self.property_calculator) - self.turbine.calculator(p_high, T_high, p_low, self.property_calculator) - self.recuperator.calculator(self.compressor.variables['outlet_state'], - self.turbine.variables['outlet_state'], 0, - self.property_calculator) - self.condenser.calculator(self.recuperator.variables['hot_outlet_state'], - self.compressor.variables['inlet_state']) - self.heater.calculator(self.recuperator.variables['cold_outlet_state'], - self.turbine.variables['inlet_state']) - - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff - - def RC(self, T_low, T_high, p_low, p_high, ploss, param = None): - if ('recompressor_eff' not in param or - 'highT_recuperator_eff' not in param or - 'x' not in param or - param == None): - param = { - 'compressor_eff': 0.98, - 'turbine_eff': 0.95, - 'recuperator_eff': 0.85, - 'recompressor_eff': 0.98, - 'highT_recuperator_eff': 0.85, - 'x': 0.9 - } - self.x = param['x'] - # 定义循环组件 - self.heater = Heater(name = "Main heater") - self.condenser = Condenser(name = "Main_condenser") - self.turbine = Turbine(name = "Main Turbine", - eff = param['turbine_eff']) - main_compressor = Compressor(name = "Main compressor", - eff = param['compressor_eff']) - recompressor = Compressor(name = "Recompressor", - eff = param['recompressor_eff']) - lrecuperator = Recuperator(name = "Low Temperature recuprerator", - eff = param['recuperator_eff'], - x = self.x) - hrecuperator = Recuperator(name = "High Temperature recuperator", - eff = param['highT_recuperator_eff']) - self.concentrator = Concentrator(name = "Concentrator") - self.recuperator = [] - self.compressor = [] - - - # 计算组件进出口参数 - # 先算压缩机和涡轮 - main_compressor.calculator(p_low, T_low, p_high, self.property_calculator) - self.turbine.calculator(p_high*(1-ploss)**3, - T_high, p_low*(1-ploss)**(-3), - self.property_calculator) - - # 假定低温回热器出口参数并进行迭代 - mw = self.property_calculator.mw - T_hr_inlet = ((self.turbine.variables['outlet_state']['T'] - + main_compressor.variables['outlet_state']['T'])/2) - # T_lr_inlet = self.turbine.variables['outlet_state']['T'] * 1.01 - max_iter = 300 - relax_fac = 0.4 - for i in range(max_iter): - hr_inlet_state_hot_mol = self.property_calculator.calculate_properties(P=p_high*(1-ploss), T=T_hr_inlet) - # 单位转换成mol - hr_inlet_state_hot = { - 'P': hr_inlet_state_hot_mol['P'], - 'T': hr_inlet_state_hot_mol['T'], - 'h': hr_inlet_state_hot_mol['h']/mw, - 's': hr_inlet_state_hot_mol['s']/mw - } - - hrecuperator.calculator(hr_inlet_state_hot, - self.turbine.variables['outlet_state'], 0, - ploss, self.property_calculator) - lrecuperator.calculator(main_compressor.variables['outlet_state'], - hrecuperator.variables['hot_outlet_state'], 1, - ploss, self.property_calculator) - recompressor.calculator(p_in = p_low/(1-ploss), p_out = p_high*(1-ploss), - T_in = lrecuperator.variables['hot_outlet_state']['T'], - property_calculator = self.property_calculator) - self.concentrator.calculator(recompressor.variables['outlet_state'], - lrecuperator.variables['cold_outlet_state'], - self.x, self.property_calculator) - - Tc_outlet = self.concentrator.variables['outlet_state']['T'] - err = abs(Tc_outlet - T_hr_inlet) - if err <= 1e-5: - break - if i == max_iter - 1: - raise ValueError(f"迭代次数超过范围，当前误差{err:.4f}") - - T_hr_inlet = relax_fac * Tc_outlet + (1 - relax_fac) * T_hr_inlet - - hrecuperator.calculator(self.concentrator.variables['outlet_state'], - self.turbine.variables['outlet_state'], 0, - ploss, self.property_calculator) - lrecuperator.calculator(main_compressor.variables['outlet_state'], - hrecuperator.variables['hot_outlet_state'], 1, - ploss, self.property_calculator) - - main_compressor.variables['Wc'] *= (1-self.x) - recompressor.variables['Wc'] *= self.x - self.recuperator.append(lrecuperator) - self.recuperator.append(hrecuperator) - self.compressor.append(main_compressor) - self.compressor.append(recompressor) - self.condenser.calculator(lrecuperator.variables['hot_outlet_state'], - main_compressor.variables['inlet_state']) - self.condenser.variables['Q_out'] *= (1-self.x) - self.heater.calculator(hrecuperator.variables['cold_outlet_state'], - self.turbine.variables['inlet_state']) - # 计算循环效率 - cycle_eff = self.cycle_eff_calculator() - return cycle_eff - - def base_params_single_optimize(self, fixed_var, base_params, target_var_name, bounds): - """ - 组件性能单变量优化器 - : 固定边界条件 - : 默认参数 - : 优化变量名称 - : 变量范围 - """ - - if target_var_name not in base_params: - raise ValueError(f"参数{target_var_name}不在参数字典中") - - def opt_fun(opt_var): - # 复制变量字典 - opt_param = base_params.copy() - # 修改要优化的变量为参数 - opt_param[target_var_name] = opt_var - # 带入循环参数计算 - try: - eff = self.RC( - T_low = fixed_var['T_low'], - T_high= fixed_var['T_high'], - p_low = fixed_var['p_low'], - p_high = fixed_var['p_high'], - ploss = fixed_var['ploss'], - param = opt_param) - return -eff - except Exception as e: - return 0.0 # 遇到物性计算崩溃时返回极差值 - res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') - if res.success: - print(f"✅ 优化完成！") - print(f"👉 最佳 {target_var_name} = {res.x:.4f}") - print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") - else: - print("❌ 优化失败。") - - return res - - def fixed_params_single_optimize(self, fixed_params, target_var_name, params, bounds): - """ - 边界条件单变量优化器 - : 固定边界条件 - : 默认参数 - : 优化变量名称 - : 变量范围 - """ - if target_var_name not in fixed_params: - raise ValueError(f"参数{target_var_name}不在参数字典中") - - def opt_fun(opt_var): - # 复制变量 - opt_params = fixed_params.copy() - # 变量替换 - opt_params[target_var_name] = opt_var - # 带入循环 - try: - eff = self.RC(T_low = opt_params['T_low'], - T_high = opt_params['T_high'], - p_low = opt_params['p_low'], - p_high = opt_params['p_high'], - ploss = opt_params['ploss'], - param = params) - return -eff - except Exception as e: - return 0.0 - res = minimize_scalar(opt_fun, bounds=bounds, method='bounded') - if res.success: - print(f"✅ 优化完成！") - print(f"👉 最佳 {target_var_name} = {res.x:.4f}") - print(f"👉 此时系统最高效率 = {-res.fun:.2%}\n") - else: - print("❌ 优化失败。") - return res - - def plot_optimization_landscape(self, fixed_params, params, target_var_name, bounds, res, num_points=50): - """ - 绘制单变量优化地形图 - :param target_var_name: 要扫描和优化的变量名（如 'x'） - :param bounds: 扫描和优化的范围 (min, max) - :param num_points: 扫描的采样点数量，越大曲线越平滑，但计算越慢 - """ - print(f"开始对【{target_var_name}】进行区间扫描，共计算 {num_points} 个点...") - plt.rcParams['font.sans-serif'] = ['SimHei'] # Windows 用黑体 - plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 - # 1. 生成扫描数组 - x_vals = np.linspace(bounds[0], bounds[1], num_points) - eff_vals = [] - valid_x = [] # 记录那些没有报错的 x - if target_var_name not in fixed_params: - # 2. 遍历计算曲线上的点 - for val in x_vals: - current_param = params.copy() - current_param[target_var_name] = val - try: - # 调用你的黑盒物理模型（注意：这里取正效率用于画图） - eff = self.RC( - T_low=fixed_params['T_low'], - T_high=fixed_params['T_high'], - p_low=fixed_params['p_low'], - p_high=fixed_params['p_high'], - ploss=fixed_params['ploss'], - param=current_param - ) - # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 - # 这里确保只有成功的点才画上去 - eff_vals.append(eff * 100) # 乘以 100 转换为百分比 - valid_x.append(val) - except Exception as e: - # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 - pass - else: - for val in x_vals: - current_param = fixed_params.copy() - current_param[target_var_name] = val - try: - # 调用你的黑盒物理模型（注意：这里取正效率用于画图） - eff = self.RC( - T_low=current_param['T_low'], - T_high=current_param['T_high'], - p_low=current_param['p_low'], - p_high=current_param['p_high'], - ploss=current_param['ploss'], - param=params - ) - # 如果系统加了夹点校验且没通过，可能会返回 None 或者抛异常 - # 这里确保只有成功的点才画上去 - eff_vals.append(eff * 100) # 乘以 100 转换为百分比 - valid_x.append(val) - except Exception as e: - # 如果某个 x 导致计算崩溃，我们跳过这个点，不画它 - pass - print("扫描完成！正在使用优化器寻找精确最高点...") - - # 4. 开始绘图 - plt.figure(figsize=(8, 6), dpi=120) # 设置画布大小和清晰度 - - # 画出目标函数曲线 - plt.plot(valid_x, eff_vals, linestyle='-', color='#1f77b4', linewidth=2, label='系统热效率曲线') - - # 如果优化成功，用醒目的红星标出最优点 - if res.success: - best_x = res.x - best_eff = -res.fun * 100 - plt.scatter(best_x, best_eff, color='red', marker='*', s=200, zorder=5, label=f'最优点 ({best_x:.4f}, {best_eff:.2f}%)') - - # 画辅助虚线对齐坐标轴 - plt.axvline(x=best_x, color='gray', linestyle='--', alpha=0.6) - plt.axhline(y=best_eff, color='gray', linestyle='--', alpha=0.6) - else: - print("优化结果未输入!") - # 设置图表装饰 - plt.title(f'系统热效率随 {target_var_name} 的变化趋势', fontsize=14) - plt.xlabel(f'优化变量: {target_var_name}', fontsize=12) - plt.ylabel('循环热效率 η (%)', fontsize=12) - plt.grid(True, linestyle=':', alpha=0.7) - plt.legend(fontsize=11) - - # 显示图像 - plt.tight_layout() - plt.show() - - # 调用测试 - # plot_optimization_landscape('x', bounds=(0.6, 0.95), num_points=40) - - - def sweep_and_optimize(self, fixed_params, params, sweep_var, sweep_bounds, opt_var, opt_bounds, num_points=50): - """ - 带内部动态优化的单变量扫描器 (极度通用版) - - :param fixed_params: 固定的边界条件字典 - :param params: 组件性能参数字典 - :param sweep_var: 你要扫描/遍历的变量名 (例如 'T_low') - :param sweep_bounds: 扫描变量的范围 (min, max) - :param opt_var: 在每个扫描点下，你需要动态寻找最优值的变量名 (例如 'x') - :param opt_bounds: 优化变量的搜索范围 (min, max) - :param num_points: 扫描点数 - """ - print(f"\n🚀 开始执行嵌套扫描：") - print(f" - 扫描变量 (X轴): 【{sweep_var}】 范围 {sweep_bounds}") - print(f" - 内部动态优化变量: 【{opt_var}】 范围 {opt_bounds}") - - # 1. 生成扫描节点 - sweep_vals = np.linspace(sweep_bounds[0], sweep_bounds[1], num_points) - - # 记录数据的列表 - valid_sweep_vals = [] - best_effs = [] - best_opt_vals = [] - - # 2. 开始逐点扫描 - for s_val in sweep_vals: - - # 【核心1：每次必须使用干净的字典副本】 - current_fixed = fixed_params.copy() - current_param = params.copy() - - # 判断扫描变量是属于 fixed_params 还是 params，并赋值 - if sweep_var in current_fixed: - current_fixed[sweep_var] = s_val - elif sweep_var in current_param: - current_param[sweep_var] = s_val - else: - raise ValueError(f"找不到扫描变量: {sweep_var}") - - # 3. 定义内部优化目标函数 (闭包) - def inner_objective(guess_val): - # 将优化器猜的值赋给 opt_var - current_param[opt_var] = guess_val - - try: - # 调用黑盒计算 - eff = self.RC( - T_low=current_fixed['T_low'], - T_high=current_fixed['T_high'], - p_low=current_fixed['p_low'], - p_high=current_fixed['p_high'], - ploss=current_fixed['ploss'], - param=current_param - ) - # 【预留口：此处可加入换热器内部夹点校验】 - # 取出两个换热器进行夹点校验 - for rec in self.recuperator: - if not rec.check_pinch_point(self.property_calculator): - return 0.0 # 核心！如果交叉了，直接返回 0 效率，强迫优化器换参数 - return -eff - except Exception: - return 0.0 - - # 4. 调用一维优化器 - res = minimize_scalar(inner_objective, bounds=opt_bounds, method='bounded') - - # 5. 结果校验与存储 - if res.success and -res.fun > 0: - best_eff = -res.fun * 100 - best_opt_val = res.x - - valid_sweep_vals.append(s_val) - best_effs.append(best_eff) - best_opt_vals.append(best_opt_val) - - print(f"✔️ {sweep_var} = {s_val:.2f} | 寻得最优 {opt_var} = {best_opt_val:.4f} | 最高效率 = {best_eff:.2f}%") - else: - print(f"❌ {sweep_var} = {s_val:.2f} | 优化失败或物理无解，已跳过") - - # 6. 调用画图方法 (将画图剥离，保持代码干净) - self._plot_results(sweep_var, valid_sweep_vals, best_effs, opt_var, best_opt_vals) - - return valid_sweep_vals, best_effs, best_opt_vals - - def _plot_results(self, sweep_var, x_data, y_eff_data, opt_var, y_opt_data): - """专门用来画图的内部方法，支持双Y轴""" - if not x_data: - print("没有有效数据可供绘制！") - return - - plt.rcParams['font.sans-serif'] = ['SimHei'] - plt.rcParams['axes.unicode_minus'] = False - - fig, ax1 = plt.subplots(figsize=(9, 6), dpi=120) - - # 画左Y轴：最高效率曲线 - color1 = '#1f77b4' - ax1.set_xlabel(f'扫描变量: {sweep_var}', fontsize=12) - ax1.set_ylabel('最优循环热效率 η (%)', color=color1, fontsize=12) - ax1.plot(x_data, y_eff_data, color=color1, linewidth=2.5, label='系统热效率') - ax1.tick_params(axis='y', labelcolor=color1) - ax1.grid(True, linestyle=':', alpha=0.6) - - # 画右Y轴：对应的最优分流量走势 - ax2 = ax1.twinx() - color2 = '#d62728' - ax2.set_ylabel(f'匹配的最优动态变量: {opt_var}', color=color2, fontsize=12) - ax2.plot(x_data, y_opt_data, color=color2, linestyle='--', linewidth=2, label=f'最优 {opt_var} 值') - ax2.tick_params(axis='y', labelcolor=color2) - - plt.title(f'系统最高效率及对应的最优 {opt_var} 随 {sweep_var} 的变化', fontsize=14) - fig.tight_layout() - plt.show() - -if __name__ == "__main__": - brayton1 = BraytonCycle(name = "simple brayton cycle test") - T_high = 650+273.15 - T_low = 42+273.15 + T_high = 650 + 273.15 + T_low = 42 + 273.15 p_high = 20.0e6 p_low = 9.09e6 - fixed_var = { - 'T_high': 650+273.15, - 'T_low': 42+273.15, - 'p_high': 20.0e3, - 'p_low': 8.16e3, - 'ploss': 0.01 - } - # param = { - # 'compressor_eff': 0.9, - # 'turbine_eff': 0.93, - # 'recuprator_eff': 0.96 - # } - - # ========================================================================= - # 焓差效能参数定义 - # param = { - # 'x': 0.279, - # 'compressor_eff': 0.9, - # 'recompressor_eff': 0.855, - # 'turbine_eff': 0.93, - # 'recuperator_eff': 0.9354, - # 'highT_recuperator_eff': 0.8857 - # } - # ========================================================================= - # eff = brayton1.simple_brayton_cycle(T_low, T_high, p_low/1e3, p_high/1e3, param) - # - param = { - 'x': 0.279, - 'compressor_eff': 0.9, - 'recompressor_eff': 0.9, - 'turbine_eff': 0.93, - 'recuperator_eff': 0.94, - 'highT_recuperator_eff': 0.96 - } ploss = 0.01 - eff = brayton1.RC(T_low, T_high, p_low/1e3, p_high/1e3, ploss, param) - c = brayton1.compressor - h = brayton1.heater - cond = brayton1.condenser - conc = brayton1.concentrator.variables - t = brayton1.turbine - r = brayton1.recuperator - vl = brayton1.recuperator[0].variables - vh = brayton1.recuperator[1].variables - Wc = brayton1.compressor[0].variables['Wc'] + brayton1.compressor[1].variables['Wc'] - # Wc = brayton1.compressor.variables['Wc'] - Wt = brayton1.turbine.variables['Wt'] - Q_input = brayton1.heater.variables['Q_in'] - Q_output = brayton1.condenser.variables['Q_out'] - print(Q_input+Wc-(Q_output+Wt)) - print((Wt - Wc)/Q_input) - # res_x = brayton1.base_params_single_optimize(fixed_var = fixed_var, - # base_params = param, - # target_var_name = 'x', - # bounds = (0.1,0.5)) - # res_min_T = brayton1.fixed_params_single_optimize(fixed_params = fixed_var, - # target_var_name = 'T_low', - # params = param, - # bounds = (30+273.15, 60+273.15)) - - x_vals, eff_vals, opt_vals = brayton1.sweep_and_optimize( - fixed_params=fixed_var, - params=param, - sweep_var='T_low', # 扫描变量 - sweep_bounds=(30+273.15, 60+273.15), - opt_var='x', # 动态优化变量 - opt_bounds=(0.1, 0.5), - num_points=30) - - # refprop_path = "C:/Program Files (x86)/REFPROP 10.0+/REFPROP" - # recuprator = Recuperator(name='Main Recuprator', eff=param['recuprator_eff']) - # calculator = CO2PropertyCalculator(refprop_path) - # cold_inlet_state = brayton1.compressor.variables['outlet_state'] - # hot_inlet_state = brayton1.turbine.variables['outlet_state'] - # recuprator.calculator(hot_inlet_state, cold_inlet_state, calculator) - # print(recuprator.variables) - - - - - # self.cycle_eff = None - # self.compressors = [] - # self.turbines = [] - # self.recuperators = [] - # self.heaters = [] - # self.condensers = [] - - # def add_compressor(self, name, eff): - # """添加压缩机""" - # compressor = Compressor(name, eff) - # self.compressors.append(compressor) - # return compressor - - # def add_turbine(self, name, eff): - # """添加透平""" - # turbine = Turbine(name, eff) - # self.turbines.append(turbine) - # return turbine - - # def add_heater(self, name): - # heater = Heater(name) - # self.heaters.append(heater) - # return heater + fixed_var = { + "T_high": T_high, + "T_low": T_low, + "p_high": p_high / 1e3, + "p_low": 8.16e3, + "ploss": ploss, + } + param = { + "x": 0.279, + "compressor_eff": 0.9, + "recompressor_eff": 0.9, + "turbine_eff": 0.93, + "recuperator_eff": 0.94, + "highT_recuperator_eff": 0.96, + } - # def add_condenser(self,name): - # condenser = Condenser(name) - # self.condensers.append(condenser) - # return condenser - - # def efficiency_calculation(self): - # cycle_Wc = 0.0 - # for compressor in self.compressors: - # cycle_Wc += compressor.Wc - # cycle_Wt = 0.0 - # for turbine in self.turbines: - # cycle_Wt += turbine.Wt - # cycle_Q_in = 0.0 - # for heater in self.heaters: - # cycle_Q_in += self.heater.Q_input - # cycle_Q_out = 0.0 - # for condenser in self.condensers: - # cycle_Q_out += self.condenser.Q_output - # cycle_eff = (cycle_Wt - cycle_Wc - cycle_Q_out) / cycle_Q_in - # return cycle_eff - \ No newline at end of file + brayton.RC(T_low, T_high, p_low / 1e3, p_high / 1e3, ploss, param) + + Wc = ( + brayton.compressor[0].variables["Wc"] + + brayton.compressor[1].variables["Wc"] + ) + Wt = brayton.turbine.variables["Wt"] + Q_input = brayton.heater.variables["Q_in"] + Q_output = brayton.condenser.variables["Q_out"] + + print(Q_input + Wc - (Q_output + Wt)) + print((Wt - Wc) / Q_input) + + x_vals, eff_vals, opt_vals = brayton.sweep_and_optimize( + fixed_params=fixed_var, + params=param, + sweep_var="T_low", + sweep_bounds=(30 + 273.15, 60 + 273.15), + opt_var="x", + opt_bounds=(0.1, 0.5), + num_points=30, + ) + + return brayton, x_vals, eff_vals, opt_vals + + +if __name__ == "__main__": + run_rc_sweep()