简介：基于MATLAB的IAPWS-IF97水和水蒸气热力学性质计算包，面向能源、化工、制冷等领域需要获取蒸汽压、密度、焓、熵等参数的工程师与科研人员。包体共13个文件，以10个m函数/脚本为主，涵盖核心计算类、焓熵与密度导数测试及示例程序，并附有README说明文档、license授权文件等，压缩包仅111KB，结构简洁便于二次开发。已有336人学习下载。该资源封装了IF97标准算法，可完成从常规状态到超临界区域的热力学参数求解，适合在Matlab中开展蒸汽动力循环仿真、换热器设计或教学实验时作为底层工具；内置测试与示例脚本也有助于理解分区计算和边界条件处理思路。

1. 在 MATLAB 里写 IAPWS-IF97：热力计算为什么不能只查表

做发电、锅炉、汽轮机或者超临界二氧化碳循环仿真的人，迟早会碰到同一个问题：水和水蒸气的热力性质查表太慢，线性插值又不够准。我用 MATLAB 做过好几个热力系统的动态仿真模型，最开始的版本全靠查表插值，表格一细化，仿真慢到跑不动；表格一简化，过热区焓值能偏出十几个千焦。后来改用 IAPWS-IF97 标准直接计算，才把速度和精度同时稳住。IF97 是国际通用的水和水蒸气物性计算工业标准，覆盖 273.15 K 到 2273.15 K、压力最高到 100 MPa 的范围，从压缩水、湿蒸汽、过热蒸汽一直算到超临界区和高温区。MATLAB 里没有官方内置的 IF97 工具箱，所以绝大多数工程实践是把标准里的分区域方程自己实现成函数库。这篇文章就是把我的实现路径、分区判定逻辑、密度迭代方法和踩过的坑完整讲一遍，让新手能照着搭出第一版，让熟手能直接对照参数和边界细节做校验。

2. IF97 的五个分区：数学结构决定了代码怎么组织

拿到 IF97 标准的第一反应通常是去找公式，但真正决定代码结构的是分区的数学形式。IF97 把全部工况切成了五个区域，每个区域用不同的特征函数描述，有的直接给解析式，有的必须迭代求解。如果不先把分区逻辑理顺，后面写出的函数一定会在边界处翻车。

2.1 五个区的物理划分与基础方程类型

Region 1 覆盖低温液态水区，温度从 273.15 K 到 623.15 K，压力从饱和线向上到 100 MPa。Region 1 和 Region 2 都使用无量纲吉布斯自由能，以温度和压力作为独立变量，所以输入 T、p 可以直接算出焓、熵、比容，不需要任何迭代。Region 2 是蒸汽区和过热区，温度跨度很大，从 273.15 K 一直到 1073.15 K，而且为了拟合精度又被切成 2a、2b、2c 三个子段，每个子段有不同的系数组。Region 3 是最特殊的临界区，覆盖 623.15 K 到 863.15 K 之间靠近临界点的区域，像临界点附近这种温度和压力都能急剧变化的地方，以 T、p 为变量的吉布斯自由能形式拟合效果变差，标准改用亥姆霍兹自由能，并且把密度 ρ 作为独立变量。这意味着给到 T、p 时并不能直接求值，要先通过迭代解出密度，再算焓、熵等导出量。Region 4 是饱和线，也就是两相区的分界线，用来判断其他区域之间的边界。Region 5 是高温低压区，温度 1073.15 K 到 2273.15 K，压力上限 50 MPa，方程形式又回到吉布斯自由能。

2.2 为什么 Region 3 不能像 Region 1 那样直接算

理解 Region 3 需要迭代这件事，是设计 MATLAB 函数的第一步。Region 1 和 Region 2 的吉布斯自由能 g(p,T) 对压力求导就能得到比容，完全显式；但 Region 3 的亥姆霍兹自由能 f(ρ,T) 是密度的函数，压力是由自由能对密度求导得到的，实际工况一般给的是 (T,p)，所以要反过来解 f 的导数关系：

