From eb2717f0994b25d00a490ad294eb7258961a13d1 Mon Sep 17 00:00:00 2001 From: ljz <425868052@qq.com> Date: Fri, 8 May 2026 15:15:12 +0800 Subject: [PATCH] first commit --- README.md | 38 ++ brayton_cycle_test.py | 980 ++++++++++++++++++++++++++++++++++++++++++ ctREFTest.py | 33 ++ refdata.xlsx | Bin 0 -> 12769 bytes 4 files changed, 1051 insertions(+) create mode 100644 README.md create mode 100644 brayton_cycle_test.py create mode 100644 ctREFTest.py create mode 100644 refdata.xlsx diff --git a/README.md b/README.md new file mode 100644 index 0000000..21a6d42 --- /dev/null +++ b/README.md @@ -0,0 +1,38 @@ +### 3 分钟了解如何进入开发 + +欢迎使用云效代码管理 Codeup，通过阅读以下内容，你可以快速熟悉 Codeup ，并立即开始今天的工作。 + +### 提交**文件** + +Codeup 支持两种方式进行代码提交：网页端提交，以及本地 Git 客户端提交。 + +* 如需体验本地命令行操作，请先安装 Git 工具，安装方法参见[安装Git](https://help.aliyun.com/document_detail/153800.html)。 + +* 如需体验 SSH 方式克隆和提交代码，请先在平台账号内配置 SSH 公钥，配置方法参见[配置 SSH 密钥](https://help.aliyun.com/document_detail/153709.html)。 + +* 如需体验 HTTP 方式克隆和提交代码，请先在平台账号内配置克隆账密，配置方法参见[配置 HTTPS 克隆账号密码](https://help.aliyun.com/document_detail/153710.html)。 + +现在，你可以在 Codeup 中提交代码文件了，跟着文档「[__提交第一行代码__](https://help.aliyun.com/document_detail/153707.html?spm=a2c4g.153710.0.0.3c213774PFSMIV#6a5dbb1063ai5)」一起操作试试看吧。 + + + + +### 进行代码检测 + +开发过程中，为了更好的维护你的代码质量，你可以开启 Codeup 内置开箱即用的「[代码检测服务](https://help.aliyun.com/document_detail/434321.html)」，开启后提交或合并请求的变更将自动触发检测，识别代码编写规范和安全漏洞问题，并及时提供结果报表和修复建议。 + + + +### 开展代码评审 + +功能开发完毕后，通常你需要发起「[代码评审并执行合并](https://help.aliyun.com/document_detail/153872.html)」，Codeup 支持多人协作的代码评审服务，你可以通过「[保护分支设置合并规则](https://help.aliyun.com/document_detail/153873.html?spm=a2c4g.203108.0.0.430765d1l9tTRR#p-4on-aep-l5q)」策略及「[__合并请求设置__](https://help.aliyun.com/document_detail/153874.html?spm=a2c4g.153871.0.0.3d38686cJpcdJI)」对合并过程进行流程化管控，同时提供在线代码评审及冲突解决能力，让评审过程更加流畅。 + + + +### 成员协作 + +是时候邀请成员一起编写卓越的代码工程了，请点击左下角「成员」邀请你的小伙伴开始协作吧！ + +### 更多 + +Git 使用教学、高级功能指引等更多说明，参见[Codeup帮助文档](https://help.aliyun.com/document_detail/153402.html)。 \ No newline at end of file diff --git a/brayton_cycle_test.py b/brayton_cycle_test.py new file mode 100644 index 0000000..9de7a8c --- /dev/null +++ b/brayton_cycle_test.py @@ -0,0 +1,980 @@ +# -*- 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 + \ No newline at end of file diff --git a/ctREFTest.py b/ctREFTest.py new file mode 100644 index 0000000..958986e --- /dev/null +++ b/ctREFTest.py @@ -0,0 +1,33 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Dec 16 15:46:51 2025 + +@author: LJZ +""" + +import ctREFPROP.ctREFPROP as ct + +path = "C:/Program Files (x86)/REFPROP 10.0+/REFPROP" +R = ct.REFPROPFunctionLibrary(path) +R.SETPATHdll(path) +fluid = "CO2.FLD" +hfm = "HMX.BNC" +hrf = "DEF" +z = [1.0] +ierr, herr = R.SETUPdll(1, fluid, hfm, hrf) +R.SETUPdll(2, 'SI', 'SI', 'DEF') +R.SETREFdll("DEF",1,[1.0],0,0,0,0) +mw = R.WMOLdll(z) + +p = 20000 # kpa +T = 391.05100419252705 # K + + +result = R.TPFLSHdll(T,p,z) +s_SI = result.s / mw +s = result.s +result2 = R.PSFLSHdll(p, s, z).T +rho = result.D * mw +h = result.h / mw + +x = min(1, 50) \ No newline at end of file diff --git a/refdata.xlsx b/refdata.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..afe0d28212fa33aba7935fdd615c763cd57315b9 GIT binary patch literal 12769 zcmeHtgLftA_V$TwJL%ZA(dn>byJOq7la8HqjE-&FM#pBy=9ig!?>93&^ZN_#tvah} z)jGAGRkdsH{k$8pk|3a{0B`^#001BYB z&|Ztq#nOT>2Nam%8vywI{{L;}5Z{>?Fs8whfpg^rRm%CA3FU{F%Q~eF&C51n!I(o> z*uAEhxp^zV~kMr8IA`%GidCD^mfp@xT)q|qoQ zArPUdeN@vvlzpVF<+aH=mkEw3%wJ~|i_<}F|Ehm3gBx>-cYfX4jw-$EgM+Z0nz}p)LK=v=5TBpQ7eEn{IQtyrk^X{oSc7_)A z^mM-+|L39q#cuhhS1*l~mg!-D4LTEh4j#OlTm6bCB?GN(FyUpu)czKm4@@SCgc8jGn92u3HxZb%eDCyDG5rT@$E>YCBbh8)9Y3^q3 zHszCqD}_^AG}ZUUqHM{bbz;$(3y~_s2|5*QD3pBMAavdoU-f=zwY4wzn3O~z& zDjV3p9mP$gc`qau9K!L3a7vv_r=Shme=%Mv_ZYMwx_Q7>QZV5(t^AT{$MwNg$H1cV zN+hiV<;8>cbNZnC2NvXerio7jq~C75)oWQUMlxJ`*`a#A51tSCg%cF5y_448oh08y zkPQn206>HQ01)0a<6=SYY;9-$#oF5ZH!3St(y~rtK=sV3ebar&@QWk@OOPTWP|vq4 zskd@3R^;rbR0d+!PZgPaeJU3hZ`jR$3PHy4$k?{`n5d&)J2`^p6`P@4HAqee!yFPl z5Mp)4di7|#9A`YDv4++@6g6_8=YBSNe4Oi3PJ61Gzk>v%()HxmSsp!QP#BK2Fs}}S zf55Rcv)Y_ndC)Q^p>LL#&`{xn^INQFGCKz%Vx{kzOp8Re%l@b{(vL3?QQ8Rhs8#&K zO5}T=C-$)vlY;nGAhG65DL*cG5RxX%ll^FC2r*G3}0(leq;gRfdoQ_YQ?gE zsf@KnW}jPxJeljbPV${3+Qp31u92_F);s=ooIKS`7dq8p(W)2Bm{6}+kv?_c6?v7e z-l6DtO>#VKi$xM%oT?ZX^)H>~NyX5)Va`r!6x89cG8;qm}i233|C?_(} zPVG!ZUvbt{Zm}bHx`hpUZH0*0Q6d<>m$JY$*p=hJ@Oez%RAk~*?&pIC=SE*@R73=H z1AF@_n>2|%G`DQki?kIpwMtm8`OBJG+4bBhvo~C|743}4^zRfZokPs;bRXUbp8L*j zJ~&@=j^*XPl$Q#gdo;N~jPhMTh(YGHxs)tOF4`imv32L|+{TF1%!O4>rUl2lTC_vA zs(&Ib8Zr09)&~?ybeZT`39TMfKaj?cc#GjyxugY| z0|mPel$fs6?dyUE9PwvRM=9p#{I1lX4> z3iy3WSv_|6v}W_qE~sz5m@?Yn+T^F9TxVncq6m=xRY(bBIMV#K*<&ua)Aj(WvOTefRo@4R@JmRF4z!kl4eT?z>$A{4Fr2 zt7=6bu_1a^%)TM+9@ZmAb?D~?9s!5C5l|^9@n_Z$(M*8}#hsjLDZQNGkq7jScgs(n zA}2BrWNxQiEF$oGE}yj~m6l#^F%h4i@J zqyvVcpJ|c&q01VpV3!uR&nT{NsYER6ZYkSJeApssmf_7Tjo^EfW!9QdgUoQn#DKW^ zNwmO4^4aBM8PQFfC)dlR1I6%;ErgdhF^_ngHR-eLz2Pdu*xM)KpI5asC0`VXb>V+- zMz1dT!pB~ZvfRFrVR!p8xetv2-zr|WQzF^!urM#s`>4&}OzD|NpX1U&~!)5ZK zL`Q$TnDe**>+cefGi;$@joBnUzm27>S+&~xge>%U7KDoL%fP`5pT&LK7`mY&X`Nvi zyRY@4aq(!M$Zu&bmV^ieRGzY!T)wWXmLMeCy%+*5oQ}?|>)Wx$$G!w6nBbB2uyUt; z<=B?-z#$ry_Z1{^pkiRB_T}QQjbPEWmb(Gd4@^m7jA-(dDgDTa z!C~RhsgH-K-!zyOl9!+*$T_$JQV`-Q3!-pFhy?seCTFkdRRoiyTcCdc!UE7g82C}4 zaRt&_WiMRkiXAa1BCaMwj#9Q!WGx|_eVs^Nx(tl+z{`n=?o8~2FdVb>;sg~gy z9ig3kRJBZsYYgRm0{O6I^GkjIzSw|4iw2&*h>@!I693bm`meP5@UQr zm4mHL&=+H6pHx3*b>wj1jy3VipGRiO6f8a#UY8%mB|8ma(Kn*rYpgv$lrHw)18LDX zLDc2KH&$@`nGkmn0|$xJFU3HSj5{OCZA;Lwahz%!sWxa>glb_5I}*NBd^R(hyp)R6 z&~`+aVNyObLbQQ5W#Jx`euUNZECoYZgpLV3Ugl{F(_XM7BK$#YGTn|0Yfip1%7$nX zGF|^Hw6=ddW^5>=4{PR-?P#B-m%BR>GoS^90}@J$Z^vjbTP@pVwkUmap;tddUYyjT zANVtF>fp$-;LQ9=OnziWZbqDFICH-VL<`FOuG6(vhZW+t3S^O-++NiozsGTwAW)@4>-J{`IW6R@iO9htyF@I-xU?UJNel zRhD~`ouUOoEUVNC&|Tsmp*YlpDmj*AW?^SeD|83Rk)#T%W`Mh4p>cu|RU|`0ThbO2 z1KvVnf>5m@5=8xyjXsH!Rdu)0U)INcZn^V~vn&KP&b(rb;|yf@u>JY#EtV<5DT^cl z?3y!eS$9SSje$*k47NfSSgsjvcOt16tm~KWB4|Vr`gY>^?cAF);Yp~4akpU+X~=t} zTG#MKlj@69H#je}{AMW9JS{sKFwfgco!f>3-I_P^{;3esrI1?O~!$Rc-70l9t7 zC&Z>te|*Z=3M%e~Q*0?qt+@&t)cx2@8Yp%e^)pO(UphH$FwqvY|Et()F^(y-jEmqA z<2a);b8B0k*_l?9Mq_x{6)IgmTVFQYf|e}naRc!RAzqlFEr-@Y>PW8sW&BHQyv#%P zrPESDQ-)t3-7*eM#6}Zc)ZX~8rZ6u$t0tdgR=tf@T9*=>sO&7;M>%JL(fp+cV3${d zcUt;;B4O1c=hN{{TxzZp_?Sb8 zIjQIv8Z2%*tXV^ZHZ_^zwgvgyOHTK`ms36N9bb7ZxWUG;gdu2!Bbn883nNzK_INHl zUPa3-*QlHB5zN2R@daodlS-9Uv=MkWx`xG_9kRjJpu2G7X_==P4n$cNcCoEIQHDlu z0-aiH8JroPGdDHN=p#QKub6byaOOv~-R1f#xg3xXTp(Gq$sxjTZ|D)=S9gmepFYFx zWm)CwtGg{j5ZE@{#hoDxI3b=zZ*0;1__!o_I&1dRfBgCtk^FTPe*R@NkDxXa<$hq3 zcv#P_m2sL6I+i1n9Sdl}iE)4LvaZaMpTpg({7u*@xi;>ckqKB^^I}3ew77=e%a)Aq zaFmmqIhy*AN94Ht=1jm{$jxK6VBN}5!%lShdPB)jSe~h4P$4 zxAekeKNq#-Qd=62*K1gJVeMsG!h>^1a0|nr5T?HxElQ@d?ei*PAwcdtc@K1Z3ZeRd z5jUVkYax~4TjAS|fASeafG2+LV;in@eDHOJxH+EHbFkK}#xo04?vo`r2uu!I=^Ut} z@e32|3?K(plw1>QK(G~r{@mM~Amt~x+#Wy1Rpk6qy}7Mmz=k5jMuT54OZc5PhVx}; z2CZPm3Uawj;Ztg3Idz`RdOawmon&Js?K-YRChUeA6XI8rwPc)#tZ2Z39q66Dz5VAz z8(uX#7rRk+WAW1=*z%;{i*_u7dXG$-w$Yg@+y=~ZI;k2t1Ch3R^_H(p2=Uoj4{d2A zUS?hBY+Vl5(n1D!;1i<(fm^+kn#-r)$jOpi3ur zcBBg%1CKwQorf6%ssI`;+ff0PaD5q!zf&n;6i~Iy_GKk8;!-ZUgiIYtg}NVi2-P{^ z+n}2Vk+qwSRYzk-giqB_A$-foHoSe{7Cc7`g*;hcv5)S`wHNl4Q`Vcsous)a?u34n zkz&`z52VOYnTuXy%w&CCUtdyiWwdW(93-_t8$w7S(Z`34E6)A!hh;Hqqsf=z|w|Ht8elHatifLqm);gd8F*^n@G|>*)gS18BXx zQd-J#nN26S zVH3H>PRY|s+CVsKn9g)0zM#L-B41zN0{15ykX67@@1H2PgpK)87s9Q;@f z^g&=E>~^EFpRhy?SfX3)yd*Gys~}s-U9)+(C<+uuXRtziJa1aWzAr^K6}sHEj9#-f zpL16JYpU>K7X%w;04iN3u>Clo|MGg_!yp^PFxTeSQ^I45!Q@||MQmuW#`>51bu%q8^wkv(mu$~Pt#a>@Ri-tj%ZrsCWkBu^wni9GplQ#d$zY8)aY^ z=mZ9lgzd!Ayuck}RC{gL!c9aoCLEx}MK+cL1nZ~We@*dccX=Duj({-GKm3SAimHOR3|$nW<&U+iSE30+f+e=~9kT&m zm&S~a3HZcE8^y^5-o{9)ERnCyC`xkC2HN-S_J@Ix=;G7M6uq~{&3OX?Ul6~yJLjjY zZjn$_GhLB~_EaDZKWuJ38tMSm5W3wpq#Rm;0-Ty{qIP1>=nq7j;fXra{9d+CUyv*f zg%qv6vhILClE;mJ_6MW!l4s$UlEL{g6xFFbvib}=YIdwu+?Z2)V*7{^tD2gHm(^$= zf@94C@(J*o-EXw(EH>d@=2+8~vmfiuvEz>r=Fj~FFBv*DMUn?i5z z3*Wna1fAELjT(oe85g~+ERUzDxGZnaR5w-I)3AXo@3)*?unWpb1m-zPf-7y4A*)Iw1OyplEHqa7+~$ljfp zvOp<8lDY*dDdRZ+&*)To>M2CiaT?fUf$9KUM!lA`P0-^k2SBRS=sF6k!T&Re4 z{bco$E23F>wdWHXay?iM{C04%x>aky9&#UDf-tJW%w~(2C>`+~hKIq{qmWQuVKP$% zf4@+pLAkAJ1=2NjC&Ovw+0;i1;wF+{x{Q_(MwBJ};8n;z=E5%A4E5y|AGCtE>s$QIRtG5v6@<#djk zT;?qG;N)u4%i!+X6Tpv%7_3}^n zGOOWkPiQ|QO#qWJZl$s41Qux&!xW~+NmNWtWka<}e1i_~c+Wyo0eGTzckHHZ_(aDi z@)Pd~0gP4Ond>d)$;dK&Grb^qMcUCy`?Et~VO z^oev@Fc38H_q9;*%h9cU#U+$iZo}vnA5MS2fwFrgi45uZYWtGTrGO`dL{?rFnR0y& z@p;bU))LMNv*$NYr5aQEtOch1x=~l!ZMP>k+Bzxib^*gUXIF;_^0mdw6hw);UB4}@ z>=r)YHIsf2%-Nq5t8toj0AnPNQEVx6tKpc$3y!!mB&?sGOu);stF3fDe=KmLYQ;jC zWSs74&JB|}ocAnrZ5`-0Ev(U0)*71)1RnBeMJ?~omLyehD9dy{Q$BL!cPMEol>~`D z3CGlvGQyE_#FW#47edW63ea9kQJr#7osvCD%2pS(-GyebL_>4gz1{fOwi_(?Lomn) zH-29qt*|zwHKH7PJVn$kj=4AmNcIaoGhwuY!cw!XtC^!|3Gb&oTH*+Ynn;Fe20x5V z8t$>gC$;|Ol;Bw=nWisIICjiOkqb!XA3qc}`jU|;O!jhHJ=B>&CQ%c!m_MK#VD=2S zfHx7(YGC9jGkRJz+I9a?7>MG?1_|vuyX4j-njA7uDqagBj!Mjsi5t z5*hc8dn>JTMi33YJZ?7L9EI`+vfUZ1pJcD|*}qa+cA_^C*I*;@x!QkLe_fSdgl|r$ z)SzJP_97|*&le)LP=Wc$P@4Fp*ED;qc-qzVA-55?py3vzvHM&2$*Fd9#p?RNh(;<5 zF*-&{dltP}&&`Wx(^hi_K}nkHIQSxC1G(d?x%N9Z>7hBch) zX1!d_#(nos*toUdSHI0$qnhnnH6YGP3_pw3$RbcBILTVGKMFo+EBAK}D?AQ*0G-Y` zS^hH}S07agh=KwDqHqAf`|r1Amc4_kg`xd#iRY+_wDmF@W;?>7HzDZB!9h==F-~k4 zSEvs_OOA%FlA2bIO)iXuf-8pmgbntmXVQybMvzuGUcJtYA-_Jo_)u(|)ddf%IUTxA z?8H;+y3}9J(XUXaFRe?}Ynm@QE;PD8j!q>9Uh z3S|J5Q-}E(?GUXNb35A369HYCig_Klb@~3R<3Tdrk+2hXyiX_DT4thE3S_FxZdoI@Zi`l4&Vz$ipdn)1Hm(jVcLhr@UA+I)Glaq*C znR*&vhzDS4UBud$5QJ||~&;r>Qd^k{hl-X|+B=f_X%>o2cY zcb1}lDI>IoVRjjBMp7#sSFfur{9e!XN!@KQnOGuEIGw~5El!3`30_vIqV>k$q;YY; z)+y>I?D~aIuYNlY8bTqd3iJn&vfv+_n$L{gzTE4N0hpfEz-aNkrj(}_A0F*qAFj8U z;J{I=`LH5_NN}Zk9{LDzGbZTiri(FL_SdiR+upW@%-0Sk9iFF3znn89MpnT+9DhdV zZVSw}euk)hS}E*WbjN;NBgAui)n)0A>oVYhg)lvUz}TL8W|Hc`zdP?OtlEEU_cSoS z+?48yi;!FEFNu`Zffn@JUB}$4kFja-_KD^H$8JrDEu0_F`*;Zx6aaww?|@^XXJ=@j z;9zHJW&C^G=}la=nh`(@I*V}$p3?3L+{Tsomhe-%QACkNKE7bE832Z44txn8>B$1s=>(QfucC!^OFBd2_CJlu|P=0=g@_^Ipz$7Wzza--sFDpZV` zm^(2a;^HJGW^7@(Y?OhiRc>E<(9$_NRP&sz!8@($e$D16N?Cg{R<@*G$(BnkJvq8n z?3l~ZGv^jz{v-()F{xroRDWL|`R_@j{R_SnHdN3QT+E>)r||JXG3t)Gi1(;qUKy*( zQj8x*w20HiVr}RX=aI}H#0YD~J(Kx8+rNiOjv)KlQza?Ybc*_BoEvX4@4#seM6E|?bV#dBX^h}%x=(cV4}z8B?=dg ze+2{N7<^NE(xe#>nM3L3P0g7mZ4J)zp!Fq;DVw%xX}+`;CjVm3X^UzI5x=ll0<}YI z?zA;n`^^NUWj7G6bu;tPzwKK1Y*YKOoqbVq{n#zo-STUP81vPUlgh>wzs&V}M)>zS z>H3NiDb%|$?Y^50;=fErU(Z5c*hJ6N>NlFQh+VY`AVd^BOS{G!TyaKPw(|aP+E;Q7 z^aFxWOa*KKUieTKDP-=;b4pqj>lOc2Pjp+K#i!g3l)aP-~YDOSkU zbgf--#NK%s9`5jRn$VCS9`_OsGIGTBq$_`i{E=eer$eWiQt8Hr^b-nCjvpwFHFMRD z8S`GaCYgrrTtNT|n2$8U4eo@z!MKBnsQ|6>9(wIAVsLlZOa@2|&76e@-3O4~4?1L^ z^U$Ses??Rrx@Da|$jmC1+2qoEilhZ*$q3C*5bK~^nk>=)=xdxOS;#@+CvJGf_Rn%- zY+9J7I&Z=@9cFQwRb!+5N|teP5Z{gw)L-9~@^`h_e}m@{eOH^>yV{W6JM0G5`m%P` zHum)T)^>)!>Fs@_>wk^)chv>N>PYr5p!ThUbqEi+rxlbb2rc^#;~Ri>Kjpl0o@gV{!ap8<0=66$H3|VGx)it+jkn9nUih-RPBfuP2KF5r z3X98EaOC-sYB%gsw5P@j)kMSZ0#_V+SSbj8o+_+r6QU>uDqWb%cHtmZ?QYh+6!T5j zy;Ba8&IT?K$_VS!r#&^qL*GewI*EOzW1aj>i^V z;cWciH%Y1IgL*BqyBAB>d(ffXJhj}+k02p zzqj<*F0&Pecd83)_n+dWKC2I?9E3KUi4_*caEfc8ieAX6Cdmt38mqaU7ACU(s<2kz z%Og6#X50ERdvba7xu-&gnRxsza!fsB$VH{{fyNk%JMCj!D8-}sn7NwDuoJ+W7L<|c}%F(+ch9u}zb*v=~|2~#$kH;M!QztOqFBMHi|Lfq zh27I{Wa@R|*|4z!#tS`&P3hq{6t2B6K`#bssMkSv;YlKF6 zD+Rk{@4Q8==HiT=8v{T1QY70n+A z0u;X^{AGRftMp%q{SR3Hz?l*N_#4gtD*o5d{AclN>OYD9J5b9?g1uW9008^GgT5Q- J0Nrm_{|~&Rswe;e literal 0 HcmV?d00001