天然气BWRS方程求解器：压缩因子计算与参数表详解

简介：这是面向石油天然气行业技术人员和高校相关专业学生的BWRS方程C++实现，用于计算高压低温条件下天然气的压缩因子与密度，可支撑管道设计、流量计算及储气设施容量评估。资源共1个文件，为一份cpp源程序，压缩包大小约2KB，代码结构完整，包含输入压力、温度及摩尔组成后，采用牛顿法等数值方法迭代求解摩尔体积，进而计算压缩因子Z和密度ρ。已有579人学习。通过阅读和运行该程序，读者可省去从零推导BWRS公式的繁琐过程，直接获得可编译的参考代码；也能基于其模块化写法调整组分和参数，用于工程验算、教学演示或二次开发，兼顾精度与实现效率。

1. 当PR方程把压缩因子算偏时，天然气BWRS方程就是那粒后悔药

做天然气计量的工程师多半都经历过这种时刻：交接报表上标况体积差了0.8%，下游单位拿着流量计检定证书来对质，压力变送器和温度探头都查过没问题，最后翻出来是压缩因子Z取错了。PR方程算得快，但有些高压富气工况它就是偏，换用天然气BWRS方程之后，模型剩余残差能降一个量级。这份资源就是把BWRS方程的完整计算链路打包好：11个纯组分参数、混合规则、Newton迭代求Z、焓熵修正，适合做管输工艺、储气库方案、计量比对的人直接拿去跑，也适合想在PR方程之外多备一套状态方程的开发者。

2. BWRS方程的底细：11个参数、混合规则和适用边界

2.1 从BWR到BWRS：温度项和第三维里项多了什么

BWR方程在20世纪40年代提出时只有8个参数，主要面向烃类蒸气压和密度的工程计算，但它在低温区域对温度变化的描述太粗糙，温度项只到C0/T²就截断了。BWRS方程把温度多项式往后多延了两级，引入D0/T³和E0/T⁴，同时把第三维里系数从「bRT − a」扩展成「bRT − a − d/T」，高次项从「a·α」改成「a(a + d/T)」。这样改的直接收益是：在天然气储运常见的240K到330K温区里，密度的温度外推行为不再是一条被T²绑死的曲线，而是能贴合实际液化点和近临界区的曲率变化。

方程的标准形式如下，这里ρ是摩尔密度，单位用kmol/m³，压力单位用kPa：

p = ρRT + (B0·R·T − A0 − C0/T² + D0/T³ − E0/T⁴)·ρ² + (b·R·T − a − d/T)·ρ³ + a·(a + d/T)·ρ⁶ + (c·ρ³/T²)·(1 + γ·ρ²)·exp(−γ·ρ²)

逐项看，ρRT是理想气体贡献；ρ²项代表分子对相互作用，系数里的B0RT是排斥贡献，A0、C0、D0、E0是不同温度尺度下的吸引贡献叠加；ρ³项对应三分子碰撞，其中d/T让它在低温下吸引效应更强；ρ⁶项是经验性的高密度修正项，a(a + d/T)保证高压下密度增大时压力不会发散；最后一项是指数衰减项，只在较高密度下起作用，负责近临界区等温线的弯曲形态。整套参数里A0、B0、C0、D0、E0、a、b、c、d共10个出现在状态方程中，第11个参数e不直接进压力方程，它出现在焓、熵残差函数的温度导数项里。所以如果你看到参数表里有e而状态方程里没有，别觉得表错了，它在算热力学性质时才出场。

2.2 混合规则与交互系数：纯组分参数怎么合到一起

纯组分的11个参数不能直接用于混合物，工程实现中用的是一套固定的混合规则。A0项由于和分子对相互作用能挂钩，采用二次交叉混合，并留了一个二元交互系数k_ij给后场调节；其余参数（B0、C0、D0、E0、a、b、c、d、e）按立方根规则混合，γ按摩尔分数线性加权。公式如下：

A0_m = ΣΣ x_i · x_j · (A0_i · A0_j)^0.5 · (1 − k_ij) B0_m = Σ x_i · B0_i C0_m = (Σ x_i · C0_i^(1/3))^3 D0_m = (Σ x_i · D0_i^(1/3))^3 E0_m = (Σ x_i · E0_i^(1/3))^3 a_m = (Σ x_i · a_i^(1/3))^3 b_m = (Σ x_i · b_i^(1/3))^3 c_m = (Σ x_i · c_i^(1/3))^3 d_m = (Σ x_i · d_i^(1/3))^3 e_m = (Σ x_i · e_i^(1/3))^3 γ_m = Σ x_i · γ_i

