用MATLAB实现火电厂热平衡计算：模型、迭代与工程实践

简介：面向火力发电厂热平衡计算场景，这套基于MATLAB的64位程序为电力工程师、热动专业学生提供了一套可运行的建模与分析工具，能将锅炉、汽轮机、回热系统等能量转换过程程序化，替代繁琐的手算与查表。资源共23个文件，压缩包仅2.74MB，以12个mexw64编译模块和4个m源文件为核心，辅以3个dll动态库及IAPWS-IF97水蒸气性质手册，分别承担数值计算、模型脚本与工质物性查询；另含txt使用说明、jpg热力系统图和ico图标，便于快速上手。已有2045人学习下载。以600MW八级回热抽汽机组为示例，程序支持设定燃料与蒸汽参数后直接输出发电煤耗、热耗率等关键指标，可用于不同工况模拟与设计优化，并预留了二次开发空间。

1. 热平衡计算为什么值得自己写程序

火力发电厂热平衡计算，说白了就是把全厂蒸汽、水、烟气、电、热量这些进出项一条一条对平。听起来像是教科书里的基础内容，但真正做过一次全厂热平衡的人都知道，手算能算到你怀疑人生。尤其到了 1000MW 级超超临界机组，回热级数七八级、再热两段、还有给水泵汽轮机、轴封系统、辅助蒸汽联箱，模型复杂程度完全不是一个量级。更关键的是，工程上一旦调整某个参数——比如给水温度变化 5℃、凝汽器背压从 4.9kPa 改到 5.5kPa——整个热力系统的阀门开度、抽汽量、热耗率都会跟着变。手算一次要半天，程序几分钟就出结果，这才是写热平衡程序最大的价值：快速、可复用、能对比多工况。

用 MATLAB 来干这件事，不是因为 MATLAB 是万能的，而是它在工程热力计算里有几个天然优势。首先是矩阵和数组操作内置得非常顺，热平衡方程本身天然就是一组线性或弱非线性方程，用矩阵表达最直观；其次是 MATLAB 的迭代求解器和数据可视化配套齐全，算完后直接 plot 热耗率随负荷变化曲线，不需要额外接绘图库；再有就是大部分热动专业的学生、电力设计院和电科院的人在校和工作时都接触过 MATLAB，团队协作起来沟通成本低。

这篇文章的目标读者，我认为有三类：一类是正在做毕业设计或课程设计的热动、能环专业学生，需要交一个完整的热平衡计算程序；一类是电厂或设计院的技术人员，想用 MATLAB 做一个内部校核工具，替代以前塞满 Excel 公式的工作簿；还有一类就是单纯想搞清楚热平衡的迭代逻辑和控制策略，希望有一个干净开源的基础框架能二次开发。无论你属于哪类，看完这篇文章，你应该能搭出一个至少可以跑通朗肯循环加一级回热的骨架程序，并且明白怎么一步步扩展成完整的多级回热模型。

2. 程序框架与热力模型拆解

2.1 系统级模型：先把边界划清楚

热平衡计算最忌讳一开始就陷入细节。写程序之前，必须先在纸上把系统边界画明白。一个典型的凝汽式机组热力系统，主边界就是锅炉、汽轮机、凝汽器、回热加热器、除氧器、给水泵、给水泵汽轮机，再加一个发电机。至于辅助蒸汽、厂用蒸汽是否纳入，取决于你算的是全厂热平衡还是汽轮机侧热平衡。

我在实际项目中一般建议建立三个层级的模型：

• 汽轮机本体模型：计算各级组进出口焓、熵、温度、压力，确定抽汽点参数；
• 回热系统模型：计算各级加热器进出水参数、抽汽量，这是整个热平衡里迭代量最集中的地方；
• 全厂系统模型：把锅炉效率、管道效率、厂用电率、机械效率带进去，最终得到发电煤耗、供电煤耗、热耗率这些考核指标。

这三个层级之间是单向依赖的关系：先算汽轮机各级组的膨胀过程和外置参数，然后算回热系统的抽汽量，最后汇总成全厂指标。如果发现全厂指标不满足预期（比如热耗率比厂家 THA 值高不少），再回溯去校核汽轮机各级组的效率取值，形成一个完整的计算循环。

