李-凱斯勒方程的牛頓-二分求解法

張帥， 何應付， 倫增珉， 計秉玉

(中國石化石油勘探開發研究院提高采收率所， 北京 100083)

多組分混合流體的壓強-體積-溫度(P-V-T)性質是求解流動控制方程時不可或缺的重要參數[1-2]，其精確程度直接決定了模擬結果的可靠性[3]。為了準確地描述混合流體的P-V-T性質，Lee等[4]在對比了大量實驗數據的基礎上提出了基于對比態原理的三參數狀態方程，即李-凱斯勒(Lee-Kesler)方程。相較于SRK[5]、MSRK[6]、PR[7]和GCEOS[8]狀態方程而言，李-凱斯勒方程在描述烴類流體P-V-T性質方面具有更高的準確性[9]。因此，李-凱斯勒方程在油氣藏開發[10]、二氧化碳埋存[11]、油氣儲運[12-13]和化工工藝設計[14-15]等領域得到了極為廣泛的應用。在稠油化學復合冷采過程中，李-凱斯勒方程也經常被用于計算稠油及其各主要成分的物性[16-17]。

求解李-凱斯勒方程的算法主要有牛頓法和二分法兩種。牛頓法的優點是收斂速度快，缺點是對初值敏感性高，選取不恰當的初值會使算法收斂至局部最優而非全局最優解。潘孝忠等[18]利用牛頓法對李-凱斯勒方程進行了迭代求解，提出將原方程中的對比比容Vr變形成對比密度Dr的倒數會更有利于進行迭代計算。謝太浩[19]通過引入臨界對比溫度的概念將李-凱斯勒方程的定義域分為氣、液兩個相區，并分別提出了相應的牛頓迭代初值計算公式，大量計算表明選用該公式生成的初值可以取得令人滿意的結果。此外，劉芙蓉等[20]也提出了一套初值生成規則，并將該規則運用于石油液化氣分離模擬過程中，指導了相關裝置的設計。與牛頓迭代法不同，二分法的優點是易于判斷解的存在性，適用于大范圍搜尋潛在最優解，缺點是計算步驟繁瑣。劉天增[21]將二分法與牛頓迭代法進行了對比，得出了二分法計算過程冗長，收斂速度慢的結論。

因此，現通過對牛頓法與二分法進行優勢互補，基于最小自由能原理的牛頓-二分混合算法用來求解李-凱斯勒方程，以期能在算法魯棒性與收斂速度之間取得更好的平衡，為各類模擬軟件中烴類混合流體P-V-T性質的求解提供理論基礎。

1 李-凱斯勒方程的牛頓-二分混合求解算法

1.1 李-凱斯勒方程

李-凱斯勒方程是一組基于對應態原理提出的三參數狀態方程，其核心假設是混合流體的各項熱力學量均可由簡單流體與參考流體的同一熱力學量線性插值得到，即

(1)

式(1)中：X為熱力學量，它既可以是壓縮因子，也可以是熵、焓和吉布斯自由能等參數；X0和Xr分別為簡單流體和參考流體的熱力學量；ω為混合流體的偏心因子，它可寫作各組分的偏心因子按該組分所占比例加權求和的形式，這里ωr的取值為0.397 8。

定義混合流體的無量綱對比壓強和對比溫度分別為

(2)

(3)

式中：Pc為臨界壓強，Pa；Tc為臨界溫度，K。它們可由各組分的臨界參數按下列規則混合產生，即

(4)

(5)

(6)

(7)

式中：n為混合流體的總組分數；xi為第i種組分在混合流體中的摩爾百分數；ωi為該組分的偏心因子；Vci為該組分的特征比容，m3/mol；R為摩爾氣體常數，取值為8.314 J/(mol·K)。對處于對比壓強Pr和對比溫度Tr條件下的簡單流體與參考流體而言，其壓縮因子Z0和Zr應滿足方程：

(8)

