基于改進二分法的立方型狀態(tài)方程的求解

＊

（蘭州理工大學(xué) 石油化工學(xué)院 甘肅 730050）

1.問題提出

流體的p-V-T關(guān)系是最基本的熱力學(xué)性質(zhì)，是化工熱力學(xué)的重要教學(xué)內(nèi)容之一[1]，也是工程應(yīng)用計算流體性質(zhì)的主要工具。流體的p-V-T性質(zhì)可以直接用于管道或容器的設(shè)計，也是其他不能直接通過實驗測量的熱力學(xué)性質(zhì)（如焓、熵、內(nèi)能、Gibbs自由能、化學(xué)位、逸度系數(shù)等）的數(shù)據(jù)基礎(chǔ)，在化工過程能量分析和相平衡研究中具有重要的地位[2]。立方型狀態(tài)方程作為流體p-V-T關(guān)系的數(shù)學(xué)模型，具有形式簡單、精度高、適用范圍廣等特點，在工程上獲得廣泛的應(yīng)用。立方型狀態(tài)方程是摩爾體積的三次方形式的方程，準確求取體積根是其使用的一個重要環(huán)節(jié)。工程上，普遍的求解方式是直接迭代法或者牛頓迭代法等，這兩種迭代法均存在初值的選取問題。當初始值選擇不合適時迭代不收斂；而程序求解前難以識別出初始值是否合適，需要根據(jù)迭代結(jié)果來判斷；此外手動調(diào)節(jié)初值沒有明確的調(diào)整方向，具有一定盲目性。相對于直接迭代法和牛頓迭代法，二分法（亦稱對分法）可以規(guī)避迭代法可能不收斂的缺點，同時其收斂速度也很快，是一種簡單易行的非線性方程數(shù)值求解方法[3-5]。其借助動態(tài)更新求解區(qū)域的邊界，即雙參數(shù)調(diào)節(jié)法逼近真值實現(xiàn)方程的求解。本文在傳統(tǒng)二分法的基礎(chǔ)之上，提出了改進的二分法，即單參數(shù)調(diào)節(jié)法，將其應(yīng)用于立方型狀態(tài)方程根的求解，以促進其在化工熱力學(xué)的教學(xué)與工程實踐中的應(yīng)用。

2.改進的二分法

傳統(tǒng)二分法的基本原理為：函數(shù)f(x0)在(x1,x2)內(nèi)連續(xù)，且f(x1)·f(x2)＜0，則f(x)在(x1,x2)內(nèi)有解。基于此原理，在方程求解之前，即可確定方程的有解區(qū)間。在有解區(qū)間，計算x1和x2的中值x0(x0=(x1+x2)/2)，并計算f(x0)；如果f(x0)與f(x1)異號，則方程的根必在(x1,x0)區(qū)間，反之，方程的根則在(x0,x2)區(qū)間。在有解區(qū)間內(nèi)，不斷二分，更新有解區(qū)間的上限和下限，可將方程的解確定在更小的區(qū)間，最終確定方程的解。由此，第n次迭代時，有解區(qū)間確定到初始區(qū)間的1/2n；并且有解區(qū)間為(x1,xn-1)或(xn-1,x2)。即不斷更新邊界值x1和x2（兩參數(shù)調(diào)節(jié)），促使其中值滿足方程。

通過對傳統(tǒng)二分法的分析可知，第n次迭代時，有解區(qū)間確定到初始區(qū)間的1/2n，第n次的有解區(qū)間的一個邊界是第n-1次迭代后有解區(qū)間的中點，第n次迭代的中點x0到邊界的距離為1/2n-1，第n-1代有解區(qū)間中點與第n代間距為初始區(qū)間的1/2n+1，因此將傳統(tǒng)二分法的兩參數(shù)調(diào)節(jié)改為有解區(qū)間中點的單參數(shù)調(diào)節(jié)。另外傳統(tǒng)二分法在迭代過程中x1和x2相互對照，并不斷刷新x1和x2。每次迭代都要重新計算f(x1)和f(x2)，并計算并判斷f(x1)·f(x2)的正負號。改進的二分算法以x1或x2作為參照來調(diào)節(jié)x0，且全過程f(x1)和f(x2)均為定值，只需迭代前計算一次。

