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

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

0bwrs數(shù)學(xué)模型在計(jì)算輸氣管道的水力熱量時(shí)，需要計(jì)算天然氣的物理參數(shù)。天然氣密度是一個(gè)重要的參數(shù)。這是解壓縮因子、氣體熵、體積、比熱、泡沫破裂等參數(shù)的基本數(shù)據(jù)。BWRS狀態(tài)方程是一個(gè)非線(xiàn)性多參數(shù)高次方程,如何精確快速的計(jì)算出密度值,對(duì)于工藝計(jì)算具有重要意義。為此,給出了應(yīng)用BWRS方程求解天然氣密度數(shù)值解的兩種快捷方法。1各參數(shù)表1BWRS狀態(tài)方程是一個(gè)多參數(shù)狀態(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——系統(tǒng)壓力,kPa;T——系統(tǒng)溫度,K;ρ——?dú)庀嗷蛞合嗝芏?kmol/m3;Rm——通用氣體常數(shù),8.3143kJ/(kmol·K)。在此,為了在求解各參數(shù)中以國(guó)際單位制為基本單位制,將BWRS方程改寫(xiě)為如下形式:Ρ=ρ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——系統(tǒng)壓力,kPa;T——系統(tǒng)溫度,K;ρ——流體密度,kg/m3;R——?dú)怏w常數(shù),kJ/(kg·K)。上述式中,A0,B0,C0,D0,E0,a,b,c,d,α,γ為狀態(tài)方程的11個(gè)參數(shù)。對(duì)于純組分的這11個(gè)參數(shù)可由臨界參數(shù):臨界溫度Tci,臨界密度ρci及偏心因子wi的下列關(guān)聯(lián)式中求得:ρ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——通用常數(shù)(i=1,2,3,…,11)。對(duì)于混合物,BWRS方程應(yīng)采用如下混合規(guī)則進(jìn)行計(jì)算: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——?dú)庀嗷蛞合嗷旌衔镏械趇組分的摩爾分?jǐn)?shù);Kij——第i、j組分間交互作用系數(shù)(Kij=Kji)。Kij表示和理論混合物所發(fā)生的偏差,Kji越大,說(shuō)明偏差越大,對(duì)于同一種組分,Kij=0。Staring給出了18種常見(jiàn)組分間的Kij數(shù)據(jù)。2參數(shù)單元推理方法中的參數(shù)2.1摩爾體積值計(jì)算單位換算的通用算法是,平均分子量(單位是kg/kmol)乘以用BWRS基本方程算得的摩爾體積值(單位是kmol/m3),計(jì)算結(jié)果精確。該方法需要在計(jì)算后進(jìn)行單位換算。2.2各單位的量綱分析在文獻(xiàn)中的熱力參數(shù)計(jì)算,均沒(méi)有將單位制統(tǒng)一。在此,介紹單位制統(tǒng)一的直接算法,取式(2)中的第一項(xiàng)和第二項(xiàng)作對(duì)比(臨界密度ρci的單位取為kg/m3,通用氣體常數(shù)R的單位取為kJ/(kg·K),溫度T的單位取為K),這兩項(xiàng)的量綱分析結(jié)果應(yīng)該與壓力P相同,于是可以得到采用式(2)時(shí)的B0的單位為m3/kg,將B0代入式(3)第一項(xiàng)可得到A1,B1為無(wú)量綱參數(shù)。同樣分析可得到其它參數(shù)的單位。從以上分析可知,按照該方法得到的密度單位就是kg/m3。3求解密度3.1求解參數(shù)的確定將已知壓力P,溫度T代入式(2)或(1)進(jìn)行計(jì)算,求得氣體密度ρ。為方便求解,將式(2)或(1)改寫(xiě)為如下形式,并用弦截法求解。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)弦截法求解需設(shè)兩個(gè)初值xk,xk-1;在求解式(5)時(shí),可設(shè)為ρ1=0,ρ2=ΡRΤ。迭代到|ρk-ρk-1|≤ε為止,取ε=10-6時(shí),一般迭代次數(shù)在6次左右即能收斂。3.2xk近自然算法仍將方程改寫(xiě)為式(5)的形式。拋物線(xiàn)法的迭代公式為: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的二階差商。弦截法求解需設(shè)三個(gè)近似根xk,xk-1,xk-2,其中自然假定xk為更接近所求的根x*;在求解式(5)時(shí),可設(shè)這三個(gè)根初值分別為ρ0=0?ρ1=ΡRΤ?ρ2=Ρ0.9RΤ。迭代到|ρk+1-ρk|≤ε為止,取ε=10-6時(shí),一般迭代次數(shù)在2次左右即能收斂。弦截法和拋物線(xiàn)法都是超線(xiàn)性收斂的,拋物線(xiàn)法的收斂速度比弦截法更接近于牛頓法,收斂速度快;但是與拋物線(xiàn)法線(xiàn)相比,它需要多設(shè)一個(gè)初值并且需要計(jì)算一階差商和二階差商。4密度優(yōu)化算法已知:0℃,101.325kPa,氮?dú)獾拿芏圈褳?.2507kg/m3,壓縮因子Z為0.9995;求解:該標(biāo)準(zhǔn)狀態(tài)下,氮?dú)獾拿芏圈押蛪嚎s因子Z;按照直接算法求得:ρ=1.250687kg/m3,Z=0.9993。結(jié)果分析:計(jì)算相對(duì)誤差:e*ρ=ρ*-ρρ*=1.2507-1.25071.2507=0.000%密度計(jì)算結(jié)果在精度范圍內(nèi)沒(méi)有誤差,精度較高。算例2求