这里x_i是摩尔分数。需要注意k_ij的取值：同族烃类之间k_ij通常取0或极小值（比如甲烷和乙烷可以当成0），但含CO2、N2的时候不能简单取0。以甲烷和CO2为例，若不设交互系数，高压力下混合物的密度会偏大，压缩因子偏小，直接造成管存计算偏差。资源包里如果带交互系数表，优先用包内数据；如果没有，我一般先按惯用值放一组进去，再用第4章的标定流程重调。

2.3 适用边界：温度、压力和组分窗口

用BWRS方程前先判断工况是否落在它的甜区里，不然算出来再准也是假的。我用一张表给边界做个整理，数值是这类状态方程文献里常见的经验范围，具体以你手上的参数表来源为准：

这里说的“不适合”不是不能算，而是误差会超过工程可接受范围。比如水含量高的湿气，分子间氢键作用会让BWRS的经验参数失效，计算密度和实测偏差很大。C7+重烃的问题更具体：BWRS参数表通常只覆盖到庚烷附近，更重的组分没有可靠参数，硬用虚拟组分拟合会引入额外误差。写完这套边界判断，再进入代码复现环节，心里就有底了。

3. 手写天然气BWRS方程求解器：参数表、Newton迭代与选根

3.1 参数表结构和单位体系：先让数据能对上

BWRS方程复现第一步不是写迭代，而是先把单位体系统一好。常见的坑是参数表用英制单位（psia、°R、lbmol/ft³），代码里却按公制算，出来的压力差好几个数量级。我处理的思路是先把参数表转成一套自洽的公制单位：压力kPa、温度K、摩尔密度kmol/m³、气体常数R取8.314 kPa·m³/(kmol·K)。注意这里的R不带质量，是摩尔基准，和参数表里B0、b这些项的量纲是一一对应的。

参数表用CSV存，每行一个纯组分，列名按固定顺序：A0, B0, C0, D0, E0, a, b, c, d, e, gamma。读取代码大致这样：

import csv def load_bwrs_params(path): params = {} with open(path, 'r', encoding='utf-8') as f: reader = csv.DictReader(f) for row in reader: name = row['component'] params[name] = { 'A0': float(row['A0']), 'B0': float(row['B0']), 'C0': float(row['C0']), 'D0': float(row['D0']), 'E0': float(row['E0']), 'a': float(row['a']), 'b': float(row['b']), 'c': float(row['c']), 'd': float(row['d']), 'e': float(row['e']), 'gamma': float(row['gamma']) } return params

这段逻辑很简单，但有两点值得说。第一，CSV表头里最好带component这一列，方便后续扩展组分；第二，读进来之后统一用float转换，避免字符串参与运算时Python抛TypeError。拿到参数表后，写一个混合函数，按第2章公式把纯组分参数合成混合参数：

def mix_params(params, x, kij_func=None): names = list(params.keys()) comps = [params[n] for n in names] n = len(names) A0 = 0.0 for i in range(n): for j in range(n): kij = 0.0 if kij_func: kij = kij_func(names[i], names[j]) A0 += x[i] * x[j] * (comps[i]['A0'] * comps[j]['A0']) ** 0.5 * (1 - kij) B0 = sum(x[i] * comps[i]['B0'] for i in range(n)) C0 = sum(x[i] * comps[i]['C0'] ** (1/3) for i in range(n)) ** 3 D0 = sum(x[i] * comps[i]['D0'] ** (1/3) for i in range(n)) ** 3 E0 = sum(x[i] * comps[i]['E0'] ** (1/3) for i in range(n)) ** 3 a = sum(x[i] * comps[i]['a'] ** (1/3) for i in range(n)) ** 3 b = sum(x[i] * comps[i]['b'] ** (1/3) for i in range(n)) ** 3 c = sum(x[i] * comps[i]['c'] ** (1/3) for i in range(n)) ** 3 d = sum(x[i] * comps[i]['d'] ** (1/3) for i in range(n)) ** 3 e = sum(x[i] * comps[i]['e'] ** (1/3) for i in range(n)) ** 3 gamma = sum(x[i] * comps[i]['gamma'] for i in range(n)) return {'A0': A0, 'B0': B0, 'C0': C0, 'D0': D0, 'E0': E0, 'a': a, 'b': b, 'c': c, 'd': d, 'e': e, 'gamma': gamma}