為簡化參數(shù)調(diào)節(jié)過程，對傳統(tǒng)二分法進行改進為單參數(shù)調(diào)節(jié)法，即在確定有解區(qū)間(x1,x2)之后，保持兩個邊界值不變，反而調(diào)節(jié)x0逐次增加或減少有解區(qū)間的1/2n+1，促使其滿足f(x0)=0。改進二分法的算法流程圖（如圖1所示）。比較發(fā)現(xiàn)，兩種算法迭代n次的精確度均為初始區(qū)間的1/2n+1，但改進的二分法由原來的雙參數(shù)調(diào)節(jié)改進單參數(shù)調(diào)節(jié)，具有步驟更少、計算更快的特點。

圖1 改進的二分法算法流程圖

改進二分法的算法基本步驟：

（1）根據(jù)f(x1)·f(x2)＜0的判據(jù)，確定二分法的初始求解區(qū)間。在立方型狀態(tài)方程根的求解時，初始求解區(qū)間的確定策略如圖2和圖3所示。

在初始求解區(qū)間內(nèi)，采用改進的二分迭代策略：若中間值x0滿足f(x0)=0則輸出精確解x0；用f(x0)·f(x1)＜0判斷兩函數(shù)值是否異號，若f(x0)與f(x1)異號，則x0調(diào)小，第一次調(diào)節(jié)幅度為初始區(qū)間長度的1/4，之后每次調(diào)節(jié)幅度為前一次調(diào)節(jié)幅度的一半，第i次迭代的調(diào)節(jié)幅度為丨x1-x2丨/2（i+1）；同理，若f(x0)與f(x1)同號，則x0調(diào)大，第j次迭代的調(diào)節(jié)幅度丨x1-x2丨/2（j+1）。

（2）當達到規(guī)定的迭代次數(shù)時，輸出結(jié)果。

3.立方型狀態(tài)方程的求解原理

由改進的二分法的原理可知，有解區(qū)間的選擇非常關(guān)鍵。采用二分法求解流體p-V-T方程時，首先選取有解區(qū)間(V1,V2)計算有解區(qū)間的中點V0，將溫度T和中點V0作為已知條件代入立方型狀態(tài)方程中，計算出壓力P0，則存在迭代方程為P(T,V)=P0-P。若P0＜P，則V0比方程的根Vs偏大，否則就是偏小；調(diào)整V0直至迭代方程收斂。

然而，立方型狀態(tài)方程是關(guān)于體積V的三次方程，在T＜TC時，方程有3個實數(shù)解，其中最大解為飽和汽相體積VSV，最小解為飽和液相體積VSL。中間解滿足，即表示“壓脹”，而現(xiàn)實中物質(zhì)只會被壓縮，因此該解無意義。因此，采用二分法求解vdW、RK、SRK和PR方程等立方型狀態(tài)方程時，存在以下問題：

（1）對于兩個有效解的情況，需確定兩個區(qū)間，且區(qū)間內(nèi)僅有一個解，若區(qū)間包含多個解將會出現(xiàn)丟解；

對于已知溫度T和壓力P求解體積V的問題，由于大多數(shù)情況下，存在：參數(shù)b＜最小根Vmin＜臨界體積Vc＜無意義根＜最大根Vmax，針對PR方程如表1所示。因此，主要是避開V=b和無意義根。

表1 幾種物質(zhì)對應(yīng)PR方程的參數(shù)及根

當已知溫度T和壓力P求解飽和汽相體積Vsv時，通過選取理想氣體體積Vig=RT/P，將Vig代入立方形狀態(tài)方程計算出壓力（記為Pig），Pig與條件中的給定壓力P做比較。因為，所以Pig偏小時必有Vig偏大，反之Vig偏小。若Vig偏大則Vig作為上限，否則作為下限，另一個界限用預(yù)估理想體積偏差程度k來確定。