2.2 核心数学模型：热平衡方程的矩阵表达

回热系统的数学本质，就是每一个加热器都满足能量守恒方程：加热器蒸汽侧放热量加上疏水放热量，等于给水侧吸热量。以典型表面式加热器为例，稳态能量平衡可以写成：

Q_给水吸热 = D_s × (h_s - h_sd) + D_d × (h_sd - h_drain)

其中 D_s 是抽汽流量，h_s 是抽汽比焓，h_sd 是疏水比焓，D_d 是上级疏水流量，h_drain 是排出的疏水比焓。

把系统中所有加热器、除氧器、凝汽器的方程联立起来，整理成 D 为未知数的矩阵方程，形式非常统一。这也是为什么我说要用矩阵来写程序——一旦模型从三个加热器扩展到八个加热器，你只需要修改矩阵的维度和系数，不用重写求解逻辑。

在 MATLAB 里我最常用的做法是：把每个加热器定义成一个 struct，里面存入口、出口的给水焓、抽汽焓、疏水焓、端差、流量等字段。后续所有核函数只用 struct 的字段做计算。这样模型可读性远高于几大坨零散的数组，而且后期添加新加热器只需往数组里 push 一个新的 struct。

2.3 汽水参数计算：蒸汽表的合理调用

热平衡另一个绕不开的基础问题，就是怎么获得水和水蒸气的热力参数。工程上教科书里会给出焓熵表，但实际编程完全靠插值表既不准确，也要写大量冗余代码。比较稳妥的技术方案有两种：

一是调用现成的 IAPWS-IF97 标准计算库。MATLAB 环境下常用的是 XSteam 工具包，m 文件直接加入路径就能用，函数接口简单，比如 XSteam('h_pT', p, T) 可以返回对应压力和温度下的比焓。这个包在山寨电脑上跑起来也没问题，只依赖 MATLAB 基础环境，不需要额外编译。更严谨的做法是直接调用 CoolProp，它支持 Python、C++、MATLAB 等多个接口，IAPWS-IF97 只是它的一个流体库选项。但 CoolProp 在 MATLAB 上的安装对 64 位环境相对友好但也偶尔出现库加载不上的问题，如果你不想在环境配置上耗费时间，XSteam 是更轻的选择。

二是自己实现一个简化版的 IF97 函数。IF97 标准把水和水蒸气分成五个区，每个区有自己的方程形式。这个工作量并不小，但如果你的项目需要脱离工具箱部署，这种方式是最可控的。我的建议是：校园项目或企业内部工具直接用 XSteam，不需要重复造轮子。

这里要特别提醒一个我踩过的坑：建议不要直接用自由度很高的通用拟合多项式来计算焓值，比如用 30 阶多项式拟合某压力下的焓温关系，表面看误差很小，但在临界区附近和湿蒸汽区误差会被放大得非常明显，而且不同压力区间需要切换不同多项式系数，边界处容易出现导数不连续。热平衡程序本身是收敛迭代的，参数微小的跳变都可能导致迭代震荡，轻则增加收敛时间，重则直接算不下去。用标准库，省心得多。

3. MATLAB 核心函数与 64 位环境适配

3.1 求解器的选型与对比

热平衡方程组中，如果假设汽轮机各级效率和加热器端差已知，方程近似线性，直接用矩阵左除就行。但实际中往往有两个非线性来源：一是汽轮机各级组效率随级组流量和水蒸气状态变化，二是湿蒸汽区的排汽焓需要迭代确定。因此程序不能只做一次矩阵求解，而要做内外嵌套迭代。

我建议采用如下两层迭代：

外层迭代变量是各级抽汽流量 D_i（也就是各级加热器蒸汽侧流量），内层是汽轮机级组出口焓的计算。MATLAB 里 fsolve 可以解非线性方程组，但它对初值敏感，你给一个离谱的初值它可能闪退或者返回错误解。我在工程中更倾向自己写迭代循环，而不是直接套 fsolve，原因是热平衡求解的物理意义非常明确，自己写循环能掌控收敛过程，排查问题方便。