p = ρ² · (∂f/∂ρ) 在给定 T 下成立的那个 ρ

这就是一个非线性方程求解问题。我见过不少人在这个环节想偷懒，用饱和蒸汽密度或者临界密度直接当近似值用，结果在超临界区算出的焓值误差大到没法看。正确的做法是写一个牛顿迭代求解器，先把 ρ 的初值定好，再不断修正。初值选得好，迭代通常 3 到 6 步就能收敛到 1e-8 的相对精度；初值选得随意，临界点附近很容易发散。另一个数学上的重要细节是，Region 3 的自由能表达式中包含了理想气体项和残差项，两类项对压力计算的贡献在临界点附近会相互抵消，导致数值上出现“两个大数相减”的精度损失，这也是为什么不能简单用双精度直接莽过去。

Region 1 和 Region 2 的吉布斯自由能表达式在标准里是带 30 到 40 项系数的大多项式，Region 3 更复杂。写代码的时候不要尝试自己重新拟合系数，直接采用标准发布的数值表，把系数作为常量数组塞进函数里。下面这张表是我在实现时用来组织代码结构的依据：

这个表的重点是：只有 Region 3 需要密度迭代。因此 MATLAB 代码里的核心函数大致可以分成区域判定、Region 1/2/5 直接求值、Region 3 的密度求解、Region 4 饱和线计算四类。不要把这四类逻辑揉进一个巨型函数，否则后续换单位、改边界条件时会非常痛苦。

3. 从零搭 IF97 计算模块：区域判定函数与 Newton-Raphson 迭代密度

明确了分区结构之后，就可以开始写代码了。我的习惯是先做两个最基础的函数：一个是判定输入状态点落在哪个区，另一个是 Region 3 的密度迭代求解器。前者是所有后续计算的入口，后者是整套库里最容易出问题的部分。

3.1 先写区域判定函数：边界逻辑必须独立成函数

区域判定的输入是温度和压力，输出是区域编号。判定顺序很重要，我一般先处理 Region 5 和 Region 1 的范围，再处理 Region 2 和 Region 3 的边界。需要注意的一个细节是 Region 3 的上边界由 B23 方程定义，这是一条从 (623.15 K, 16.529 MPa) 延伸到 (863.15 K, 100 MPa) 的曲线，而不是简单的“温度大于 623.15 就归 Region 3”。下面这个函数是我常用的判定骨架：

function region = if97_region(T, p) % 判定 IAPWS-IF97 状态点所在区域 % 输入: T 温度(K), p 压力(MPa) % 输出: region 区域编号(1~5), 返回0表示超出范围 Tc = 647.096; % 临界温度 K pc = 22.064; % 临界压力 MPa T_min = 273.15; % 标准下限温度 % 先做整体范围检查 if T < T_min || T > 2273.15 || p < 0 || p > 100 region = 0; return; end % Region 5: 高温区 if T > 1073.15 if p <= 50 region = 5; else region = 0; % 高温高压超出标准覆盖 end return; end % Region 1: 低温液态水区 % 压力严格高于该温度下饱和压力，且温度不超过 623.15K psat = if97_psat(T); % 调用饱和压力方程，Region 4 if T <= 623.15 if p >= psat region = 1; return; end end % Region 3 判定取决于 B23 边界曲线 if T > 623.15 && T < 863.15 p_b23 = if97_b23_p(T); % B23 边界压力 if p > p_b23 region = 3; return; end end % 剩余情况: 低于饱和线的低压区归 Region 2 region = 2; end

