李-凯斯勒方程的牛顿-二分求解法

张帅， 何应付， 伦增珉， 计秉玉

(中国石化石油勘探开发研究院提高采收率所， 北京 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的取值进行了讨论。基于最小自由能原理，只有当Z0取0.291 8时系统才会处于平衡态。