自己写迭代有一个技巧：把 D_i 看成一组变量，用上一次迭代的 D_i 求解出新的级组焓、加热器出口水温，再求解新的 D_i，如此反复。因为热力系统本身的物理特性决定了这种 Picard 迭代在回热系统里收敛性很好，通常几十次就能到达 1e-6 级别的残差。

3.2 64 位环境下的内存、精度与工具箱兼容性

很多人忽略了 64 位和 32 位 MATLAB 在工程计算上的差异，实际上在热平衡程序里主要体现在这几个方面：

一是 double 精度。MATLAB 默认 double 是 64 位存储，32 位版本同样支持 double，但 64 位版本在处理大型稀疏矩阵时的内存上限高得多。热平衡程序本身规模不大，单个工况的矩阵撑死几百阶，内存不会成为瓶颈。但如果你用 fsolve 处理带约束的优化问题，或者跑多工况批量计算时把每个工况的结果都存成了 cell 数组再加自变量的网格，内存增长会非常快。64 位环境下建议打开大型数组的预分配，把一维数组预分配到最大工况数，避免循环中动态增长。

二是外部接口兼容性。热平衡程序经常需要读 Excel 里的设计参数、汽轮机厂家热平衡图上的抽汽参数。32 位 MATLAB 在调用 Excel COM 接口时没问题，但 64 位 MATLAB 对老版本 Office（比如 Office 2007 之前）支持的不好，经常报“服务器抛出异常”或找不到 ActiveX 组件。如果你还在用 Win7 64 位系统配老 Office，第一件事就是确认 MATLAB 版本和 Office 版本是否匹配，否则直接在 MATLAB 里用 readtable 读取 CSV 格式，可以避开 COM 接口问题。

三是.mex文件的编译。如果某些热力计算需要调用自编的 C 语言蒸汽表计算代码，64 位 MATLAB 要求 C 编译器生成的 mex 文件也必须是 64 位。老机器上如果既有 32 位又有 64 位的 MATLAB，注意千万别把 32 位编译的 mex 文件直接拖到 64 位环境里，必然报错。解决办法是统一用 MinGW-w64 编译器链重新编译。

3.3 主程序架构示例

一个足够清晰的主程序，骨架大致是：

% 主程序：热平衡计算 % 设计工况数据读取 para = readtable('design_case.csv'); p_boiler = para.p_boiler; T_main = para.T_main; D_main = para.D_main; T_reheat = para.T_reheat; % 汽轮机膨胀过程计算（各级组） stage = calc_turbine_stage(p_boiler, T_main, D_main); % 回热系统迭代求解 D_extract = ones(1, 8) * 30; % 初值：每级抽汽30t/h for iter = 1:200 h_feed = calc_heater_enthalpy(stage, D_extract); D_new = solve_extract_flow(stage, h_feed); err = max(abs(D_new - D_extract) ./ D_extract); D_extract = 0.5 * D_extract + 0.5 * D_new; % 阻尼迭代 if err < 1e-5 break; end end % 全厂指标计算 result = calc_plant_performance(stage, D_extract, para); disp(result);

这段代码里有几个值得展开的设计选择：

• 阻尼系数取 0.5，即新值只吸收一半。这是为了抑制迭代初期可能出现的振荡。如果阻尼系数取 1，即直接用新值，往往前几次迭代会出现流量负数，物理上不可行导致后续计算崩溃。
• 残差用的是相对误差abs(D_new - D_extract) ./ D_extract，而不是绝对误差。因为各级抽汽流量数量级差异很大，比如高压加热器几十一百多吨每小时，但轴封加热器可能就两三吨每小时，用绝对误差判断收敛会导致大流量级早已满足条件，小流量级迟迟不满足。
• 主循环只迭代了 200 次，结合阻尼系数和相对误差判断，在绝大多数情况下 30 轮以内就能收敛。如果 200 次还不收敛，说明模型有问题，直接输出中间变量检查，而不是盲目加大迭代次数。

