简介这是一份面向油气储运、管道设计及化工工艺从业者的天然气物性计算工具采用程序设计语言实现了BWRS状态方程用于精确求解高压工况下的压缩因子和密度。代码基于四位学者提出的状态方程模型修正理想气体假设在高压低温场景下的偏差可服务于管道输送计算、储气库容量评估及流量计量等实际需求。程序完整覆盖初始化、摩尔体积非线性迭代、压缩因子换算与密度计算流程采用牛顿法或二分法迭代求解方程用户输入压力、温度和气体组成即可获得对应结果。资源包内共包含一个源代码文件整体体积约两千字节文件结构简单、便于阅读和移植。已有五百七十九人浏览学习适合油气行业工程师、科研人员以及希望深入理解BWRS方程编程实现的高校学生使用。1. 天然气物性计算先别急着套 PR 方程做天然气集输、储气库或深冷工艺设计时压缩因子 Z、密度、焓熵偏差这些参数直接决定后续的压缩机功率、换热面积和管径核算。不少人在前期用 PR 方程算了个大概进入高压工况或含重烃气样后密度误差能到百分之十几焓值更是偏离得没法用。本文要讲的是天然气领域里特别能打的 BWRS 方程全称 Benedict-Webb-Rubin-Starling 状态方程它用一个带指数项的十一参数方程把气相、液相和超临界区一起覆盖在宽温度宽压力范围内密度和焓熵计算都更贴近实验数据。文章从方程结构、参数来源和混合规则讲起再用 Python 一步步实现压缩因子与焓熵偏差的计算最后给出收敛、单位和参数匹配上常见的踩坑记录。适合正在写气体物性模块、做工艺包或要填热力学接口的工程人直接照着复现。2. BWRS 方程的数学结构多一项少一项差在哪2.1 从 BWR 到 BWRS指数项为什么被保留BWRS 方程的本质是建立在密度幂级数上的显式压力表达式。一般形式写作P ρRT (B₀RT - A₀ - C₀/T² D₀/T³ - E₀/T⁴)ρ² (bRT - a - d/T)ρ³ α(a d/T)ρ⁶ (cρ³/T²)(1 γρ²)exp(-γρ²)这里面 ρ 是摩尔密度T 是开尔文温度R 是通用气体常数。对照就能发现它比原始 BWR 方程多加了 D₀、E₀ 两个与温度高次幂相关的项还把原来的单一指数项扩展成带 γρ² 的指数修正。之所以这么干是因为原始 BWR 方程在低温重烃区和临界点附近的密度预测经常偏大加进 D₀、E₀ 后方程对温度倒数的高阶敏感度提升气液相边界拟合得更平滑。指数项的作用也不只是装饰它在高密度区压制多出的六次方项避免方程在液相密度下产生不合理的压力反弹。实际工程中BWRS 常用于温度范围约 -160℃ 到 200℃、压力最高 30~50MPa 的天然气体系。相比 PR 方程它最大的优势是密度计算在液相区和超临界区仍然稳定不会出现 PR 方程那种一到高密度就明显偏离实验值的毛病。代价也很明确参数多、混合规则复杂、求根迭代时对初始区间敏感。不能把它当黑匣子用至少要知道每个参数贡献在哪一段密度区间。2.2 十一参数表与混合规则参数一致性是底线BWRS 方程有十一个单组分参数A₀、B₀、C₀、D₀、E₀、a、b、c、d、α、γ。它们各自的物理意义大致可以这样理解A₀、B₀ 是最基础的分子吸引和排斥贡献C₀、D₀、E₀ 描述温度依赖的校正a、b、c、d 对应高次密度项α 和 γ 则控制指数修正项的形态。不同文献给出的甲烷参数会略有差别因为拟合时使用的实验数据温度压力窗口不同这不是谁对谁错而是适用范围不同。以甲烷为例在 kPa、L、mol、K 单位体系下的一套常用拟合值如下参数甲烷拟合值演示用说明A₀2.3921e2kPa·(L/mol)²B₀4.5200e-2L/molC₀2.2160e4kPa·(L/mol)²·K²D₀5.0400e3kPa·(L/mol)²·K³E₀0.0甲烷有时取 0a1.1010e-2kPa·(L/mol)³b1.5760e-3(L/mol)²c2.6200e3kPa·(L/mol)³·K²d1.5800e2kPa·(L/mol)³·Kα1.0100e-1无量纲γ7.4000e-2(L/mol)²看到不同来源值相差 1% 到 3% 不必惊讶选表的原则是整套参数必须来自同一次拟合或同一个参数库不能 A₀ 从这篇文献拿、C₀ 从那本手册抄。BWRS 方程参数与混合规则是成套的混搭会让结果变得没法解释。混合规则直接决定多组分天然气能不能算对。常见做法是B₀ 和 E₀ 用摩尔分数线性平均A₀、C₀、D₀ 用双重重叠的几何平均并引入二元交互系数 kᵢⱼa、b、c、d、α 用三分之一次方平均后三次方γ 用平方根平均后平方。写成公式就是A₀ ΣᵢΣⱼ xᵢxⱼ(A₀ᵢA₀ⱼ)^0.5(1 - kᵢⱼ)B₀ Σᵢ xᵢB₀ᵢC₀ ΣᵢΣⱼ xᵢxⱼ(C₀ᵢC₀ⱼ)^0.5(1 - kᵢⱼ)³D₀ ΣᵢΣⱼ xᵢxⱼ(D₀ᵢD₀ⱼ)^0.5(1 - kᵢⱼ)³a [Σᵢ xᵢaᵢ^(1/3)]³γ [Σᵢ xᵢγᵢ^0.5]²二元交互系数 kᵢⱼ 通常较小甲烷与乙烷、丙烷之间一般在 0.005 到 0.03 量级。对于没有实验回归数据的体系取 0 也能算出可用的气相性质但液相计算误差会明显变大。因此对含凝析油的天然气kᵢⱼ 一定要通过实验数据或可靠文献值补齐这是 BWRS 工程应用里被忽略最多的一环。3. 用 Python 实现 BWRS 压缩因子求解单组分到混合气3.1 单组分最小求解器二分法求密度压缩因子 Z 的定义是 P/(ρRT)因此只要从状态方程中解出密度 ρZ 就顺理成章得到。这里不推荐对密度直接做牛顿迭代因为 BWRS 方程在临界区和液相区容易出现平台段导数接近零牛顿法经常直接越过根。更稳的做法是用二分法先确定一个压力包络区间再缩小区间。以下是我常用的最小可运行实现import math # 单位约定T[K], P[kPa], rho[mol/L] # R 取 kPa·L/(mol·K)与参数表单位保持一致 R 8.314462618 # 甲烷 BWRS 单组分参数演示用正式项目先校核参数来源 CH4 { A0: 2.3921e2, B0: 4.5200e-2, C0: 2.2160e4, D0: 5.0400e3, E0: 0.0, a: 1.1010e-2, b: 1.5760e-3, c: 2.6200e3, d: 1.5800e2, alpha: 1.0100e-1, gama: 7.4000e-2, } def bwrs_p(rho, T, prm): 已知密度和温度用 BWRS 方程计算压力。 g prm[gama] p rho * R * T p (prm[B0] * R * T - prm[A0] - prm[C0] / T**2 prm[D0] / T**3 - prm[E0] / T**4) * rho**2 p (prm[b] * R * T - prm[a] - prm[d] / T) * rho**3 p prm[alpha] * (prm[a] prm[d] / T) * rho**6 p prm[c] * rho**3 / T**2 * (1.0 g * rho**2) * math.exp(-g * rho**2) return p def solve_density(T, P, prm, liquidFalse): 二分法求密度。liquidTrue 时偏向液相大密度根。 rho_ig P / (R * T) if liquid: rho_lo, rho_hi rho_ig, 40.0 else: rho_lo, rho_hi 1e-9, max(rho_ig * 4.0, 0.5) # 确保右端点计算压力已超过目标 P否则不断拓宽上界 f_hi bwrs_p(rho_hi, T, prm) - P while f_hi 0 and rho_hi 500: rho_hi * 2.0 f_hi bwrs_p(rho_hi, T, prm) - P if f_hi 0: raise RuntimeError(搜索区间上限仍低于目标压力) # 二分迭代 for _ in range(300): mid 0.5 * (rho_lo rho_hi) f_mid bwrs_p(mid, T, prm) - P if f_mid 0.0: return mid if f_mid 0: rho_hi mid else: rho_lo mid return 0.5 * (rho_lo rho_hi) # 单组分算例300K、5MPa 甲烷 T, P 300.0, 5000.0 rho solve_density(T, P, CH4) Z P / (rho * R * T) print(frho {rho:.5f} mol/L, Z {Z:.5f})这段代码里最关键的是 bwrs_p 函数的单位一致性R 取 8.314462618压力参数 A₀、C₀、D₀、c、d 都按 kPa 体系填写。solve_density 里先算理想气体密度 rho_ig用它作为二分区间下界的基准再把上界按 4 倍扩展。如果右端点压力仍小于目标压力就循环翻倍 rho_hi这个兜底逻辑能避免高压工况下因初始上界不足而直接报错。运行后会发现 Z 在 1 附近5MPa 下的甲烷压缩因子大约 0.9 出头符合工程直觉。真正要留意的是 liquid 参数当低温液化或高压凝析时状态方程在同一个压力下可能存在多个密度根气相根的密度小、液相根的密度大默认二分法只会收敛到区间内的某个根不区分物理分支。所以在明确知道是液相时要把 liquid 置为 True让搜索区间覆盖到大密度区域。3.2 多组分天然气混合规则代入按摩尔分数拼装参数实际天然气不是纯甲烷而是甲烷、乙烷、丙烷、氮气、二氧化碳等组分的混合物。多组分计算首先要把单组分参数按混合规则合成一组等效参数。下面这段代码实现常见的 BWRS 混合规则并把二元交互系数矩阵作为可选输入def mix_params(comps, z_frac, kij_matrixNone): comps: 单组分参数字典列表 z_frac: 摩尔分数列表长度与 comps 一致 kij_matrix: n x n 二元交互系数矩阵没有则全取 0 n len(comps) if kij_matrix is None: kij_matrix [[0.0] * n for _ in range(n)] m {} # B0 与 E0 按线性平均 m[B0] sum(z_frac[i] * comps[i][B0] for i in range(n)) m[E0] sum(z_frac[i] * comps[i][E0] for i in range(n)) # A0、C0、D0 用几何平均并引入二元交互系数 for key in [A0, C0, D0]: m[key] 0.0 for i in range(n): for j in range(n): kij kij_matrix[i][j] base math.sqrt(comps[i][key] * comps[j][key]) power 1.0 if key ! A0: power (1.0 - kij) ** 3 else: power 1.0 - kij m[key] z_frac[i] * z_frac[j] * base * power # a、b、c、d、alpha 按 1/3 次方平均再三次方 for key in [a, b, c, d, alpha]: s 0.0 for i in range(n): s z_frac[i] * (comps[i][key] ** (1.0 / 3.0)) m[key] s ** 3 # gamma 按 1/2 次方平均再平方 s 0.0 for i in range(n): s z_frac[i] * math.sqrt(comps[i][gama]) m[gama] s ** 2 return m # 二元交互系数矩阵甲烷-乙烷-丙烷三组分Kij 为对称矩阵 kij [ [0.0, 0.010, 0.018], [0.010, 0.0, 0.006], [0.018, 0.006, 0.0], ] # 演示用的三组分参数非完整物理库 C2H6 {k: v for k, v in CH4.items()} C3H8 {k: v for k, v in CH4.items()} C2H6[A0], C3H8[A0] 5.8e2, 1.1e3 z [0.85, 0.10, 0.05] mix mix_params([CH4, C2H6, C3H8], z, kij) rho_mix solve_density(320.0, 8000.0, mix) Z_mix 8000.0 / (rho_mix * R * 320.0) print(f混合气 rho {rho_mix:.5f} mol/L, Z {Z_mix:.5f})这里要特别说明两点。第一A₀ 的交互项指数是 (1 - kᵢⱼ)而多数资料中 C₀、D₀ 是三次方E₀ 有时也按三次方处理这一点不同文献写法不同混合规则必须和参数表配套不能看哪个顺眼用哪个。第二C2H6 和 C3H8 只是用甲烷参数做形变示意真实乙烷丙烷参数必须查手册或可靠的物性库否则算出来的 Z 只有教学意义不能进正式项目。混合气算例里如果交互系数从 0 改到 0.02密度结果通常只差千分之几但在液相区和临界区附近这个差异会被放大到百分之几甚至百分之十几。所以组分越重、工况越靠近临界kᵢⱼ 的准确性就越重要。4. 从 Z 到工程参数焓偏差与熵偏差怎么算才不虚4.1 偏差函数直接套热力学关系式压缩机功率计算需要焓值节流过程需要熵值但 BWRS 方程给出的是 P、ρ、T 关系不能直接积分得到焓熵。工程上常用偏离函数法先算出真实气体相对理想气体的焓偏差和熵偏差再加上理想气体基准值。推导不展开只看结果。定义密度 ρ、温度 T、压力 P压缩因子 Z P/(ρRT)那么焓偏差和熵偏差可写成ΔH RT(Z - 1) E(ρ, T)ΔS -R ln Z E(ρ, T)/T其中 E(ρ, T) 是对密度从 0 到 ρ 的解析积分结果E (A₀ 3C₀/T² - 4D₀/T³ 5E₀/T⁴)ρ (a 2d/T)ρ²/2 - α(a 2d/T)ρ⁵/5 (3c/(γT²))[(1 γρ²/2)exp(-γρ²) - 1]这套公式成立的前提是参数表单位一致且所有推导基于前文给出的 BWRS 压力表达式。如果把原始 BWR 的压力项混进来E 的系数就对不上焓熵偏差会系统性偏移。所以代码实现时最好把压力计算和偏差计算放在同一个模块里共用同一套参数避免两处来源不一致。4.2 代码落地求解密度的同时输出焓熵偏差把上面公式落到 Python就是在前面的基础上增加一个偏差计算函数def enthalpy_entropy_deviation(rho, T, P, prm): 返回 (焓偏差 J/mol, 熵偏差 J/(mol·K))。 Z P / (rho * R * T) g prm[gama] # 积分中的一次性温度系数 A (prm[A0] 3.0 * prm[C0] / T**2 - 4.0 * prm[D0] / T**3 5.0 * prm[E0] / T**4) B prm[a] 2.0 * prm[d] / T C prm[alpha] * (prm[a] 2.0 * prm[d] / T) g_rho2 g * rho * rho if g_rho2 300: exp_term 0.0 else: exp_term math.exp(-g_rho2) E_integ A * rho 0.5 * B * rho**2 - 0.2 * C * rho**5 E_integ (3.0 * prm[c] / (g * T * T)) * ((1.0 0.5 * g_rho2) * exp_term - 1.0) dH R * T * (Z - 1.0) E_integ dS R * math.log(Z) - E_integ / T return dH, dS T, P 310.0, 6000.0 rho solve_density(T, P, CH4) dH, dS enthalpy_entropy_deviation(rho, T, P, CH4) print(fdH {dH:.2f} J/mol, dS {dS:.4f} J/(mol·K))这里出现了一个工程上容易忽视的单位换算问题R 取 8.314462618 kPa·L/(mol·K)而 1 kPa·L 等于 1 J所以焓偏差的自然单位是 J/mol。如果后续要 kJ/kmol需要在这个函数外自行除以 1000 或乘上分子量不要在函数里写死因为不同模块需要的单位体系不一样。exp_term 的溢出保护很有必要。在液相高密度下γρ² 可能达到几十甚至几百指数项趋于零直接调用 math.exp 在极端参数下可能数值下溢但不会报错可一旦参数表里的 γ 因为混合规则算得偏大指数项在中等密度下就趋零会导致压力曲线出现异常拐点。设一个 g_rho2 大于 300 时直接赋值 0 的阈值是最简单的防护手段。焓熵偏差的精度验证有个笨办法但很有效把压力降到常压附近理想气体偏差应趋近于零把同一点位的结果和手册里的焓熵表对比偏差应小于百分之几。如果对不上九成不是积分公式错而是参数表或单位体系串了。5. 工程落地避坑单位混乱、参数张冠李戴、迭代不收敛5.1 二分法区间不够高压液相直接报错现象计算低温液相或者压力超过 20MPa 的气相时solve_density 抛出 RuntimeError或者返回一个明显偏大的密度。原因默认的 rho_hi 初始值只有 rho_ig 的 4 倍而液相密度可能比理想气体密度大几十倍初始上界根本压不到目标压力虽然循环中有翻倍逻辑但如果 BWRS 压力曲线在大密度区先增后减出现多个极值点翻倍后的右端点可能恰好落在第二个极值下方循环会一直翻到 500 还触发不了条件。解决对液相工况直接把 rho_hi 设为 40 到 60 mol/L 的量级并在循环里同时检查右端点函数值的单调趋势。一个小技巧是先用 PR 方程粗算一个密度作为二分区间参考再套进 BWRS 缩小范围能省很多迭代次数。5.2 单位体系串了R 值错得离谱现象同样一组输入参数别人算出来 Z 是 0.92你算出来 1.35而且压力越高偏差越大。原因这是踩坑率最高的一类参数表来自 kPa、L、mol、K 体系但代码里 R 用了 0.08314 bar·L/(mol·K)或者 A₀ 的值从手册抄来是 bar 单位却直接填进 kPa 体系。解决在代码最上方用一行注释写清楚单位约定并把 R 定义成常量所有压力参数在进入状态方程前统一换算成 kPa。我在做物性模块时会把参数表源文件的单位也读进程序运行时检查一遍一旦发现压力单位与参数单位不一致直接报错而不是让它算出一个看似合理的坏数。5.3 混用两套参数表A₀ 和 C₀ 张冠李戴现象气相 Z 还算正常到了高压或低温液相密集区密度偏差大到难以置信焓熵偏差更是完全不可用。原因许多手册提供了修正型 BWR 方程的参数表与 BWRS 虽然长得像但 D₀、E₀ 的幂次和指数项的系数并不完全相同。把原始 BWR 的 A₀ 拿到 BWRS 的公式里用表面看量纲一致实际会破坏与 C₀、D₀ 的温变平衡。解决参数表必须来自同一套方程的同一来源混用前先核对压力表达式特别是 D₀、E₀ 是除以 T³ 还是 T⁴α 项是否存在指数中是 γρ² 还是 γρ⁶。这个核对过程看起来很机械但值得每次换参数库时做一遍。5.4 混合规则只对了一半重烃组分密度全线跑偏现象甲烷占 90% 的干气算得不错换成含 10% 丙烷丁烷的湿气后液相密度算出来比实验值高 8% 到 12%。原因混合规则里 A₀ 用了几何平均加二元交互系数但 C₀、D₀、a、b、c、d、α 全部偷懒用了线性平均重烃组分贡献被严重低估。BWRS 对重烃的敏感度恰恰集中在这些非线性的混合项上线性平均会让混合参数落在拟合区域的边缘。解决严格按照 Starling 体系使用配套混合规则a、b、c、d、α 用 1/3 次方平均γ 用平方根平均C₀、D₀ 按交互系数三次方修正。交互系数没有实验值可用时宁可用 0 也别乱填填错比不填更危险。6. 验证技巧与参数回归习惯最后一个环节也是我建议你长期保留的工作习惯先验证再信结果。每次写完 BWRS 计算模块我至少做三件事。第一跑一个压力趋于零的算例比如 0.1MPa、300K 的甲烷Z 应该无限接近 1偏差应在 0.1% 以内。第二跑文献或手册中已有的天然气混合物密度点把模拟值和实验值放在同一张表里算相对偏差而不是只看单个点。这样做能快速暴露单位或混合规则的问题。一旦发现偏差系统性随密度增大而增大优先检查参数表和混合规则而不是怀疑迭代代码。第三件事更进阶一些如果你的项目有实验密度数据可以考虑自己回归 BWRS 参数。做法是以实验密度为基准构造目标函数拟合时让 BWRS 压力与实验压力残差最小同时加入饱和液相密度点约束。最小二乘回归需要小心十一个参数之间强相关直接自由拟合容易陷入局部极小得到一组在统计上很好用但物理上荒谬的参数。我的习惯是先用文献值做初值固定 γ 和 α只回归 A₀、B₀、C₀、a、b、c、d最后再放开 γ 微调。这样得到的参数组更容易收敛也更方便与同类流体横向对比。这套验证习惯救过我很多次特别是从别人手里接过计算结果时先用极限工况打点再谈精度。BWRS 不是万能方程但在天然气领域只要参数源一致、混合规则正确、求根区间合理它的稳定性和精度确实对得起那一点实现成本。希望这次梳理能帮你在自己的物性模块里少走几步冤枉路。本文还有配套的精品资源点击获取