用弦截法和引力法求解bwrs狀態(tài)方程中密度根的簡易途徑

用弦截法和引力法求解bwrs狀態(tài)方程中密度根的簡易途徑

0bwrs數學模型在計算輸氣管道的水力熱量時，需要計算天然氣的物理參數。天然氣密度是一個重要的參數。這是解壓縮因子、氣體熵、體積、比熱、泡沫破裂等參數的基本數據。BWRS狀態(tài)方程是一個非線性多參數高次方程,如何精確快速的計算出密度值,對于工藝計算具有重要意義。為此,給出了應用BWRS方程求解天然氣密度數值解的兩種快捷方法。1各參數表1BWRS狀態(tài)方程是一個多參數狀態(tài)方程,其基本形式為:Ρ=ρ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,α,γ為狀態(tài)方程的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方程應采用如下混合規(guī)則進行計算: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表示和理論混合物所發(fā)生的偏差,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密度優(yōu)化算法已知:0℃,101.325kPa,氮氣的密度ρ為1.2507kg/m3,壓縮因子Z為0.9995;求解:該標準狀態(tài)下,氮氣的密度ρ和壓縮因子Z;按照直接算法求得:ρ=1.250687kg/m3,Z=0.9993。結果分析:計算相對誤差:e*ρ=ρ*-ρρ*=1.2507-1.25071.2507=0.000%密度計算結果在精度范圍內沒有誤差,精度較高。算例2求