式(8)中：Vr為混合流體的對比比容；B、C、D均為關于Tr的函數，其具體形式為

(9)

(10)

(11)

值得注意的是，式(8)是一個一元高次方程，它可能存在多個實根Vri，i=1，2，…。為了從這些實根中挑出最有物理意義的解，需要進一步對流體的焓、熵和吉布斯自由能進行計算

(12)

(13)

(14)

式(14)中：函數G*和E的形式分別為

(15)

(16)

表1 常量取值表Table 1 The value of constants

根據最小自由能原理，當一個體系處于平衡態時，該體系的吉布斯自由能最小。因此，流體在平衡態下的對比比容Vr應當選取dG最小時所對應的實根。該實根即為李-凱斯勒方程最有物理意義的解。

1.2 牛頓-二分混合求解算法

考慮到二分法易于判斷解的存在性和牛頓法收斂速度快的優點，在對兩種算法進行優勢互補的基礎上首次提出了牛頓-二分混合求解算法。該混合算法的基本思路是：首先在李-凱斯勒方程的定義域內進行撒點，將整個定義域分為若干個區間；然后利用二分法的思路搜尋出存在實根且該實根最有可能滿足最小自由能原理的目標區間；在目標區間內，利用牛頓迭代公式求出實根Vr及其所對應的吉布斯自由能dG；最后，反復迭代這一過程直至dG不再下降。此時所對應的Vr即為滿足最小自由能原理的對比比容。算法的流程圖如圖1所示，具體計算步驟如下。

步驟1給定對比壓強Pr和對比溫度Tr，并將式(8)改寫為

(17)

同時，設定Vr的初值為1，并代入式(14)中計算對應的吉布斯自由能dG。值得注意的是，與謝太浩[19]的研究不同，此處Vr的初值是一個與Pr和Tr無關的常數。定義變量start=1。

步驟2定義一維常數數組，即

α=[-8;-6;-4;-3;-2.5;-2;

-1.5;-1;-0.5;-0.01; 0.01; 0.5; 1;

1.5; 2; 2.5; 3; 4; 6; 8]

(18)

并計算出相應的f(exp(αi)Vr)，其中i=1，2，…，20，代表元素在數組中的編號。沿用二分法的思路來判斷解的存在性，若

f(exp(αi)Vr)f(exp(αi+1)Vr)≤0

(19)

則意味著在[exp(αi)Vr，exp(αi+1)Vr]區間內必然存在實根x使得f(x)=0。將i的取值從1循環至19，求出滿足條件的區間[exp(αi1)Vr，exp(αi1+1)Vr]至[exp(αin)Vr，exp(αin+1)Vr]。

圖1 牛頓-二分混合求解算法流程圖Fig.1 The flow chart of hybrid Newton-Bisection algorithm

步驟3針對各區間邊界exp(αi1)Vr，exp(αi1+1)Vr，…，exp(αin)Vr，exp(αin+1)Vr計算出相應的吉布斯自由能dG。找到具備最低邊界自由能min(dG)的目標區間[exp(αi*)Vr，exp(αi*+1)Vr]，為防止重復計算，當變量start的值取0時，i*應等于不為10的整數。在極端情況下，可能存在兩個或多個區間均具備最低的邊界自由能，此時任取其一即可。

步驟4將變量start的值修改為0。在目標區間[exp(αi*)Vr，exp(αi*+1)Vr]內，從區間中點出發，按牛頓迭代公式

(20)

(21)

上述計算步驟不僅適用于簡單流體，也同樣適用于參考流體，兩者之間的區別僅為式(17)中的常數不同，具體取值如表1所示。在分別求得兩種流體的對比比容Vr后，可根據式(8)～式(14)計算其各項熱力學量，并最終代入式(1)中得到混合流體的P-V-T性質。

1.3 對氣、液相區邊界的考慮