3.4 多工况批量计算的实现

热平衡程序另一个刚需就是变工况计算。设计院最常用的一种计算，是 100% THA、75% THA、50% THA、30% THA 甚至 20% THA 工况下全厂热耗率的计算，用来给汽轮机性能曲线做校核。如果每次改一次主蒸汽流量就手动跑一遍程序，人能累死。

高效做法是写一个外层循环，把主蒸汽流量、主蒸汽压力、再热蒸汽温度、背压都定义为数组：

loads = [100, 75, 50, 30]; % 百分比 results = zeros(length(loads), 4); % 热耗率、发电煤耗、供电煤耗、排汽干度 for k = 1:length(loads) para.D_main = D_THA * loads(k) / 100; para.p_main = p_THA * (loads(k)/100) * 0.9 + p_THA * 0.1; % 注意：滑压运行时的主蒸汽压力不是线性下降，应按实际运行方式设定 r = run_heat_balance(para); results(k, :) = [r.heat_rate, r.coal_rate_g, r.coal_rate_s, r.exhaust_x]; end

这里刻意注释了滑压运行的问题，因为这恰恰是初学者最容易犯的错误。实际机组在 100% 负荷和 30% 负荷的主蒸汽压力差异非常大，是定压还是滑压直接决定了汽轮机入口的节流损失。你如果全工况都拿额定压力来算，低负荷热耗率偏差能到 2%～3%。所以变工况程序里的压力-负荷对应关系，最好的来源是汽轮机厂家热平衡图上的标注值，其次是运行规程上的实际滑压曲线。

4. 工程数据准备与误差控制

4.1 输入参数表设计与单位统一

工程计算里绝大多数 bug 都出在单位上，这话一点都不夸张。热平衡程序的数据源来自不同渠道，设计参数来自汽轮机厂家热平衡图，锅炉效率来自锅炉厂性能计算书，辅助系统耗汽量来自管道图。厂家给的数据单位各不相同：有的压力用 MPa，有的用 bar；温度用摄氏度没问题，但焓值可能用 kJ/kg，也可能用 kcal/kg。我在程序一开始就写了一个统一单位模块：

% 统一单位：压力->MPa 温度->℃ 焓->kJ/kg 流量->t/h p = p_bar * 0.1; % bar -> MPa h = h_kcal * 4.1868; % kcal/kg -> kJ/kg D = D_tph; % t/h 保持不变

然后全程代码里只允许出现这四种单位。注释里也要明确标注，因为自己写的程序几个月后自己看也未必记得清楚。另外，厂家文件里的压力值很多是绝对压力，但有些图表上标注的是表压，这个差异也要在读取时注意。工程上热平衡计算一律使用绝对压力，如果原始材料没有明确说明，默认是绝对压力，但最好找厂家确认一次。

4.2 收敛判据与残差监控

判断收敛是不是只看抽汽流量就够了？不够。我一般会同时监控三个物理量：抽汽流量变化量、加热器端差计算值和设定值的偏差、整机功率计算值和目标功率的偏差。三个量都收敛才认为这个工况算完了。

给一个我之前项目中实际使用的判据设定：

有人会觉得 0.01% 的目标太苛刻了，但工程上这台程序是拿来出热耗率保证值校核的，精度不够厂里不认。0.01% 对 MATLAB 的 double 精度来说毫无压力，唯一要保证的是蒸汽表算法本身的精度足够稳定。XSteam 的 IF97 计算在正常参数范围内误差在千分之一以下，完全够用。

迭代过程中我会用一个全局结构体记录每轮的中间值，方便最后画收敛曲线。有一次程序在 30% 负荷工况下怎么都不收敛，最后就是把中间值画出来后才发现是第 7 级加热器抽汽压力低于除氧器工作压力，导致抽汽量反复在正负值之间震荡——这属于模型边界设定问题，不是求解器问题。

4.3 热平衡验算：程序算完怎么确认结果靠谱

