# -*- coding: utf-8 -*- """ Created on Mon Dec 22 15:35:07 2025 """ import ctREFPROP.ctREFPROP as ct from scipy.optimize import minimize_scalar import numpy as np import matplotlib.pyplot as plt 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'] # ============================================================================= # 采用温差效能的方式来计算换热器进出口参数 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 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 # 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