注意A0循环里用了双重求和，n=5时是25次，效率不是问题。kij_func传入一个回调，返回两个组分名的交互系数，这样交互系数表可以单独维护，不需要改混合函数本体。这是我在实际工程里比较喜欢的组织方式：参数表只管纯组分，交互系数单独建表，混合函数保持纯净。

3.2 压力函数与解析导数：Newton迭代的骨架

给定温度T、密度ρ、混合参数prm，压力函数直接按第2章公式翻译成代码。同时要把导数也写出来，Newton迭代用解析导数收敛更快，数值差分虽然省事，但在高密度区步长选不好容易把指数项抖出噪声。我一般两边都留，默认走解析导数：

import math R = 8.314 # kPa·m3/(kmol·K) def bwrs_pressure(T, rho, prm): A0 = prm['A0']; B0 = prm['B0']; C0 = prm['C0'] D0 = prm['D0']; E0 = prm['E0'] a = prm['a']; b = prm['b']; c = prm['c'] d = prm['d']; gamma = prm['gamma'] r = rho term1 = r * R * T term2 = (B0 * R * T - A0 - C0 / T**2 + D0 / T**3 - E0 / T**4) * r**2 term3 = (b * R * T - a - d / T) * r**3 term4 = a * (a + d / T) * r**6 term5 = (c * r**3 / T**2) * (1 + gamma * r**2) * math.exp(-gamma * r**2) return term1 + term2 + term3 + term4 + term5 def bwrs_dpdrho(T, rho, prm): A0 = prm['A0']; B0 = prm['B0']; C0 = prm['C0'] D0 = prm['D0']; E0 = prm['E0'] a = prm['a']; b = prm['b']; c = prm['c'] d = prm['d']; gamma = prm['gamma'] r = rho term2 = 2 * (B0 * R * T - A0 - C0 / T**2 + D0 / T**3 - E0 / T**4) * r term3 = 3 * (b * R * T - a - d / T) * r**2 term4 = 6 * a * (a + d / T) * r**5 g = gamma u_prime = (3 * r**2 + 3 * g * r**4 - 2 * g**2 * r**6) * math.exp(-g * r**2) term5 = (c / T**2) * u_prime return R * T + term2 + term3 + term4 + term5

这段代码里term5的导数我手推过，u = r³(1 + γr²)e^(−γr²)，对r求导展开后中间几个γr⁴项会抵消，最后落在3r² + 3γr⁴ − 2γ²r⁶这个形状上。如果你用数值差分验证，能对上说明导数没问题。注意term4里a(a + d/T)的d和term3里的d是同一个参数，写代码时不要拆成两个变量。

有了压力函数和导数，Newton迭代本体就很短了：

def solve_rho(T, p, prm, rho0, tol=1e-10, max_iter=100): rho = rho0 for it in range(max_iter): f = bwrs_pressure(T, rho, prm) - p df = bwrs_dpdrho(T, rho, prm) drho = f / df rho -= drho if abs(drho) < tol * rho: return rho, it raise RuntimeError(f'Newton not converged at T={T}, p={p}, rho0={rho0}')

收敛判据用的是密度相对变化量tol * rho，不是压力残差绝对值。原因是压力单位kPa量级差很大，低压下残差1Pa就算很多了，高压下100Pa也不算事，用相对密度变化做判据更稳。迭代上限给100次，正常5到10次就收敛，到50次还没收敛基本就是初值问题，再来100次也没用，直接抛异常比死循环好排查。

3.3 汽液双根与初值策略：别让迭代选错根

天然气压缩因子计算里最隐蔽的问题，是同一组T、p下方程组可能有多于一个的物理根。低压高温区只有一个气相根，但高压低温区会出现三个根：液相根、不稳定根、气相根。Newton迭代从不同初值出发，会收敛到不同的根，而工程上你要的是哪个根完全取决于工况判断。

我常用的策略是用两个初值分别迭代：气相初值用理想气体密度ρ_gas = p/(RT)，液相初值用一个经验密度，比如15到20 kmol/m³（温度越低取越大，240K时可以取到25）。这两个初值跑完会得到两个候选根，然后比较哪个根更符合物理直觉。对于纯气相管输，直接取气相根；对于管存计算、储气库注采，往往涉及两相区，需要进一步判断：

