基于改進二分法的立方型狀態方程的求解

＊

（蘭州理工大學 石油化工學院 甘肅 730050）

1.問題提出

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

2.改進的二分法

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

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

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

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

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

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

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

（2）當達到規定的迭代次數時，輸出結果。

3.立方型狀態方程的求解原理

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

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

（1）對于兩個有效解的情況，需確定兩個區間，且區間內僅有一個解，若區間包含多個解將會出現丟解；

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

表1 幾種物質對應PR方程的參數及根

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

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

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

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

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

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

4.案例實現

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

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

5.結論

本文在傳統二分法的基礎上，提出了改進的二分法、立方型狀態方程飽和汽相體積和飽和液相體積有解區間的確定策略。主要結論如下。

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

（2）改進的二分法將原來的雙參數調節改進為單數值參數調節，具有步驟更少、計算更快的特點。

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