程序跑完拿到热耗率、煤耗这些指标后，最关键的验证工作就是拿结果和厂家 THA 工况热平衡图对比。如果一款 1000MW 超超临界机组厂家给出的设计热耗率是 7350 kJ/kWh，你程序算出来 7350±20 左右基本可以认为是准的。偏差超过 1%，先回去查汽轮机各级组的等熵效率取值和管道压损设得对不对，再查加热器端差有没有设错。

验算的另一种方式是做整体能量平衡校核：把锅炉输入热量减去发电机输出功率、排烟损失、散热损失、厂用电等所有损失，看剩下的热量是否平衡。如果有超过 0.5% 的差额，说明程序里有个别能量项没算进去或者算错方向了。

我自己的项目里，程序最后会输出一张热平衡汇总表，格式类似：

把这张表严格控制到能量守恒，后续审计和报告输出都会方便得多。

5. 常见问题与排查实录

5.1 常见报错排查速查表

热平衡程序运行中，最典型的故障集中在 steam table 调用、矩阵维度、收敛失败这三类。以下是我实际项目中遇到过的典型问题：

排查这类问题时，别一上来就翻喷代码。我建议先把参数范围打印出来，对比设计值和读入值，多数问题出在数据读取阶段而不是数学求解阶段。

5.2 几个必须说明的边界情况处理

热平衡程序对边界情况处理得好不好，直接决定了程序的实用性。第一是湿蒸汽区排汽焓的迭代。低压缸排汽通常处于湿蒸汽区，干度在 0.88 到 0.93 之间，直接查 XSteam 的压力-温度点会出问题，因为饱和压力下温度是固定的，而你给定的排汽压力往往不是饱和压力。正确做法是由排汽压力和等熵过程计算排汽焓，然后用干度公式 x = (h - h')/(h'' - h') 计算排汽干度，并检查干度是否在合理范围内。如果干度小于 0.85，多半是低压缸效率设置偏低，或者再热温度设置不对。

第二是凝汽器背压过低的情况。冬季循环水温低，凝汽器背压可能低至 2.5kPa，此时饱和温度只有二十多摄氏度，蒸汽比容大、容积流量大。某些蒸汽表库在低压区拟合精度会下降，如果计算结果跳动，建议对该工况单独调用高精度的 IF97 分区判断函数，不要统一走默认插值路径。

第三是机组跳闸或深度调峰降到极低负荷的情况。机组可能无法维持额定主蒸汽压力，热平衡模型中的很多假设（比如加热器端差恒定）不再成立。我的程序里直接加了判断，当负荷低于 40% 时弹出提示，提示用户确认运行方式，而不是闷头算出一个数字出来误导判断。

5.3 程序性能优化与部署心得

热平衡程序虽然计算频率远不如实时仿真系统高，但批量计算几十个工况时，性能差距也能到几分钟和几十秒的差别。性能优化最有效的手段有两个：一是把蒸汽表调用次数降下来，同一压力和温度下不需要重复调用 XSteam；二是把迭代循环改为向量化，多个加热器并列计算时，能一次算一组就一次算一组。

我项目里做的优化是：把 XSteam 每次调用放到一个带缓存的自定义函数里，用容器 Map 缓存相同输入参数的输出。实际测试下来，在 1000 次蒸汽参数调用的情况下，运行时间从 45 秒降到了 8 秒左右。这个优化非常简单但收益极其明显。

部署方面，我的程序最后打包为 MATLAB Compiler 生成的独立 exe，脱离了 MATLAB 环境也能运行。但注意 MATLAB Compiler 生成的安装包非常大，而且首次启动会释放运行时组件，在旧电脑上可能要一两分钟。如果只是内部自用，直接保留 MATLAB 脚本开发环境就够了，没必要打包。

最后再分享一个经验：热平衡程序的核心不是代码写得花哨，而是模型假设一定要在代码注释和数据源文档里写清楚。比如再热压损默认取 8%，管道压损默认取 5%，这些数值必须有出处，否则半年之后回头看自己写的程序，你根本没法判断结果里包含了哪些假设。对工程计算程序来说，结果可追溯比性能重要得多。

本文还有配套的精品资源，点击获取