在氣、液相區邊界附近，一種可能出現的情況是簡單流體自由能最低的狀態是氣(液)態，而參考流體自由能最低的狀態是液(氣)態。在此情形下，由于式(1)中Xr-X0的絕對值較大，因此混合流體的熱力學量X會高度依賴于它的偏心因子ω，這與實驗觀測數據[4]不符。

為解決這一問題，可在迭代結束之后，引入經驗性限制條件，即

(22)

(23)

表2 氣、液相區邊界附近迭代計算得到的最優解和次優解Table 2 The optimal and sub-optimal solutions obtained by iterative calculations around the boundary between the gas and liquid phase regions

2 結果分析

在求解李-凱斯勒方程時，牛頓-二分混合求解算法具有對初值不敏感和迭代次數相對較少的優點。表3展示了牛頓法、二分法和牛頓-二分混合法計算結果的對比情況。為了提升計算結果的可比性，這里三種算法選用了相同的初值和相同的收斂精度要求。對比三種算法的迭代次數可以看出，牛頓-二分混合法需要的迭代次數要遠小于二分法，整個收斂過程可以加速1～2倍。這主要是因為其在迭代過程中運用了牛頓法的思想，提升了收斂效率。盡管直接使用牛頓法所需的迭代次數更少，如表3所示，但牛頓法容易收斂至局部最優而非全局最優解，使得最終結果與Lee等[4]提供的參考值之間存在較大偏差。特別是對于Pr=0.05、Tr=0.50這種處在氣相區的狀態點而言，當Vr的初值取1時，很難收斂至正確的解。

牛頓法對初值的高度敏感性使它必須搭配一套非常復雜的初值生成規則使用，并且即使這樣也難以徹底避免收斂至局部最優解。圖2展示了在固定Tr或Pr的條件下，簡單流體的壓縮因子Z0隨Pr或Tr的變化曲線。圖2中藍線代表牛頓法的計算結果，參考了謝太浩[19]提出的初值生成規則，紅線代表牛頓-二分混合算法的計算結果，黑色星號是Lee等[4]提供的參考值。從圖2中可以看出，盡管牛頓法在特定狀態點上能夠求得與參考值完全一致的解，但是在其他區域仍可能出現局部不連續的跳變。相比之下，牛頓-二分混合算法可以提供連續、準確的結果，證明該算法具有更好的魯棒性。除此之外，還復算了Lee等[4]提供的6 000條參考數據，除臨界點Pr=1.00、Tr=1.00外全部吻合。

圖2 壓縮因子隨對比溫度的變化曲線Fig.2 The change of compression factor with contrast temperature

圖3 相對吉布斯自由能隨壓縮因子Z0的變化曲線Fig.3 The change of relative Gibbs free energy with increasing compression factor Z0

表3 牛頓法、二分法和牛頓-二分混合法計算結果對比Table 3 The comparison of results of Newton method， bisection method and hybrid Newton-Bisection algorithm

3 結論

(1)通過對牛頓法與二分法進行優勢互補，首次提出了牛頓-二分混合算法用于求解李-凱斯勒方程。該算法的主要思路是：首先參照二分法的思路搜尋出存在實根且該實根最有可能滿足最小自由能原理的目標區間；然后，在目標區間內，利用牛頓迭代公式求出實根Vr及其所對應的吉布斯自由能dG；最后，迭代上述過程直至dG不再下降。此時所對應的Vr即為滿足最小自由能原理的李-凱斯勒方程的解。

(2)牛頓-二分混合算法的主要優勢在于：①對初值不敏感，無需配合復雜的初值生成規則使用；②收斂速度快，相較于二分法而言可以加速1～2倍；③適用范圍廣，可同時用來求解氣相與液相流體的P-V-T性質。該算法提供了一種不依賴初值的一元高次方程求解思路，可推廣至求解BWRS[23]等其他狀態方程。

(3)對臨界點處簡單流體壓縮因子Z0的取值進行了討論?；谧钚∽杂赡茉恚挥挟擹0取0.291 8時系統才會處于平衡態。