这个函数的逻辑顺序有三个关键点。第一，Region 3 的判断必须放在 Region 2 前面，因为 B23 曲线下方仍然是 Region 2 的范围，不能简单用温度一刀切。第二，饱和压力计算 if97_psat 本身是一个独立的 Region 4 函数，如果 Region 1 的判断依赖它，那么模块的初始化顺序就是先实现饱和压力函数。第三，对于温度在 623.15 K 到 863.15 K 之间、压力低于 B23 的工况，代码会走到最后的region = 2，这是正确的。

3.2 Region 3 密度迭代：牛顿法的工程化写法

Region 3 的密度求解，本质是在给定 T 和 p 的条件下，找到使亥姆霍兹自由能导出压力等于输入压力的密度值。用 MATLAB 写最直接的牛顿迭代，形式如下：

function rho = if97_rho_region3(T, p, rho0) % 求解 Region 3 密度 % 输入: T 温度(K), p 压力(MPa), rho0 密度初值(kg/m3) % 输出: rho 密度(kg/m3) % 依赖: p_calc_region3(T, rho) 根据亥姆霍兹自由能计算压力 % 默认初值: 临界密度的 0.8 倍，工程上常用 if nargin < 3 || isempty(rho0) rho0 = 0.8 * 322; % 临界密度约 322 kg/m3 end rho = rho0; tol = 1e-10; % 相对压力误差容忍度 maxIter = 50; for k = 1:maxIter p_calc = p_calc_region3(T, rho); % 当前密度下的计算压力 resid = p_calc - p; % 残差，单位 MPa % 中心差分求 d(p_calc)/d(rho) drho = max(abs(rho) * 1e-6, 1e-6); % 差分步长 p_plus = p_calc_region3(T, rho + drho); p_minus = p_calc_region3(T, rho - drho); dp_drho = (p_plus - p_minus) / (2 * drho); % 防止导数过小导致步长爆炸 if abs(dp_drho) < 1e-6 error('if97:rhoSolver', '密度导数接近零，检查是否在临界点附近'); end drho_step = -resid / dp_drho; rho = rho + drho_step; if abs(drho_step) / rho < tol break; end end if abs(p_calc - p) / p > 1e-8 warning('if97:rhoSolver', '迭代未完全收敛，相对误差 %e', abs(p_calc - p)/p); end end

这里的 p_calc_region3 是核心物性函数，它把标准里 Region 3 亥姆霍兹自由能对密度求偏导后包装成压力的计算函数。我刻意用中心差分去算导数，而不是写解析导数，原因只有一个：Region 3 的表达式太复杂，解析求导容易在抄写时出错，中心差分配合 1e-6 的相对步长，精度已经足够让牛顿迭代收敛到 1e-8。差分步长不能取得太大，否则临界点附近会引入过多的截断误差。初值的选取上，不少工程代码会用标准推荐的区域 3 初始密度估计公式，我自己的习惯是先用临界密度附近的固定比例起步，收敛失败再换一个初值，工程够用。

3.3 调用示例与单位约束

写完求解器后，调用方式很简单。所有内部接口统一用 K 和 MPa，这是 IF97 标准公式里最常用的单位制，可以避免 Pa 与 MPa 混用导致数量级错误。焓和熵这类导出量在标准里是以 kJ/kg 为单位的，但中间计算时要注意无量纲化的因子还原。

% 示例: 计算 650 K, 25 MPa 的密度(超临界工况) T = 650; % K p = 25; % MPa rho = if97_rho_region3(T, p); fprintf('密度 = %.4f kg/m3\n', rho);

运行这段代码前，需要先实现 if97_region 和 p_calc_region3。如果 p_calc_region3 还空着，可以先用一个简单的理想气体压力近似占位，把迭代链路跑通，再逐步替换成标准方程。这种分层替换的做法能减少调试时的排查范围。

4. 区域判定与边界方程：B23 曲线与饱和线判据

区域判定是整个 IF97 模块的入口，也是最容易写错的部分。标准的边界曲线并不都是一条简单的等温线或者等压线，Region 2 和 Region 3 之间的 B23 边界、Region 4 的饱和压力方程，都有各自的数学形式。判定逻辑错了，后面所有计算都失去意义。