即飽和汽相體積的有解區(qū)域V1，V2為：

大多數(shù)情況下最大根比無意義解大上幾十到上百倍，所以不必考慮無意義解落入初始區(qū)間的情況。氣體的摩爾體積越大(也即高溫低壓)，分子間距越大，氣體與理想氣體偏差程度越小，此時k值應(yīng)較小。多數(shù)情況下，理想氣體與真實氣體偏差程度不會超過20%，因此程序中k值的默認值設(shè)置為0.2即可，初始區(qū)間無解時可調(diào)大k值，絕大多數(shù)情況下k值不會超過1。有解邊界初始值確定算法如圖2所示。

圖2 已知T和P求解飽和汽相體積Vsv時，有解邊界的確定算法

求解飽和汽相體積Vsl時，當條件中所給溫度高于臨界溫度時，該物質(zhì)不能被液化，此時不必求解。若存在飽和液相解時，根據(jù)P-V圖可知，臨界體積為飽和液相體積的最大值，因此臨界體積就是液相體積的上限。由于所以參數(shù)b是體積的下限。V=b處是一個無窮間斷點，該點不能參與運算，可以把下限設(shè)定為V=1.0001b，即表示下限從左端很接近b。即有解區(qū)間為（b，VC）。考慮到無實際意義的中間解可能落進該區(qū)間，造成無法迭代，因此引入?yún)?shù)m，將上限變?yōu)椤Ｈ裟J區(qū)間包含兩個解，則必然不滿足迭代條件，可通過調(diào)大m來排除無意義解，但m過大可能排除掉兩個解，導(dǎo)致迭代往邊界收斂，可通過程序逐步調(diào)大，調(diào)節(jié)程序與k值調(diào)節(jié)類似，調(diào)節(jié)步幅應(yīng)小一些，如圖3所示。由于大多數(shù)情況下，存在：參數(shù)b＜最小根Vmin＜臨界體積Vc＜無意義根＜最大根Vmax，因此m默認值設(shè)為0，即默認初始區(qū)間為（b，VC），m越大上限越靠近參數(shù)b。

圖3 已知T和P求解飽和液相體積Vsl時，有解邊界的確定算法

4.案例實現(xiàn)

以已知溫度和壓力，求異丁烷摩爾體積的案例為例（該案例選自馮新等編著的化工熱力學(xué)第二章例題2-6）[6]。由于改進的二分法是單參照的單值調(diào)節(jié)，容易用Excel實現(xiàn)該程序，借助WPS Excel迭代10次計算數(shù)據(jù)見表2，其中理想氣體偏離度為0.2，飽和液相解的初始區(qū)間為（1.0001b，VC），同時根據(jù)二分法的精確度與迭代次數(shù)的關(guān)系，計算出根值的誤差區(qū)間。計算結(jié)果為VSV=(6.0685±0.0011)L/mol（教材的求解值6.068L/mol，兩種方法的結(jié)果一致），Vsl=(0.10048±0.00019)L/mol。

表2 基于改進二分法PR方程求解異丁烷的摩爾體積

5.結(jié)論

本文在傳統(tǒng)二分法的基礎(chǔ)上，提出了改進的二分法、立方型狀態(tài)方程飽和汽相體積和飽和液相體積有解區(qū)間的確定策略。主要結(jié)論如下。

（1）對于立方型狀態(tài)方程的求解，二分法迭代能夠通過程序一步識別初始求解是否合適，初值選擇方便快捷。

（2）改進的二分法將原來的雙參數(shù)調(diào)節(jié)改進為單數(shù)值參數(shù)調(diào)節(jié)，具有步驟更少、計算更快的特點。

（3）針對立方型狀態(tài)方程飽和汽相體積和飽和液相體積的有解區(qū)域的確定是合理有效的，方便確定氣液相的摩爾體積初值，以便于準確、高效求解立方型狀態(tài)方程。