def solve_roots(T, p, prm): rho_gas = p / (R * T) try: rho_v, _ = solve_rho(T, p, prm, rho_gas) except RuntimeError: rho_v = None try: rho_l, _ = solve_rho(T, p, prm, 18.0) except RuntimeError: rho_l = None return {'vapor': rho_v, 'liquid': rho_l}

判断哪个根是稳定相态的常用做法是比较吉布斯自由能，但工程上有个更糙但有效的简化：把两个根各自反带回bwrs_pressure，看哪个根的残差更小，同时结合工况判断。如果气相根和液相根算出来的压力都能对上目标压力（说明两个都是数学解），那就按工艺状态选。这个环节没有一劳永逸的自动判断，和PR方程拿三次方程解析解就能选根不同，BWRS的指数项让选根必须保留人工确认的口子。这也是我把选根逻辑单独抽成函数的原因，后续接界面调用时，允许外部传入phase参数，强制指定取气相还是液相根。

4. 误差分析与交互系数标定：把计算值拉到贴近实测

4.1 最小二乘标定kij：三条实测点就能调一版

交互系数k_ij是BWRS混合规则里唯一可以自由调节的旋钮。它没有严格的理论推导值，来源是回归实验数据。工程上最常见的场景：你拿到一组现场取样组分和对应的实测密度（比如在线色谱仪和流量计配套的历史数据），发现用默认kij算出来的密度系统性偏小，这时就该启动一次小规模标定。

原理很简单：定义一个目标函数，让所有实验点上的计算压力或密度与实测值的偏差平方和最小。下面这段代码用scipy的least_squares做回归，实验数据格式是每行一个点（T, p, x数组, 实测密度rho_exp）：

import numpy as np from scipy.optimize import least_squares def fit_kij(exp_data, params, comps, kij0): # exp_data: list of dict, each with T, p, x(list), rho_exp def residuals(kij_vec): # 只调节第一个交互系数，假设是CH4-CO2 errs = [] for pt in exp_data: prm = mix_params(params, pt['x'], kij_func=lambda a, b: kij_vec[0] if {a,b} == {'CH4','CO2'} else 0.0) rho, _ = solve_rho(pt['T'], pt['p'], prm, pt['p'] / (R * pt['T'])) errs.append((rho - pt['rho_exp']) / pt['rho_exp']) return errs res = least_squares(residuals, kij0) return res.x

注意residuals里return的是相对误差而不是绝对误差，这样做的好处是把高压和低压点的权重拉平，避免低压点因为密度绝对值小而被淹没。least_squares默认的损失函数是平方，基本够用。三点经验：第一，标定数据至少要覆盖高压、中压、低压三个压力段，只在一个压力点回归出来的kij拿到别的压力段会翻车；第二，kij的变化范围一般限定在0到0.2之间，超过这个范围说明参数表本身有问题，别硬调；第三，每调一个kij，都要回到第3章的选根逻辑重新确认一下，因为kij变化会影响液相根的稳定性。

4.2 灵敏度查看：温度、压力、组分各自的影响路径

调完参数后要做一次灵敏度排查，搞清楚计算结果的误差到底从哪来。我用一张表来总结BWRS方程在这三个方向上的典型行为，方便你对照自己的工况判断：

这轮排查的重点不是算得多精细，而是确定你的输入数据里哪个量最不可靠。比如站场压力变送器的检定精度是±0.5%，那压力方向的误差传导就没必要花力气去抠；如果组分分析仪隔几天就漂一次，那才是误差的主要来源。我习惯在代码里把温度、压力、组分三个输入分别加一个微小扰动，各算一遍密度，观察输出的变化幅度，这比凭感觉找原因高效得多。

4.3 和商业物性库对表的步骤

做完标定以后，最后一个验证步骤是拿BWRS计算结果和商业物性库或者权威数据表做交叉对比。对表的目的不是要求两个软件完全一致，而是确认偏差在合理范围，并且偏差的方向是系统性的还是随机的。我推荐按这套流程走：