4.1 B23 方程：Region 2 与 Region 3 的分界曲线

B23 边界是把 (T, p) 平面分成 Region 2 和 Region 3 的关键曲线。标准给出的 B23 方程是一个以温度为自变量、压力为因变量的三次多项式，压力和温度都取 MPa 和 K。这条曲线从低温侧的 (623.15 K, 16.529 MPa) 出发，一路延伸到 (863.15 K, 100 MPa)。

function p_b23 = if97_b23_p(T) % B23 边界: 给定温度计算边界压力(MPa) % 适用范围: 623.15 <= T <= 863.15 K % 系数来自 IAPWS-IF97 标准附录(略写) n1 = 0.348051856289696e+3; n2 = -0.116718598799752e+1; n3 = 0.101929700393261e-2; n4 = -0.152879418971373e-5; p_b23 = n1 + n2 * T + n3 * T.^2 + n4 * T.^3; end

这段函数很短，但边界逻辑的影响面很大。如果 T 刚超过 623.15 K、p 在 16.529 MPa 以下，按 B23 判断会落在 Region 2，但实际上这个温度压力下的状态非常接近临界区的延伸，直接把方程类型切过去会导致焓值在边界面两侧出现不连续。我在某锅炉仿真项目里就碰到过：省煤器模型在 623 K 附近来回波动时，物性输出一跳一跳的，原因就是判定逻辑里用了小于号而不是小于等于，边界点被重复归类。后来把所有边界判断都改成单侧包含，跳变就消失了。

4.2 Region 4 饱和线：区域判定与两相区处理的争议

Region 4 饱和线的作用有两层。第一层是作为 Region 1 和 Region 2 的判定边界；第二层是当你需要计算两相区的湿度或者干度时，需要调用饱和压力方程。标准里的饱和压力方程是一个隐式形式，最常用的是反推式，也就是给定温度直接算饱和压力的高次多项式。

实际工程中一个常见困惑是：湿蒸汽两相区到底算哪个区？从 IF97 的角度看，Region 4 是区域之间的边界，本身不覆盖整个两相区内部。当你给一个压力和一个干度时，需要先用饱和线求出饱和温度和饱和焓，再按干度线性组合得到两相区的焓值。这种“干度加权”的做法是工程上常用的方式，我在第 3 章的区域判定函数中没有单独处理两相区，而是把低于饱和线的工况都归入 Region 2，让上层调用时再判断是否需要两相修正。

区域判定中有个容易忽略的点是临界点。647.096 K、22.064 MPa 在标准里恰好是 Region 1/2/3/4 的交汇处，在临界点附近，不同方程的数值行为都很极端，密度对压力的导数趋近无穷大，温度微小变化会引起密度大幅波动。实际编码时，我通常会在判定函数里显式加一个临界点附近的标记：

% 临界点邻近工况特殊处理 if abs(T - 647.096) < 0.01 && abs(p - 22.064) < 0.01 warning('if97:region', '临界点附近工况，物性计算误差可能增大'); end

这段预警代码不改变区域判定结果，但能在仿真跑到临界区时给出提示，帮我们快速把问题定位到物性模型上，而不是误以为控制算法出了 bug。

4.3 判定函数在实际仿真循环里的性能注意点

MATLAB 仿真的时间步进了区域判定函数的调用频率极高，如果每个时间步在循环里逐点调用完整判定函数，性能会很难看。常见的做法是向量化：把一批温度压力同时传入区域判定函数，让 MATLAB 的数组运算替代 for 循环。但 B23 方程和饱和压力方程里都有条件分支，向量化时要注意用逻辑索引分批处理，不能对整块数组直接套用同一个多项式。

