简介这是面向石油天然气行业技术人员和高校相关专业学生的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³压力单位用kPap ρ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方程前先判断工况是否落在它的甜区里不然算出来再准也是假的。我用一张表给边界做个整理数值是这类状态方程文献里常见的经验范围具体以你手上的参数表来源为准维度适用窗口边界行为温度约240K ~ 330K管输天然气常见范围低于200K后液相根稳定性变差接近临界点附近误差放大压力常压 ~ 40MPa高于40MPa的高压气田工况需要检查是否超出参数表回归范围组分甲烷、乙烷、丙烷、CO2、N2、H2S等常见天然气组分强极性组分水蒸气和C7重烃比例高时不推荐误差预期压缩因子Z典型误差0.2% ~ 0.5%近临界区和富气工况误差可能到1%以上这里说的“不适合”不是不能算而是误差会超过工程可接受范围。比如水含量高的湿气分子间氢键作用会让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, encodingutf-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_funcNone): 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循环里用了双重求和n5时是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, tol1e-10, max_iter100): 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(fNewton 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_expimport 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_funclambda 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方程在这三个方向上的典型行为方便你对照自己的工况判断输入方向影响路径我的一般处理方式温度升高密度降低Z升高温度探头误差超过±1K时压缩因子偏差可能到0.05%计量场景要注意压力升高密度升高Z先降后升高压段对交互系数更敏感标定数据必须覆盖这一区重组分摩尔分数增大密度增大气相根收窄先确认参数表覆盖该组分再确认kij表覆盖该二元对这轮排查的重点不是算得多精细而是确定你的输入数据里哪个量最不可靠。比如站场压力变送器的检定精度是±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错了。原因大部分情况下不是方程本身的问题而是交互系数没设置或单位体系不统一。排查顺序应该是这样的排查步骤操作判断标准1检查参数表单位参数表如果是英制换算错误会导致密度级联偏差2检查kij表CO2、N2相关二元对是否赋了非零交互系数3检查气相初值是否从理想气体密度出发还是误用了液相初值4检查压力函数打印输出逐项数值看哪一项异常大第4步是最快定位问题的手段把bwrs_pressure里的term1到term5分别打印出来正常情况下term2和term3占主导term4和term5只有在液相高密度区才会显著。如果term5在气相工况下大得离谱那就是指数项里的γ参数读错了行。这一条排查顺序我每次上线BWRS计算服务都要过一遍能省下大量排查时间。6. 快速验证与批处理把BWRS算成一组管线数据6.1 用一个实测点验证代码正确性拿到这套BWRS资源后我建议你做的第一件事不是直接上管线数据而是先找一个纯甲烷的实测点验证代码整体正确性。纯甲烷的实验数据最完整也最不容易被组分误差干扰。选一个典型工况比如T300K、p5MPa纯甲烷的压缩因子应该在0.93附近具体值以你参考的数据表为准。跑一下密度反算Z如果偏差超过2%先检查单位体系如果偏差在0.5%内再进混合组分测试。这一步花不了几分钟但能帮你确认参数表、混合规则、求解器三个模块是通的。我自己吃过亏拿到参数表后跳过了单组分验证直接上混合气数据结果计算偏差大来回排查了两天才发现是参数表里有一行数据多了个小数点。从那以后我每次拿到状态方程代码都强制先跑一遍纯甲烷单点验证再进管线数据。希望帮到你。本文还有配套的精品资源点击获取