第一步，选一组标准组分，比如摩尔分数CH4 0.9、C2H6 0.05、C3H8 0.03、N2 0.02，固定不变。第二步，在250K、270K、290K三个温度下从1MPa到10MPa各取5个压力点，用你的脚本算出密度和Z。第三步，把同一组工况输入商业物性库（或者用公开的实验数据表），记录结果。第四步，把两组数做差值，如果Z的相对偏差都小于0.5%，说明你的复现没有原则性问题；如果偏差集中在某一压力段，回头查交互系数；如果全部偏差符号相同，先怀疑单位体系。

这一章做完，你手里的BWRS代码就从一个“能跑”的状态升级成“能用于工程判断”的状态。接下来聊聊我实际使用中最容易翻车的几个场景。

5. 常见问题排查：错根、发散、C7+和单位混用

5.1 压缩因子突然翻倍：根选错了

现象：某一组工况算出来的Z从0.85变成1.7，密度直接掉一半，而且只发生在压力升到某个值之后。原因：压力升高后方程组出现了第三个根，Newton迭代从气相初值出发时被“吸引”到了不稳定根或液相根的另一个分支上。解决：把压力函数画出来，看给定T下等温线的形状。我一般会在选根函数里加一个密度扫描兜底：从ρ=1e-6开始逐步增大到30 kmol/m³，每步做一次Newton迭代，把收敛结果按ρ排序，然后取与理想气体密度同侧的第一个根作为气相根。这个兜底虽然多花几十次迭代，但能避免报表上出现0.8%的交接偏差。

5.2 迭代发散的“玄学”时刻：初值太糙

现象：同一个求解器，250K时一切正常，240K时第一次迭代后密度变成负数，报错退出。原因：低温下指数项exp(−γρ²)在ρ较大时衰减很快，导数项接近零，Newton步长过大直接把密度推到了物理边界外。解决：给Newton迭代加一个阻尼因子，每步的步长限制在当前密度值的50%以内，超了就按边界走：

step = f / df max_step = 0.5 * abs(rho) if abs(step) > max_step: step = math.copysign(max_step, step) rho -= step

这里max_step取0.5倍当前密度，防止迭代点跨过零点和越过液相高密度区。加了阻尼之后收敛速度会慢两三轮，但稳定性提升明显，这段代码值得写进你的求解器里。

5.3 C7+轻烃不好处理：参数表根本覆盖不到

现象：组分表里C7+的摩尔分数给了2%，你从参数表里找不到对应行，硬把C7+的物性按C7H16近似，结果计算密度和现场实测差了3%。原因：BWRS参数表通常只回归到庚烷、辛烷附近，C7+是一个混合虚拟组分，其平均分子量和沸点都和纯庚烷不一致，直接替换是不合理的。解决：资源包内一般会附带重馏分特性化说明，按平均沸点和密度把C7+拆解成若干个虚拟单组分，再分别查参数。如果包内没有这层处理，我的建议是换用PR方程加重馏分特性化，别在BWRS上硬撑。

5.4 和PR方程结果差太多：先查交互系数再查代码

现象：同一组工况，BWRS算的密度比PR方程大2%，你觉得是BWRS错了。原因：大部分情况下不是方程本身的问题，而是交互系数没设置或单位体系不统一。排查顺序应该是这样的：

第4步是最快定位问题的手段：把bwrs_pressure里的term1到term5分别打印出来，正常情况下term2和term3占主导，term4和term5只有在液相高密度区才会显著。如果term5在气相工况下大得离谱，那就是指数项里的γ参数读错了行。这一条排查顺序我每次上线BWRS计算服务都要过一遍，能省下大量排查时间。

6. 快速验证与批处理：把BWRS算成一组管线数据

6.1 用一个实测点验证代码正确性

拿到这套BWRS资源后，我建议你做的第一件事不是直接上管线数据，而是先找一个纯甲烷的实测点验证代码整体正确性。纯甲烷的实验数据最完整，也最不容易被组分误差干扰。选一个典型工况，比如T=300K、p=5MPa，纯甲烷的压缩因子应该在0.93附近（具体值以你参考的数据表为准）。跑一下密度，反算Z，如果偏差超过2%，先检查单位体系；如果偏差在0.5%内，再进混合组分测试。

这一步花不了几分钟，但能帮你确认参数表、混合规则、求解器三个模块是通的。我自己吃过亏：拿到参数表后跳过了单组分验证，直接上混合气数据，结果计算偏差大，来回排查了两天才发现是参数表里有一行数据多了个小数点。从那以后，我每次拿到状态方程代码都强制先跑一遍纯甲烷单点验证，再进管线数据。希望帮到你。

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