我习惯把区域判定做成两个入口：一个是标量入口，给单点调试用；另一个是向量入口，内部用arrayfun或者纯逻辑索引做批量映射。这样做的好处是代码不重复，性能也能满足实时仿真需求。如果你的仿真步长很紧，还可以进一步把区域编号缓存下来，只有在温度压力变化超过一定阈值时才重新判定。

5. IF97 计算翻车的 5 个高频原因：从单位混乱到临界区发散

IF97 计算模块写出来容易，跑对难。我把这几年见过和踩过的坑总结成五条，每条都是“现象、原因、解决”的结构，方便直接对号入座。

5.1 单位制混乱导致压力相差三个数量级

现象：Region 1 算出的液体焓值正常，但到了 Region 2 的蒸汽状态点，焓值离谱偏小，密度偏大接近 1000 倍。

原因：标准里的压力在不同公式中有时用 MPa，有时用 Pa。我在早期把某个 Region 2 的系数数组直接照搬，却忘了输入压力从 MPa 转成了 Pa，导致多项式内压力项变成 10^6 倍量级。这类错误最隐蔽的地方在于：高温区的多项式对压力的依赖没有那么强，部分工况算出来的结果“看起来”还有合理性，只有对比到饱和线附近才暴露。

解决办法：所有对外接口统一用 MPa，函数内部任何用到标准公式的位置都先做一次单位换算并加注释。更稳妥的做法是把单位转换封装成独立函数，例如p_pa = p_mpa * 1e6，不要散落在各处。

5.2 临界区附近牛顿迭代发散

现象：Region 3 密度求解在 645 K、接近临界压力的工况下不收敛，迭代步数全部用尽后密度还在来回震荡，或者一步跳成负数。

原因：临界点附近 dp/dρ 趋近于零，残差函数变得非常平缓，牛顿时步长过大导致震荡。另外，如果中心差分步长 drho 固定使用绝对 1e-6，在密度只有 100 kg/m³ 级别的工况下，差分步长相对过大，导数计算失真。

解决办法：对差分步长改成相对步长，并限制最大步长。我在第 3 章的代码里用了max(abs(rho)*1e-6, 1e-6)，这能解决一部分问题，但还不够。临界点附近的工况还要加一个阻尼因子，让步长乘以一个 0.5 的衰减系数，防止一次跨出合理范围。

% 在迭代循环中加入阻尼 drho_step = -resid / dp_drho; drho_step = sign(drho_step) * min(abs(drho_step), rho * 0.3); rho = rho + drho_step;

把每一轮的最大修正量限制在当前密度的 30% 以内，收敛速度会慢一两个迭代，但稳定性好很多。温度越接近临界点，这个限制越重要。

5.3 Region 2 的子段边界被忽略

现象：过热蒸汽在跨越 2a 到 2b 的边界时，计算焓值出现约 0.1 kJ/kg 量级的跳变，虽然数值不大，但在求解器等熵线时会导致分叉或非物理震荡。

原因：Region 2 的吉布斯自由能表达式是分段定义的，2a、2b、2c 三个子段各自有独立的系数组。如果代码里只实现其中一段，或者实现了两段但边界判断条件是错误的直线而不是标准规定的温度边界，就会出现不连续。

解决办法：严格按照标准规定的子段边界划分，Region 2 内部的判定不需要用户干预，直接在求值函数里按照 T 的值选择对应子段即可。调试时可以用跨越边界的一串温度点检查焓值是否连续，这是最直观的验证。

5.4 Region 3 迭代初值离真实解太远

现象：同一套代码在 630 K、20 MPa 的工况下正常收敛，但换到 660 K、23 MPa，迭代直接发散或者收敛到一个非物理的高密度值。

原因：Region 3 的压力-密度曲线不是单调处处陡峭的，在某些区域可能存在多个密度满足同一个压力。牛顿迭代没有全局收敛性，初值落在错误的吸引域里就会收敛到非物理解。

