用弦截法和引力法求解bwrs状态方程中密度根的简易途径

用弦截法和引力法求解bwrs状态方程中密度根的简易途径

0bwrs数学模型在计算输气管道的水力热量时，需要计算天然气的物理参数。天然气密度是一个重要的参数。这是解压缩因子、气体熵、体积、比热、泡沫破裂等参数的基本数据。BWRS状态方程是一个非线性多参数高次方程,如何精确快速的计算出密度值,对于工艺计算具有重要意义。为此,给出了应用BWRS方程求解天然气密度数值解的两种快捷方法。1各参数表1BWRS状态方程是一个多参数状态方程,其基本形式为:Ρ=ρRmΤ+(B0RmΤ-A0-C0Τ2+D0Τ3-E0Τ4)ρ2+(bRmΤ-a-dΤ)ρ3+α(a+dΤ)ρ6+cρ3Τ2(1+γρ2)exp(-γρ2)(1)P=ρRmT+(B0RmT−A0−C0T2+D0T3−E0T4)ρ2+(bRmT−a−dT)ρ3+α(a+dT)ρ6+cρ3T2(1+γρ2)exp(−γρ2)(1)式中P——系统压力,kPa;T——系统温度,K;ρ——气相或液相密度,kmol/m3;Rm——通用气体常数,8.3143kJ/(kmol·K)。在此,为了在求解各参数中以国际单位制为基本单位制,将BWRS方程改写为如下形式:Ρ=ρRΤ+(B0RΤ-A0-C0Τ2+D0Τ3-E0Τ4)ρ2+(bRΤ-a-dΤ)ρ3+α(a+dΤ)ρ6+cρ3Τ2(1+γρ2)exp(-γρ2)(2)P=ρRT+(B0RT−A0−C0T2+D0T3−E0T4)ρ2+(bRT−a−dT)ρ3+α(a+dT)ρ6+cρ3T2(1+γρ2)exp(−γρ2)(2)式中P——系统压力,kPa;T——系统温度,K;ρ——流体密度,kg/m3;R——气体常数,kJ/(kg·K)。上述式中,A0,B0,C0,D0,E0,a,b,c,d,α,γ为状态方程的11个参数。对于纯组分的这11个参数可由临界参数:临界温度Tci,临界密度ρci及偏心因子wi的下列关联式中求得:ρciB0i=A1+B1wi;ρ3ciαi=A7+B7wiρciA0iRΤci=A2+B2wi;ρ2ciaiRΤ3ci=A8+B8wiρciC0iRΤ3ci=A3+B3wi;ρ2cidiRΤ2ci=A10+B10wiρ2ciγi=A4+B4wi;ρciD0iRΤ4ci=A9+B9wiρ2cibi=A5+B5wi;ρciE0iRΤ5ci=A11+B11wie-3.8wiρ2ciaiRΤci=A6+B6wi}(3)ρciB0i=A1+B1wi;ρ3ciαi=A7+B7wiρciA0iRTci=A2+B2wi;ρ2ciaiRT3ci=A8+B8wiρciC0iRT3ci=A3+B3wi;ρ2cidiRT2ci=A10+B10wiρ2ciγi=A4+B4wi;ρciD0iRT4ci=A9+B9wiρ2cibi=A5+B5wi;ρciE0iRT5ci=A11+B11wie−3.8wiρ2ciaiRTci=A6+B6wi⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(3)式中Ai,Bi——通用常数(i=1,2,3,…,11)。对于混合物,BWRS方程应采用如下混合规则进行计算:A0=∑i∑jyiyjA0.50iA0.50j(1-Κij);a=[∑iyia1/3i〗3B0=∑iyiB0i;b=[∑iyib1/3i〗3C0=∑i∑jyiyjC0.50iC0.50j(1-Κij)3;c=[∑iyic1/3i〗3D0=∑i∑jyiyjD0.50iD0.50j(1-Κij)4;d=[∑iyid1/3i〗3E0=∑i∑jyiyjE0.50iE0.50j(1-Κij)5,α=[∑iyiα1/3i〗3γ=[∑iyiγ1/2i〗2}(4)A0=∑i∑jyiyjA0.50iA0.50j(1−Kij);a=[∑iyia1/3i〗3B0=∑iyiB0i;b=[∑iyib1/3i〗3C0=∑i∑jyiyjC0.50iC0.50j(1−Kij)3;c=[∑iyic1/3i〗3D0=∑i∑jyiyjD0.50iD0.50j(1−Kij)4;d=[∑iyid1/3i〗3E0=∑i∑jyiyjE0.50iE0.50j(1−Kij)5,α=[∑iyiα1/3i〗3γ=[∑iyiγ1/2i〗2⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(4)式中yi——气相或液相混合物中第i组分的摩尔分数;Kij——第i、j组分间交互作用系数(Kij=Kji)。Kij表示和理论混合物所发生的偏差,Kji越大,说明偏差越大,对于同一种组分,Kij=0。Staring给出了18种常见组分间的Kij数据。2参数单元推理方法中的参数2.1摩尔体积值计算单位换算的通用算法是,平均分子量(单位是kg/kmol)乘以用BWRS基本方程算得的摩尔体积值(单位是kmol/m3),计算结果精确。该方法需要在计算后进行单位换算。2.2各单位的量纲分析在文献中的热力参数计算,均没有将单位制统一。在此,介绍单位制统一的直接算法,取式(2)中的第一项和第二项作对比(临界密度ρci的单位取为kg/m3,通用气体常数R的单位取为kJ/(kg·K),温度T的单位取为K),这两项的量纲分析结果应该与压力P相同,于是可以得到采用式(2)时的B0的单位为m3/kg,将B0代入式(3)第一项可得到A1,B1为无量纲参数。同样分析可得到其它参数的单位。从以上分析可知,按照该方法得到的密度单位就是kg/m3。3求解密度3.1求解参数的确定将已知压力P,温度T代入式(2)或(1)进行计算,求得气体密度ρ。为方便求解,将式(2)或(1)改写为如下形式,并用弦截法求解。f(ρ)=ρRΤ+(B0RΤ-A0-C0Τ2+D0Τ3-E0Τ4)ρ2+(bRΤ-a-dΤ)ρ3+α(a+dΤ)ρ6+cρ3Τ2(1+γρ2)exp(-γρ2)-Ρ=0(5)f(ρ)=ρRT+(B0RT−A0−C0T2+D0T3−E0T4)ρ2+(bRT−a−dT)ρ3+α(a+dT)ρ6+cρ3T2(1+γρ2)exp(−γρ2)−P=0(5)弦截法迭代公式为:xk+1=xk-f(xk)f(xk)-f(xk-1)(xk-xk-1)xk+1=xk−f(xk)f(xk)−f(xk−1)(xk−xk−1),即为:xk+1=f(xk)xk-1-f(xk-1)xkf(xk)-f(xk-1)xk+1=f(xk)xk−1−f(xk−1)xkf(xk)−f(xk−1)(6)弦截法求解需设两个初值xk,xk-1;在求解式(5)时,可设为ρ1=0,ρ2=ΡRΤ。迭代到|ρk-ρk-1|≤ε为止,取ε=10-6时,一般迭代次数在6次左右即能收敛。3.2xk近自然算法仍将方程改写为式(5)的形式。抛物线法的迭代公式为:xk+1=xk-2f(xk)ω±√ω2-4f(x)f[xk,xk-1,xk-2](7)其中ω=f[xk,xk-1]+f[xk,xk-1,xk-2](xk-xk-1)(8)f[xk,xk-1],f[xk,xk-1,xk-2]分别为xk,xk-1的一阶差商和xk,xk-1,xk-2的二阶差商。弦截法求解需设三个近似根xk,xk-1,xk-2,其中自然假定xk为更接近所求的根x*;在求解式(5)时,可设这三个根初值分别为ρ0=0‚ρ1=ΡRΤ‚ρ2=Ρ0.9RΤ。迭代到|ρk+1-ρk|≤ε为止,取ε=10-6时,一般迭代次数在2次左右即能收敛。弦截法和抛物线法都是超线性收敛的,抛物线法的收敛速度比弦截法更接近于牛顿法,收敛速度快;但是与抛物线法线相比,它需要多设一个初值并且需要计算一阶差商和二阶差商。4密度优化算法已知:0℃,101.325kPa,氮气的密度ρ为1.2507kg/m3,压缩因子Z为0.9995;求解:该标准状态下,氮气的密度ρ和压缩因子Z;按照直接算法求得:ρ=1.250687kg/m3,Z=0.9993。结果分析:计算相对误差:e*ρ=ρ*-ρρ*=1.2507-1.25071.2507=0.000%密度计算结果在精度范围内没有误差,精度较高。算例2求