解决办法：给每个典型工况准备一组初值，最常用的是把初值按饱和液密度和临界密度的加权组合来设定。工程上还有一个更稳的手段：用前一个时间步的密度作为当前时间步的初值，这在动态仿真中非常有效，因为相邻时间步的密度变化很小。静态计算时则用标准推荐的初始密度估计式。

5.5 向量化时误用条件分支导致数组维度爆炸

现象：把区域判定函数改写成向量版本时，某些输入数组返回的结果是错位的，例如温度数组和压力数组不匹配，或者结果是一列而输入是行向量。

原因：暴力使用if T > 623.15对数组求值会得到逻辑数组，对逻辑数组取 else 分支会导致维度丢失或错误索引。

解决办法：用逻辑索引分块计算，每一块先用布尔掩码选出符合条件的元素，计算完再放回原数组的对应位置。如果你不想在这上面花时间，最简单的方式是保持标量接口，在调用处用循环包装，性能和可读性都能接受，只是会慢一些。

6. 最后一道关口：用已知状态点给 IF97 实现做回归自检

代码写完不是终点，必须做一次精度回归。IF97 标准里有一组用于验证的参考状态点，这些点覆盖五个区域，可以用来检验你的实现有没有抄错系数。我给自己的代码做了一次“体检”，把常用的参考值列成表格，写了一个自动回归脚本，每次改动完系数都会重新跑一遍。

回归脚本的逻辑很简单，对每个点调用区域判定函数和密度求解，然后和预期值比较相对误差。如果误差超过 1e-6，就定位出问题区域。这里有个细节值得说一下：Region 1 和 Region 2 是显式计算，理论上精度应该和标准系数本身的精度一致；Region 3 因为经过了牛顿迭代，误差主要由迭代容忍度决定，所以我把压力残差的容忍度设置到 1e-10，以保证最终焓值的稳定。

% 回归自检脚本 testPoints = [ 300, 3, 998.9; 500, 0.5, 2.175; 900, 1, 3.948; 650, 30, 526.9; 1500, 1, 0.456 ]; for i = 1:size(testPoints, 1) T = testPoints(i, 1); p = testPoints(i, 2); rhoExpected = testPoints(i, 3); rho = if97_density(T, p); % 统一的密度求值入口 relErr = abs(rho - rhoExpected) / rhoExpected; fprintf('点 %d: 计算密度 %.4f, 参考密度 %.4f, 相对误差 %.2e\n', ... i, rho, rhoExpected, relErr); end

这个入口函数if97_density内部按区域分发：Region 1/2/5 直接由吉布斯自由能对压力求导得到比容再取倒，Region 3 走迭代求解，Region 4 则看调用方传入的是压力还是干度。分发逻辑清晰后，任何上层业务代码都只需要面对一个函数接口，不需要关心区域细节。

关于精度还有一个容易被忽视的因素：MATLAB 的浮点数运算在临界点附近会因为两个大数相减而损失有效位数。解决这个问题的方法是把 Region 3 自由能表达式里的理想气体项和残差项分别计算后合并，不要展开成单一多项式。这个调整能把临界点附近的密度误差从 1e-5 改善到 1e-7 量级，代价是代码稍微长一点。我现在的习惯是：Region 3 的实现严格保持“分项”结构，并在注释里标注哪些项来源于理想气体贡献，哪些是残差项。

最后一个工程经验：不要试图自己“优化”标准方程的形式来加快计算。IF97 的系数本身就是经过数值稳定性设计的，任何重新合并同类项、近似替换的做法，都会在某一个不明显的工况下让精度突然崩掉。我在一个余热回收模型里曾经为了让 Region 2 计算提速，把其中一项做了泰勒展开截断，结果在 850 K、20 MPa 附近出现了约 0.5% 的焓值偏差，后来花了很长时间才定位到是那次“优化”引入的。从那以后，我对已经跑通的标准实现敬而远之，只做调用层优化，不动方程层。这个教训值得每个做物性库的人记住。希望帮到你。

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