CN101657824B - 使用改进Newton-Raphson算法求解S形非线性函数的稳定方法和设备 - Google Patents

Abstract

提供了一种用于求解非线性S形函数F＝f(S)的设备和方法，该函数表示一个物理系统中的特性S，例如在储藏模拟中的饱和度。Newton迭代(T)在S^V上对函数f(s)执行以确定下一迭代值S^V+1。接着判断S^V+1是否相对S^V位于拐点S^C的另一侧。如果S^V+1相对S^V位于拐点的另一侧，那么S^V+1被设置成S¹，一个改进的新的估计。这个改进的新的估计S¹最好被设置成拐点S^C，或者S^V和S^V+1之间的平均值，即，S¹＝0.5(S^V+S^V+1)。重复上述步骤直到S^V+1在预定的收敛判据内。此外，描述了对于具有重力和毛细管压力的双相和三相流的求解算法。

Description

使用改进Newton-Raphson算法求解S形非线性函数的稳定方法和设备

技术领域

本发明一般地涉及使用Newton-Raphson算法的求解方法，更具体地，涉及用于求解诸如与储藏模拟器中的传输问题相关的S形函数的非线性方程的隐式求解算法。

背景技术

隐式求解算法被广泛应用于解决在诸如储藏模拟的计算机流体动力(CFD)应用中的非线性时间相关的传输问题。这些隐式求解算法通常被用来克服显式方法的时步大小限制。这些传输问题在新的时间层上的解常常使用著名的Newton-Raphson方法来获得。

图1简要地说明了用于求解或至少是近似函数F＝f(s)的解的Newton-Raphson方法的一次简单执行。图1描绘了下列函数的图形：

F＝f(s)    (1)

求解f(s)＝F＝0确定了一个根，即，函数f(s)与水平轴F＝0相交的地方。S的截断位置用字母A表示。如果该位置，或者到其上的令人满意的近似可被确定，那么可以认为方程(1)被求解。

Newton-Raphson方法要求对解做一个初始猜测。该初始猜测在图1中被记为S⁰。在曲线F＝f(s)上的对应点B具有坐标(S⁰，f(S⁰))。在点B处画一条与水平轴相交的切线。该切线的斜率是f′(S⁰)。

代数上，在普通点S⁰上的切线方程是：

F-f(S⁰)＝f′(S⁰)(S-S⁰)    (2)

解S是F＝0时切线的水平坐标。求解该坐标：

S＝S⁰-f(S⁰)/f′(S⁰)    (3)

该坐标随后被用作f(S)＝0的解的第二猜测或近似，并被记为S¹。

因此，

S¹＝S⁰-f(S⁰)/f′(S⁰)    (4)

对方程(1)的根的更好的估计可通过再次获得在坐标(S¹，f(S¹))，即，点C上的切线，并确定它与水平轴的交点来确定。这可是用下式来实现：

S²＝S¹-f(S¹)/f′(S¹)    (5)

可以看到S²相比S¹是对方程(1)的解的更好近似。通过获得产生在S³处的截距的在点D，(S²，f(S²))上的切线可以得到更好的估计。总的来说，如果对F＝f(S)的解的初始近似或猜测是Sⁿ，下一近似Sⁿ⁺¹可用下式确定：

Sⁿ⁺¹＝Sⁿ-f(Sⁿ)/f(Sⁿ)    (6)

方程(6)提供了Newton-Raphson方法的基本公式。通过对函数f(S)重复使用Newton-Raphson方程(6)，对于方程(1)的越来越好的近似可被确定。然而，这假设用来找寻f(S)的解的Newton-Raphson方法的使用是无条件稳定的。事实并不总是如此。

在储藏模拟中，待求解函数通常是传输问题，它在解的最优估计附近被线性化，并在每次迭代中被隐式求解。在一个典型的储藏模拟器中，该传输问题在每次Newton-Raphson迭代中对于储藏模型中的每个单元进行求解。一旦在多次Newton-Raphson迭代之后获得了一个时步的收敛解，该模拟随后继续进行下一时步。如果流量函数如图1所示是凸的(冲击)，凹的(稀疏波)，或线性的(接触间断)，这种方法简单并且无条件稳定，这是大部分计算机流体动力(CFD)问题的情况。然而，对于在大部分储藏模拟中遇到的S形流量函数，如果所选时步过大，Newton-Raphson方法不收敛。在储藏模拟中，这个问题通常用经验性时步控制技术来克服。然而，得到的时步大小常常是保守的或过大的，导致了多余的计算或时步缩减。

因此，需要一种用于求解具有S形流量函数的传输问题的隐式求解算法，它不管时步大小无条件稳定。该发明提供了对于这个需求的解决方案。

发明内容

提供了一种用于求解非线性S形函数F＝f(S)的方法，该函数表示物理系统中的特性S。S形函数具有拐点S^c和f(S)收敛到的预定值F。初始猜测S^v被选作函数F(S)的可能解。Newton迭代(T)在f(S)的S^v处被执行以确定下一迭代值S^v+1。接着判断S^v+1是否相对S^v位于拐点S^c的另一侧。如果S^v+1相对S^v的位于拐点的另一侧，S^v+1被设置成S^l，一个改进的新的估计。然后判断S^v+1是否位于预定的收敛判据内。如果不是，那么S^v被设置成S^v+1，并重复上述步骤直到S^v+1位于预定的收敛判据内。随后，S^v+1的收敛值被选作S形函数F＝f(S)的令人满意的解。改进的新的估计S^l最好被设置成拐点S^c，或者S^v和S^v+1之间的平均值，即，S^l＝0.5(S^v+S^v+1)。可替换地，除了0.5之外的松弛因子可被使用，例如从0.3至0.7中选择的值。

此外，求解算法对于具有重力和毛细管压力的双相和三相流进行描述。

本发明的目标是提供一种使用无条件稳定的改进Newton-Raphson框架求解S形非线性函数以得到基于物理系统的问题的近似解的方法和设备。

附图说明

本发明的这些及其它目标，特点和优点将参考下面的说明书，待审批的权利要求和附图变得更好理解，其中：

图1清晰地描绘了使用传统Newton-Raphson方法的函数F＝f(S)的迭代求解方法；

图2是对于传统储藏模拟器中一个时步的流程图，其中算子T表示使用Newton迭代求解线性化的传输函数；

图3对于包含注入和生产井位置的2D测试情况显示了30×30网格的渗透场k；

图4A，4B和4C显示了在几个初始时步之后，来自图3的渗透场k的饱和度分布；

图5A和5B显示了凸和凹的函数f(S)各自的收敛图；

图6A和6B描绘了S形流量函数，而图6C和6D显示了对应的使用传统Newton-Raphson框架得到的收敛图；

图7A和7B显示了根据本发明使用第一改进Newton-Raphson框架求解如图6A和6B所示的S形函数所获得的收敛图；

图8是对于一个时步的第二改进Newton-Raphson框架的流程图，其中求解线性化的传输函数用算子T表示；

图9A和9B显示了使用改进Newton-Raphson框架的第二实施例求解如图6A和6B所示的S形函数的收敛图；

图10A，10B，10C，10D，10E，10F，10G，10H和10I分别显示了对于三种粘度比μ₁/μ₂和三个时步大小，在0.5个注入的孔隙体积倍数(PVI)之后的饱和度场。

图11显示了使用第一和第二改进Newton框架以及使用具有常数松弛因子的Newton框架在0.5个PVI并且μ₁/μ₂＝1之后的收敛历史；

图12是对于一个时步具有非常数速度场的流程图，其中求解压力方程用算子P表示，而求解线性化的传输方程用算子T表示；

图13显示了对于三种模拟的一个时步的收敛历史，两种使用改进Newton-Raphson框架执行，第三种在所有迭代步骤中使用常数松弛因子；

图14是在本发明的优选实施例中执行的用以求解对于诸如饱和度的物理特性S的S形函数F＝f(S)中的步骤的流程图；

图15清晰地显示了使用根据本发明的改进Newton-Raphson框架的S形函数F＝f(S)的迭代求解过程；

图16A和16B显示了相1的部分流量曲线：(A)m＝1而(B)m＝0.2；

图17A和17B显示了相2的部分流量曲线：(A)m＝1而(B)m＝0.2；

图18A和18B描绘了在单元i和n之间的部分流量：在上侧图中，单元i是逆风的，而在下侧图中，单元n是逆风的，(A)无逆流；(B)逆流在单元i中；(C)逆流在单元n中；以及(D)逆流在单元i和n中；

图19A和19B显示了具有逆流和并流的油相(较轻的流体)的部分流量曲线；

图20A和20B显示了对于(A)相1和(B)相2：m₁＝0.2，m₃＝10并且C_g1＝C_g2＝0的没有重力影响的部分流表面；

图21A和21B显示了对于(A)相1和(B)相2：m₁＝0.2，m₃＝10.C_g1＝5，C_g2＝10的有重力影响的部分流表面；

图22A和22B描绘了对于(A)相1和(B)相2：m₁＝0.2，m₃＝10.C_g1＝0，C_g2＝0的没有重力影响的部分流表面；

图23A和23B描绘了对于(A)相1和(B)相2：m₁＝0.2，m₃＝10.C_g1＝0.5并且C_g2＝1.0的具有小的重力影响的部分流表面；

图24A和24B描绘了对于(A)相1和(B)相2：m₁＝0.2，m₃＝10.C_g1＝5并且C_g2＝10的具有大的重力影响的部分流表面；

具体实施方式

I.储藏模拟的理论背景和数值研究

考虑下面的双曲线方程：

$\frac{\∂S}{\∂t}+\▿\left(\mathrm{fu}\right)=0\mathrm{on\Ω}---\left(7\right)$

方程(7)是表示在双相表层流问题中单相传输的模型方程。因变量S是该相的饱和度(0≤S≤1)，f是S的函数，0≤f≤1，而u是无散度总速度场。为了简单而不失一般性，u被认为是常数，并且假设

f/

S≥0。

公式(7)的非线性传输问题可用如图2的流程图所示的传统Newton-Raphson方法迭代地求解。该过程得到新的时间层上的解Sⁿ⁺¹。S形函数的解的初始猜测S^v被给出。该猜测可以是来自先前时步的饱和度值Sⁿ。在每次迭代v中，线性化传输方程最好对于储藏模型或子模型中的每个单元进行求解(用算子T表示)，得到一个更新的饱和度场S^v+1。一旦从一次迭代到下一迭代的最大绝对饱和度变化没有超出一个适当定义的门限ε，这个非线性循环被认为收敛，并且S^v+1构成了在新的时间层上的解Sⁿ⁺¹。

研究方程(7)的2D离散形式：

这是基于隐式Euler和一阶展开。为了简单起见，方程(8)对于速度向量u＝{u_xu_y}^T的两个分量均大于等于零，并假设正交的均匀间隔的网格的情况进行书写。每个网格单元可用索引对i，j确定，并且在两个坐标方向上的网格间距被定义为Δx和Δy。上标n和v分别表示旧的时间和迭代级别，而Δt是时步大小。为了求解方程(8)，方程(8)中的流量fⁿ⁺¹被线性化近似

最后，Sⁿ⁺¹的近似S^v+1通过求解线性化系统被获得：

$S_{i,j}^{v+1}+c_x(f_{i,j}^{\′v}{S}_{i,j}^{v+1}-f_{i-1,j}^{\′v}{S}_{i-1,j}^{v+1})+c_y(f_{i,j}^{\′v}{S}_{i,j}^{v+1}-f_{i,j-1}^{\′v}{S}_{i,j-1}^{v+1})=S_{i,j}^n-c_x(f_{i,j}^v-f_{i-1,j}^v)$

$-c_y(f_{i,j}^v-f_{i,j-1}^v)+c_x(f_{i,j}^{\′v}{S}_{i,j}^v-f_{i-1}^{\′v}{S}_{i-1,j}^v)+c_y(f_{i,j}^{\′v}{S}_{i,j}^{\′v}-f_{i,j-1}^{\′v}{S}_{i,j-1}^v)---\left(10\right)$

它是通过将方程(9)代入方程(8)得到的。在每次求解方程(10)时强制0≤S^v+1≤1是很重要的。

对于在本说明书中提到的数值研究，30×30的网格单元被用于建模。图3显示了用于通过求解压力p的椭圆方程计算总速度u＝-λΔp的渗透场k和注入和生产井位置：

·(λ

p)＝-qonΩ    (11)

这里，λ是总流度 $\mathrm{\λ}=(\frac{K_{r1}}{{\mathrm{\μ}}_1}+\frac{K_{r2}}{{\mathrm{\μ}}_2}),$ 其中，K_rj是相j的相对渗透率，μ_j是该相的粘度。源项q对生产井是1，对注入井是-1，并且S≡0是饱和度的初始条件。对于简单二次方程的情况，典型的S形部分流量函数可被表示为

$f\left(S\right)=\frac{\frac{S^2}{{\mathrm{\μ}}_1}}{\frac{S^2}{{\mathrm{\μ}}_1}+\frac{{(1-S)}^2}{{\mathrm{\μ}}_2}}---\left(12\right)$

其中，μ₁和μ₂是两相的粘度(注入相用下标1标记)。为了说明上面参考图2解释的隐式求解算法对于过大的时步是不稳定的，给出了三个数值实验的结果，对于这些实验，μ₁/μ₂被分别选取为1，10和0.1。图4显示了在少许初始时步之后的饱和度分布。尽管时步相对注入的孔隙体积倍数(PVI)较小，该结果显示了以饱和度前端的棋盘形图案为特征的强烈的振荡行为。图2的求解算法的非线性循环不收敛，并且迭代次数被限制为200。现在研究这种不稳定性或者振荡行为的起因，并描述可得到对于Newton-Raphson方法的使用无条件稳定的框架的变型。

理解非线性传输问题的稳定性问题的关键在于流量函数的方程(9)的线性化。首先，传输问题f(S)的非线性被分离。然后判断对于这些在0≤S⁰≤1范围内的初始值，迭代框架是否满足收敛条件：

$\underset{v\→\∞}{\lim }f\left(S^v\right)=F,---\left(13\right)$

f(S^v)+f′(S^v)(S^v+1-S^v)＝F.    (14)

在S^v周围的f(S)。Newton迭代公式如下获得：

$S^{v+1}=S^v+\frac{F-f\left(S^v\right)}{f^{\′}\left(S^v\right)}---\left(15\right)$

此外，在每次迭代之后强制0≤S≤1。现在研究条件(13)满足的S⁰和F的组合。如果函数f是S的线性，凸的或凹的函数，条件(13)总是满足。由于线性情况很简单，只给出了凸的和凹的情况的例子，分别用f(S)＝S²和f(S)＝1-(1-S)²表示。

图5A和5B显示了这些凸函数和凹函数的收敛图。收敛图通过选择初始值S⁰和流量函数f的目标值F被创建。然后迭代次数被确定，它要求收敛到距目标值F预定的容许偏差的范围内。水平轴表示起始值S⁰，而垂直轴表示目标值F。该处理在目标值F处对于大量起始值0≤S⁰≤1进行重复。接着该处理对于大量0≤F≤1的目标值进行重复。随后，对于每组值(S，F)达到收敛所需的迭代次数的确定被画出以产生收敛图。

快速收敛被表示为较暗的阴影，而白色显示了参数空间中迭代方程(15)不收敛的区域。在图5A和5B中不存在白色区域，因此该迭代框架对于凸函数和凹函数f(S)是无条件稳定的。

然而，如果流量函数f是S形函数，情况就不同了。S型函数的特征是具有在拐点S^C处相连的凸部分和凹部分，如图6A和56B所示。拐点S^C的特征是函数f(S)的二阶导数f’(S)在其周围改变符号的点。

图6A和6B中，上侧图分别显示了对于μ₁/μ₂＝1和μ₁/μ₂＝0.1由方程(12)定义的f。对应的收敛图被显示在图6C和6D中。大的白色区域的存在清楚地显示了应用于方程(12)的流量函数f的迭代框架对于这些μ₁/μ₂值不是无条件稳定的。如果S⁰的初始猜测很接近最终解，即，如果|f(S⁰-F)|不太大，则能够达到收敛。如果求解传输问题的话，这与限制时步大小或最大饱和度变化是一致的。

图6C和6D中的收敛图显示了如果S⁰＝S^c，迭代方程(15)对于所有可能目标值F均收敛，其中S^c是f的拐点。这个发现表明如果f(S^v)和f(S^v+1)的二阶导数具有相反的符号，设置S^v+1＝S^c。如果这在Newton迭代的最后被完成，改进的Newton-Raphson框架被产生，它对于S形流量函数是无条件稳定的。

图7A和7B显示了对于图6A和6B的S形流量函数的收敛图。不管怎样，如果(S^v)和(S^v+1)的二阶导数具有相反的符号，设置S^v+1＝S^v，这就提供了第一改进Newton-Raphson框架。这可通过计算下式进行校验：

f″(S^v)f″(S^v+1)＜0    (16)

如果二阶导数的乘积是负的，那么S^v和S^v+1在拐点S^c的两侧。

水平轴表示起始值S⁰，而垂直轴是目标值F。暗色区域表示快速收敛，而不收敛对于白色区域中的参数对被获得。注意，在图7A和7B中不存在白色区域，因此第一改进Newton-Raphson框架是无条件稳定的。将该改进应用于传输问题的求解算法上是直观的，并且允许任意大的时步。

确定拐点S^c的精确位置是很难的。如果求解部分流量函数依赖于多个变量的传输方程的系统，这就更加真实。在某些迭代中，查找表可被用来表示部分流量函数f(S)，并且拐点不能被简单地解析地计算。此外，找到复部分流量函数的拐点S⁰是计算复杂的。

因此，可选择地，并且最好是提出在现实设置中更容易实现的第二改进框架。第二改进框架在每次Newton迭代结尾时S^v+1的改进方式上与第一框架不同，其中S^v和S^v+1位于曲线f的拐点S^c的两侧。代替将S^v+1设置为S^c，可使用一种有选择性的低松弛法。例如，如果二阶导数，即f″(S)＝

²f/

S²的符号在新旧迭代级别上不相同，S^v+1可用0.5(S^v+S^v+1)代替。其它级别的松弛也可被使用。也就是说，低松弛可由下式推广：

S^l＝αS^v+(1-α)S^v+1    (17)

其中，S^l＝S的新的估计；而α＝低松弛系数。

第二增强框架的流程图在图8中描绘。由于确定f的二阶导数相对容易，与当S^v和S^v+1在S^c的两侧时将S^v+1设置为S^c的第一改进框架相比，第二可替换框架可以在现有的储藏模拟器中被容易地实现。

与图6C，6D和图7A，7B中显示的收敛图相似，使用可选低松弛的第二可替换改进对应的收敛图被显示在图9A和9B中。水平轴表示起始值S⁰，而垂直轴表示目标值F。暗色区域表示快速收敛，而不收敛对于白色区域中的参数对被获得。与图7A和7B中的收敛图一样，没有白色区域显示了第二可替换改进得到了一个无条件稳定的框架。

如果被应用于非线性传输方程的求解算法，第一和第二改进均用作预测。第一和第二改进框架几乎是等效的。此外，如果时步大小在未改进或传统Newton-Raphson系统的稳定范围内，由本发明的改进引起的效率上的损害没有在数值研究中被发现。另一方面，如果通过在每次Newton迭代结尾的任何地方均使用低松弛来获得稳定性，而不是如上所述在改进框架中可选地使用，那么效率会被严重损害。对于使用由方程(12)定义的S形部分流量函数的上述测试情况，数值结果被显示在图10A-I中。所有模拟用第一和第二改进框架来执行，并且证实任意大的时步可被使用。图10A-I分别显示了在时间＝0.5PVI之后当μ₁/μ₂＝1(顶行——图10A，10B和10C)，μ₁/μ₂＝10(中间行——图10D，10E和10F)，以及μ₁/μ₂＝0.1(底行——图10G，10H和10I)时获得的饱和度场。在这些图左侧列中的结果使用100倍小的时步，或Δt＝0.005PVI来计算，对此未改进的Newton-Raphson框架证明也是稳定的。注意这些时步是用来获得如图4所示的不收敛结果的大小的一半。图的中间和右侧列分别显示了Δt＝0.1PVI和Δt＝0.5PVI时计算的结果。

图11显示了使用上述第一和第二改进Newton-Raphson框架的结果，其中只有当二阶导数f″(S^v)和f″(S^v+1)的符号在迭代步骤中改变时，S^v+1被分别可选地设置为S^c和0.5(S^v+S^v+1)。进一步，研究在其中低松弛被应用于每次迭代的结尾的第三框架，即单纯低松弛框架。垂直轴显示了在迭代步长中的最大饱和度变化，而垂直轴显示了迭代的次数。图11对于时步大小是0.5PVI，并且粘度比率μ₁/μ₂是一的情况进行绘制。对于执行的所有模拟，第一和第二改进框架均被测试，并且在所需迭代方面没能观察到明显的差异。虽然前两个框架的收敛率基本相同，使用简单低松弛的第三框架被显示更为低效。

如果流度因子λ是饱和度的函数，改进Newton-Raphson框架也可被用作求解算法的积木。这对于石油行业特别重要。对于二次相对渗透的情况，流度因子λ被定义为：

$\mathrm{\λ}=k(\frac{S^2}{\mathrm{\μ}1}+\frac{{(1-S)}^2}{\mathrm{\μ}2})---\left(18\right)$

其中k是渗透率；

S＝第一相的饱和度

μ₁＝第一相的粘度；而

μ₂＝第二相的粘度。

作为λ是饱和度的函数的结果，速度场u不再是常数，并且如果λ变化，需要通过求解压力方程(18)进行更新。该算法通过图12中的流程图进行概述。注意图8的流程图再次出现在更新整个速度场u所需循环中的虚线框内。这个收敛解与通过完全隐式的框架获得的解是一样的，其中压力和传输方程被同时求解。

图13显示了对于两种模拟的一个时步的收敛历史，一个由使用可选低松弛的改进Newton-Raphson求解算法执行，另一个由使用简单或单纯低松弛的方法执行。测试情况与前面的研究相同，时步大小是0.5PVI，而粘度比μ₁/μ₂＝1。使用改进算法，求解线性化传输方程需要大约60次，而椭圆压力方程只需要九次(收敛图中的峰的个数)。简单或单纯低松弛在另一方面被证明是更低效的。

用于求解对于物理特性的S形函数的典型步骤

本发明包含用于稳定求解非线性S形曲线F＝f(S)，以确定物理系统的特性S的表示的方法和设备。图14显示了可被用来依照上面章节描述的原理实现本发明的步骤的优选示范性流程图。图15清晰地描述了改进Newton框架。

在改进Newton-Raphson算法中，在S形函数的二阶导数f″(S^v)和f″(S^v+1)的符号在相继迭代之间改变了的情况下，下一迭代值S^v+1被设置成新的估计或稳定性限制S¹。改进的新的估计，或稳定性限制S¹被理想地选择以使得改进Newton-Raphson方法保证稳定收敛。

第一步骤110获得表示物理系统中的特性S的S形函数f(S)。例如，特别感兴趣的物理系统是地下储藏中的流体流量。S形函数f(S)可用数学函数表示，例如上面的方程(12)，或使用查找表的方式获取S形函数的要素。像储藏模拟领域的技术人员所知的那样，诸如部分流量函数的S形函数可根据相对渗透率度量从实验室测试中得到。

求解对于物理系统的非线性问题领域的技术人员将会认识到使用改进Newton-Raphson框架的本发明可被用来求解表示除储藏模拟之外的物理系统的特性的特征的S形函数。

S形函数f(S)通过包含至少一个在拐点S^c处相连的凸部分和凹部分来表征。函数f(S)被求解以确定满足条件F＝f(S)的物理特性S，其中F是目标值。例如，在储藏模拟中，当求解在线性化传输问题的网格单元中的饱和度值S的S形部分流量函数f(S)时，来自先前时间层的流量值F被理想地用作目标值F。可替换地，目标值F也可以通过数值模拟中的非线性迭代来估计。

初始值或猜测S^v在步骤120中被选作满足条件F＝f(S)的物理特性S的可能解。S^v的初始值最好用初始条件或来自先前时步的结果进行估计。尽管次选的，S^v+1的初始值也可用内插方法或随机化方法得到。

在图8中用算子(T)表示的Newton迭代在步骤130中使用当前迭代的值S^v在函数f(S)上被执行以确定下一迭代值S^v+1。该Newton迭代对应于将方程(15)应用于S形函数f(S)以确定下一迭代值S^v+1。在图15中，初始猜测是S⁰。

可选的检测可在步骤140中进行以查看S^v+1是否在可接受的范围内。例如，如果在诸如与方程(12)中类似的函数的S形函数的分母接近零，对于S^v+1的估计值可能相当大。如果S^v+1超过某一边界限制，例如饱和度S小于零或大于一，那么可在步骤150中将S^v+1设置为诸如0.0或1.0的值以保持物理特性S在现实范围内。在图15的例子中，可以看到点(S⁰，f(S⁰))的正切与表示流量值F的水平线的截距没有出现在饱和度上限S^b＝1.0内。因此，S^v+1在图15中第一迭代的情况下被设置为上限S^b＝1.0。

然后在步骤160中判断S^v和S^v+1是否位于S形函数f(S)的拐点S^c的两侧。进行这种判断最好的方式是查看f(S)二阶导数f″(S)在S^v和S^v+1处的符号是否彼此相反。这可用公式(16)计算f″(S^v)和f″(S^v+1)的乘积来检验。

如果二阶导数的乘积是负的，那么S^v和S^v+1位于拐点的两侧。尽管是次优的，其它进行这种判断的方式是明确计算拐点S^C的位置，然后查看S^v和S^v+1是否在的S^c的同一侧。拐点S^c可通过求解f″(S)＝0来计算。这在f(S)足够简单时可解析地实现，而在f(S)比较复杂时在数值上实现。在使用查找表定义S形函数f(S)的情况下，f(S)，f′(S)和f″(S)的初始值表可被产生。然后查找表可被检查以确定拐点S^C的位置。在图15中，S⁰和S^b的值位于拐点S^C的两侧。

如果S^v和S^v+1位于拐点S^c的两侧，那么为了保证在求解S时无条件稳定，在步骤170中，S^v+1被设置成新的估计或稳定性限制S^l。该改进估计S^l被选择以保证Newton-Raphson迭代的收敛。在第一优选实施例中，这个改进估计S^l被设置成等于拐点S^c。在如图15所示的第二或更好的方式中，低松弛方法被用来计算S^l。理想地，S^l根据下面的方程进行计算：

S^l＝αS^v+(1-α)S(^v+1)    (19)

其中0＜α＜1。

优选地，α的范围是0.3＜α＜0.7，更优选地，0.4＜α＜0.6，而最优选地，α＝0.5。使用α＝0.5，S^v+1或S¹被设置成0.5(S^V和S^V+1)。

如果判断S^v和S^v+1位于拐点S^c的同一侧，那么S^v+1最好保持不变而不是被设置成改进估计S^l，即S^c或S^l＝0.5(S^v+S^v+1)。

在步骤180中，判断当前迭代值S^v+1是否达到预定的容许条件。例如，可以要求在S^v和S^v+1的连续迭代值之间的差足够小。

数学上，这可被表示为：

ε＞|S^v+1-S^v|    (20)

其中，ε是预定的容许极限。

可替换地，为了考虑到解S^v+1可以足够接近真实解，也可以要求：

ε＞|F-f(S^v+1)|    (21)

如果S^v+1不被认为是一个令人满意的解，那么在步骤190中S^v被设置成S^v+1，并且继续步骤130-180直到S^v+1的值收敛到预定的容许极限内。图15举例说明了对于饱和度S¹，S²和S³等的连续估计如何收敛到饱和度S的真实值。如果达到充分的收敛，物理值S被认为是令人满意的解。在储藏模拟器的情况下，当对于在其上执行改进Newton迭代的每个网格单元达到收敛之后，对于下一时间层Sⁿ⁺¹在网格单元中的饱和度S随后被设置成S^v+1。

本发明还包括设备或机器可读的确实包含指令程序的程序存储设备。这些指令被机器执行以实现稳定求解表示物理系统中的特性S的非线性S形函数F＝f(S)的方法步骤。该S形函数具有拐点S^c以及f(S)收敛到的预定值F。所述方法包括下列步骤。选择F＝f(S)的可能解的初始值。在函数f(S)上在S^v处执行Newton迭代以确定下一迭代值S^v+1。然后判断S^v+1是否相对S^v位于拐点S^c的另一侧。如果S^v和S^v+1位于相对侧，那么S^v+1被设置成等于改进的新估计S^l。接着进行检查以通过判断S^v+1是否在预定收敛判据内来判断S^v+1是否是特性S的令人满意的解。如果S^v+1没有令人满意地收敛，那么S^v被设置成等于S^v+1。重复上述步骤直到S^v+1在预定的收敛判据内。当达到收敛时，S^v+1被选作问题F＝f(S)的令人满意的解。

II.对于具有重力和毛细管压力的双相流的算法

重力和毛细管压力在Newton-Raphson方法的数值稳定性上的影响将对于双相流的非线性传输方法进行测试。重力和/或毛细管压力在储藏模拟中的影响可导致储藏单元中的逆流。这种逆流在传统的储藏模拟器中会造成严重的数值不稳定。通过扩展上面在第I章中描述的方法，将得到对于具有重力和毛细管压力的双向流的非线性传输方程无条件稳定的算法。

首先，将在部分流量曲线上研究重力的影响。然后用特征参数描述部分流量曲线的特征。接着描述对于具有重力的双向流的无条件稳定算法。改进也可被提供给无条件稳定算法以包含毛细管压力的影响。

重力影响

对于具有重力的双向流的Darcy判据可被写作

${\mathrm{\μ}}_1=\frac{k{k}_{r1}}{{\mathrm{\μ}}_1}(\▿p-{\mathrm{\ρ}}_{1}g\▿D)---\left(22\right)$

${\mathrm{\μ}}_2=\frac{k{k}_{r2}}{{\mathrm{\μ}}_2}(\▿p-{\mathrm{\ρ}}_{2}g\▿D)---\left(23\right)$

这里，k是渗透率，u是相速度，k_r是相对渗透率，μ是粘度，D表示深度而ρ是密度。下标表示相特性。

从方程(22)和(23)，可以容易地推导出：

u₁＝f₁(S)u-f₃(S)C_g    (24)

u₂＝f₂(S)u+f₃(S)C_g    (25)

其中S表示相1(即，水)的饱和度。总速度被定义为

u＝u₁+u₂    (26)

$f_1=\frac{k_{r1}}{k_{r1}+m{k}_{r2}}---\left(27\right)$

$f_2=\frac{m{k}_{r2}}{k_{r1}+m{k}_{r2}}---\left(28\right)$

$f_3=\frac{m{k}_{r1}{k}_{r2}}{k_{r1}+m{k}_{r2}}---\left(29\right)$

这里，m是粘度比率μ₁/μ₂。

重力引起的速度可用变量C_g表示，

$C_g=\frac{\mathrm{kg}({\mathrm{\ρ}}_1-{\mathrm{\ρ}}_2)}{{\mathrm{\μ}}_1{\mathrm{\μ}}_c}\▿D---\left(30\right)$

其中，u_c是特征速度而D是深度。注意，方程(24)和(25)中的速度u，u₁和u₂关于u_c也是无量纲的。

从方程(24)和(25)分析部分流量曲线是有意义的。在图16和17中，部分流量分别对于相1和相2进行绘制，

$F_1=\frac{u_1}{u}---\left(31\right)$

$F_2=\frac{u_2}{u},---\left(32\right)$

具有二次相对渗透率，

$k_{r1}=S_1^2---\left(33\right)$

$k_{r2}=S_2^2={(1-S_1)}^2---\left(34\right)$

如果0＜mC_g≤1，部分流量曲线是单调的，反之，它对于mC_g＞1或mC_g＜0可能变成多值函数。注意部分流量可以小于零或大于一，这取决于重力和速度的方向。如果Newton非线性迭代不能考虑这些曲线的模式，这会导致数值上的困难。注意来自负C_g和正u的部分流量曲线等价于来自正C_g和负u的情况。在实际的具有固定坐标的模拟中，C_g可保持恒正，但是平均速度可正可负。注意不失一般性，具有负平均速度的情况等价于通过选择合适的坐标使平均速度为正的具有负C_g正平均速度的情况。

如果油(较轻的流体)方程被选作初级方程，从方程(25)，油相的部分流变成零S₀ ^*或一S₁ ^*的状态：

$S_0^{*}=1-\sqrt{\left|u\right|/C_g}\mathrm{foru}<0,$ 以及    (35)

$S_1^{*}=1-\sqrt{\left|u\right|/m{C}_g}\mathrm{foru}>0.---\left(36\right)$

如果对于u＞0，mC_g/|u|≤1，部分曲线变成单调曲线；反之如果u＜0，或对于u＞₀，mC_g/|u|＞1，部分流量曲线包括部分曲线变成多值(非单调)的区域。当部分流量曲线具有最小值或最大值点时，那么这里存在在非线性迭代中饱和度并不穿过的拐点(参见图19)。如果相对渗透率类似方程(33)和(34)中的二次方，该拐点可被解析地计算。

$0=m{S}_r^{*3}-{(1-S_r^{*})}^3+\frac{u}{C_g}.---\left(37\right)$

三次方程可被解析求解，并且总能获得至少一个实数根：

$S_r^{*}=\mathrm{\&alpha;}+{(-\frac{q}{2}+\frac{\sqrt{q^2+4p^3}}{2})}^{1/3}+{(-\frac{q}{2}-\frac{\sqrt{q^2+4p^3}}{2})}^{1/3};---\left(38\right)$

其中

$\mathrm{\&alpha;}=\frac{1}{1+m};---\left(39\right)$

$p=\frac{1}{3}({3\mathrm{\&alpha;}}^2+2\mathrm{b\&alpha;}+c);---\left(40\right)$

q＝(α³+bα²+cα+d)；    (41)

$b=-\frac{3}{1+m};---\left(42\right)$

$c=\frac{3}{1+m};---\left(43\right)$

$d=-\frac{1-u/C_g}{1+m}.---\left(44\right)$

饱和度可被分割成两个范围：(1)S＞S_r ^*和S＜S_r ^*。总的来说，饱和度不穿过这个反射点，但如果该单元在无流边界附近，或者挨着具有其它范围的饱和度的单元，该状态应该被允许穿越。由于存在允许流体在没有逆流的情况下进行累加的边界，部分流量曲线并不适用。

有限差分近似

图16和17中的部分流量曲线表明双相流可以是并流和逆流，取决于特征变量C_g，粘度比率m和饱和度S。如图18所示，当单元i和邻近单元n中的部分流没有引起逆流时，普通的有限展开差分近似可被明确地使用。然而，如果逆流在单元i或单元n中发生，如图18所示，上行流是由相决定的。

作为结果，从在压力解决方案中计算的总速度分离的相带有在单元i或n之间的并流或逆流的详细信息。在图18中，描述了能够在单元i和它的邻近单元n之间发生的所有可能情况。假设来自压力解决方案中的总速度u_f是在单元i和n的分界面上的速度，并且在单元i和n内的总速度是一致的，称之为u。

在单元i中，并向和对向的相速度分别用u_p ⁱ和u_r ⁱ表示其总速度u。这些相速度可从方程(24)和(25)中获得。该方程可被表示为

$u_p^i=f(u,S^i,C_g,m),---\left(45\right)$

$u_r^i=u-u_p^i.---\left(46\right)$

类似地，单元n中的相速度可通过下式获得

$u_p^n=f(u,S^n,C_g,m),---\left(47\right)$

$u_r^n=u-u_p^n.---\left(48\right)$

如果单元i是上行流，

$u_f=u_p^i+u_r^n---\left(49\right)$

如果单元n是上行流，

$u_f=u_p^n+u_r^i.---\left(50\right)$

注意假设某一相速度可以是并向或对向的，这依赖于流和重力方向。此外，如果0＜mC_g≤1，逆流并不出现。

在顺序隐式算法中，饱和度方程使用固定的总速度来迭代求解，以保持质量守恒(即，发散自由度)。该算法可用总速度和部分流量公式直接实现。然而，对于逆流的情况，依赖于单个饱和度的原始部分流必须被改进以保持质量守恒。如果相速度u_i从上行流α_i单元中确定，该单元α_i的部分流可被表示为：

$u_i^{{\mathrm{\&alpha;}}_i}=k\frac{k_{\mathrm{ri}}^{\mathrm{\&alpha;i}}}{{\mathrm{\&mu;}}_i^{\mathrm{\&alpha;i}}}(\mathrm{\&Delta;p}-{\mathrm{\&rho;}}_{i}g\&dtri;D)---\left(51\right)$

那么总速度可被写为：

$u=\underset{l}{\mathrm{\&Sigma;}}{u}_l^{\mathrm{\&alpha;l}}=\mathrm{k\&Delta;p}\underset{l}{\mathrm{\&Sigma;}}\frac{k_{\mathrm{rl}}^{{\mathrm{\&alpha;}}_l}}{{\mathrm{\&mu;}}_l^{{\mathrm{\&alpha;}}_l}}-\mathrm{kg}\&dtri;D\underset{l}{\mathrm{\&Sigma;}}\frac{k_{\mathrm{rl}}^{\mathrm{\&alpha;l}}{\mathrm{\&rho;}}_{l}g}{{\mathrm{\&mu;}}_l^{\mathrm{\&alpha;l}}}(\mathrm{\&Delta;p}-\&dtri;D)---\left(52\right)$

$u_l^{{\mathrm{\&alpha;}}_i}=\frac{{\mathrm{\&lambda;}}_i^{{\mathrm{\&alpha;}}_i}}{{\mathrm{\&lambda;}}_T}u+\mathrm{kg}\&dtri;D\frac{{\mathrm{\&lambda;}}_i^{{\mathrm{\&alpha;}}_i}}{{\text{\&lambda;}}_T}\underset{l}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}^{{\mathrm{\&alpha;}}_l}({\mathrm{\&rho;}}_l-{\mathrm{\&rho;}}_i)---\left(53\right)$

${\mathrm{\&lambda;}}_l^{{\mathrm{\&alpha;}}_l}=\frac{k_{\mathrm{rl}}^{{\mathrm{\&alpha;}}_l}}{{\mathrm{\&mu;}}_l^{{\mathrm{\&alpha;}}_l}}---\left(54\right)$

${\mathrm{\&lambda;}}_T=\underset{l}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_l^{{\mathrm{\&alpha;}}_l}.---\left(55\right)$

稳定性检查

非线性传输问题的稳定性问题在于流量函数的线性化。注意饱和度方程现在是该单元的饱和度和所有上行流单元的饱和度的函数。数值框架关于该框架能否对于所有单元从初始饱和度估计0≤S⁰≤1收敛到满足如下收敛判据的最终饱和度S^v进行检查：

$\underset{v\&RightArrow;\&infin;}{\lim }F\left(S^v\right)=F_0,---\left(56\right)$

其中F₀是右矢量。在方程(56)中，饱和度用矢量符号表示。随着F(S)在S^v附近线性化，

$F\left(S^v\right)+{\left(\frac{\&PartialD;F}{\&PartialD;S}\right)}_v\&CenterDot;(S^{v+1}-S^v)=F_0.---\left(57\right)$

可以获得迭代方程

$S^{v+1}=S^v+{(F_0-F)}_v\&CenterDot;{\left(\frac{\&PartialD;F}{\&PartialD;S}\right)}_v^{-1}.---\left(58\right)$

此外，在每次迭代之后，对于所有j强制 $0\&le;S_j^{v+1}\&le;1.$

在有限差分模拟的三维模型中，重力参数C_g依赖于流的方向(即，对于水平流F_h，C_g＝0，而对于垂直流F_v，C_g≠0)。作为结果，必须对于水平流F_h和垂直流F_v考虑两条部分流量曲线。对于mC_g＞1，垂直流F_v可能变成非单调曲线。较轻的流体饱和度(S₀)的方程被求解，而另一相由饱和度限制(S_w＝1-S₀)进行计算。

重力的考量

重力的其它考量的影响有可能改变部分流量曲线F＝f(S)，方程(31)和(32)，以使得这些曲线可能包含反射点S_r ^*。如上所述，C_g，方程(30)被计算，并在最终估计该部分流量曲线，F(S)时在方程(24)和(25)中使用。

在没有逆流的情况中，部分流量曲线将仍像第I章所述没有任何反射点S_r ^*。然而，如果存在逆流，那么诸如图19A-B所示的曲线可能包含一个反射点S_r ^*。

区域I和II可被定义为由反射点S_r ^*分开的区域，区域I是在拐点之下0＜S＜S_r ^*的区域，而区域II是在S_r ^*＜S＜1.0之上的区域。

在正在估计S^v的网格单元靠近无流边界的情况下，上述第I章的方法被使用，而不考虑任何反射点S_r ^*。这是因为逆流不可能在这种条件下发生。可替换地，如果在相邻单元i和n中的S^v位于不同区域I和II，那么上述第1章的方法被使用以允许饱和度S^v+1穿过反射点S_r ^*。

网格单元i是(1)不靠近无流边界，或者(2)单元i和n位于相同区域的情况下，那么S^v+1对于S^v所在的特定区域被确定。例如，如果S^v在区域I中，那么S^v+1将继续位于区域I中。如果f″(S^v)f″(S^v+1)≥0，那么无需低松弛来确保稳定性。尽管如此，在S^v+1位于区域I外侧的情况下，S^v+1将被设置成0或者S_r ^*。如果f″(S^v)f″(S^v+1)＜0，即，在该区域中拐点S^c的两侧，那么低松弛将被使用以确保稳定性。再次，如果S^v+1位于S^v所在的区域I的外侧，在计算S^v+1的更新或低松弛值，即0.5(S^v+S^v+1)之前，S^v+1将被设置成为下限或上限，0或S_r ^*中的一个。

类似地，在区域II，S^v和S^v+1被检查以查看S^v和S^v+1是否位于区域II中拐点S^c的同一侧。如果是的话，无需松弛来保持稳性定——只需要进行检查以确保S^v+1在区域II的范围内。如果S^v和S^v+1位于该区域中拐点S^c的两侧，即f″(S^v)f″(S^v+1)＜0，那么低松弛被使用，即S^v+1＝0.5(S^v+S^v+1)。再次，在进行后续的低松弛计算之前S^v+1必须落在区域II内。

毛细管压力可以很容易的并入到上述步骤中。被简单地设置成C_g-C_a，就像下面在方程(65)中描述的那样。

在该点上，已对于这个特定的网格单元i确定了本次迭代v+1的S^v+1。如第I章所述，收敛检查在模型或子模型的每个网格单元上进行以查看收敛是否在所有网格单元中均已达到以使得一个时步可被认为已完成，并且在该时步无需进一步的Newton迭代。

在网格单元中更新S^v-＞S^v+1的算法

下面的算法将用于具有重力的双相流。

1.检查相速度和重力的方向。

2.如果相速度是并向的，使用常规的稳定性分析；

·  计算定向的部分流量及其二阶导数(F_h ^*和F_v ^*)；

·  计算对于垂直流的C_g(对于水平流，C_g＝0)；并且

·  使用饱和度限制0≤S^v+1≤1。

·  如果F_h″(S^v+1)F_h″(S^v)＜0或F_v″(S^v+1)F_v″(S^v)＜0，对于饱和度迭代使用低松弛，即：S^v+1＝0.5(S^v+S^v+1)。

如果相速度是对向的；

·计算反射饱和度S_r ^*；并且

·使用饱和度限制：如果该单元不靠近边界或者相邻单元并不属于不同的饱和度区域，那么0≤S^v+1≤S_r ^*或S_r ^*＜S^v+1＜1，否则0≤S^v+1≤1。

·如果F_h″(S^v+1)F_h″(S^v)＜0或F_v″(S^v+1)F_v″(S^v)＜0，对于饱和度迭代使用低松弛，即：S^v+1＝0.5(S^v+S^v+1)。

毛细管压力

在多孔介质内的多相流中，由于对岩石的润湿性能和在流体之间的表面张力，流体趋向于具有不同的压力。这种毛细管压力经常用饱和度的简单函数来模拟：

p_c(S)＝p₂-p₁    (59)

假设相1是润湿的，而相2是非润湿的。通常，毛细管压力控制小孔中微小尺寸流，但是它在储藏模拟尺寸中扮演很小的角色。因此，不希望毛细管压力改变传输方程的稳定性计算的任何特征。

具有毛细管压力的Darcy判据可被写作：

$u_1=-\frac{{\mathrm{kk}}_{r1}}{u_1}(\&dtri;p_1-{\mathrm{\&rho;}}_{1}g\&dtri;D)---\left(60\right)$

$u_2=-\frac{{\mathrm{kk}}_{r2}}{u_2}(\&dtri;p_2-{\mathrm{\&rho;}}_{2}g\&dtri;D)---\left(61\right)$

从方程(59)-(61)，

u₁＝f₁(S)u-f₃(S)(C_g-C_a)    (62)

u₂＝f₂(S)u+f₃(S)(C_g-C_a)    (63)

其中

$C_a=\frac{k}{{\mathrm{\&mu;}}_1}\&dtri;p_c---\left(64\right)$

传输方程不会因为毛细管压力的出现改变其特性。毛细管压力的影响可通过改进重力参数C_g被包含。对于总速度和重力具有不同方向的情况(即，C_g＞0的相1)，毛细管压力减小重力的影响。毛细管压力包含毛细管压力的空间梯度(

p_c)，并且希望该项相对粘度项μ或重力项C_g保持更小。如果存在任何毛细管压力，相同的稳定性算法可通过改进C_g来使用：

${\tilde{C}}_g=C_g-C_a.---\left(65\right)$

III.对于具有重力和毛细管压力的三相流的算法

现在将描述对于具有重力和毛细管压力的三相流的无条件稳定算法。在储藏模拟中，三相流区域在模型中非常有限，并且通常模型中出现的三相区域是由两个双相流的结合引起的。结果，并不清楚在实验室的三相相对渗透率测量可被直接应用于储藏模拟。此外，三相相对渗透率的测量在技术上非常困难，M.J.Oak，L.E.Barker，andD.C.Thomas，Three-phase relative permeability of Brea sandstone，J.Pet.Tech.，pages 1057-1061，1990。对于在实验测量中给出的不确定性以及在储藏模型中网格块的相对渗透率的尺度问题，对于稳定性分析的相对渗透率的简单幂律模型被使用。对于三相流的部分流量曲线的一般特征将被检查。根据部分流量曲线的特征，无条件稳定的算法对于具有重力和毛细管压力的三相流进行设计。

三相流

由于物理环境的复杂性以及在测量三相相对渗透率中的实际困难，多孔介质中的三相流并没有被很好地理解。然而，在储藏模拟中，网格单元中的三相流通常被发现主要通过双相流的合并而出现。尽管如此，这使得储藏模拟中的三相相对渗透率可以是水油或油气系统的两个相对渗透率的平均值。

改进的Stone II方法常常在许多模拟器中被使用，其中油的相对渗透率用复杂的函数形式模拟：

$k_{\mathrm{ro}}=k_{\mathrm{rocw}}[(\frac{k_{\mathrm{row}}}{k_{\mathrm{rocw}}}+k_{\mathrm{rw}})(\frac{k_{\mathrm{rog}}}{k_{\mathrm{rocw}}}+k_{\mathrm{rg}})-k_{\mathrm{rw}}-k_{\mathrm{rg}}]---\left(66\right)$

这里，k_ro，k_rw和k_rg分别是油，水和气的相对渗透率。k_rocw是在自然的水的饱和度下的油的相对渗透率，k_row是油水系统中油的相对渗透率，而k_rog是油气系统中油的相对渗透率。水和气的相对渗透率与双相系统给出的相同。改进Stone II方法的油的相对渗透率被发现对于一些系统变成负的。不足的实验证据不允许对于三相相对渗透率使用诸如Stone II的复杂非线性函数形式。A.Heiba，Three-phaserelative permeability，COFRC Technical Memorandum，1987，提出了一种k_row和K_rog的直线内插。对应实验测量中给出的不确定性和储藏模拟中网格块的相对渗透率的尺度问题，对于稳定性分析的简单幂律模型被使用。该分析可对于任意合理的三相相对渗透率的模型(无负值或非单调曲线)被进一步扩展。

具有重力的三相流将首先被检查。与上面在双相流稳定性分析中的一样，毛细管压力将通过改进重力数值在稍后被并入。

对于三相流的Darcy判据可被写作：

$u_1=\frac{{\mathrm{kk}}_{r1}}{u_1}(\&dtri;p-{\mathrm{\&rho;}}_{1}g\&dtri;D).---\left(67\right)$

这里，k是渗透率，u是相速度，k_r是相对渗透率，μ是粘度，D表示深度，而ρ是密度。下标i表示相特性(即，i＝1，2，3分别表示水相，油相和气相)。从方程(67)，可以容易地推出：

u₁＝f₁(S₁，S₂，S₃)u-f₄(S₁，S₂，S₃)C_g1-f₅(S₁，S₂，S₃)C_g2，    (68)

u₂＝f₂(S₁，S₂，S₃)u+f₄(S₁，S₂，S₃)C_g1-f₆(S₁，S₂，S₃)C_g3，and    (69)

u₃＝f₃(S₁，S₂，S₃)u+f₅(S₁，S₂，S₃)C_g2-f₆(S₁，S₂，S₃)C_g3.    (70)

总速度被定义为：

u＝u₁+u₂+u₃    (71)

而部分参数由下式给出

$f_1=\frac{k_{r1}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(72\right)$

$f_2=\frac{m_2{k}_{r2}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(73\right)$

$f_3=\frac{m_3{k}_{r3}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(74\right)$

$f_4=\frac{m_2{k}_{r1}{k}_{r2}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(75\right)$

$f_5=\frac{m_3{k}_{r1}{k}_{r3}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(76\right)$

$f_6=\frac{m_2{m}_3{k}_{r2}{k}_{r3}}{k_{r1}+m_2{k}_{r2}+m_3{k}_{r3}}---\left(77\right)$

其中，m_i是关于相1的粘度的粘度比率μ₁/μ_i。由重力引起的速度用变量C_gi表示：

$C_{g1}=\frac{\mathrm{kg}({\mathrm{\&rho;}}_1-{\text{\&rho;}}_2)}{{\mathrm{\&mu;}}_1{\mathrm{\&mu;}}_c}\&dtri;D,---\left(78\right)$

$C_{g2}=\frac{\mathrm{kg}({\mathrm{\&rho;}}_1-{\text{\&rho;}}_3)}{{\mathrm{\&mu;}}_1{\mathrm{\&mu;}}_c}\&dtri;D,---\left(79\right)$

C_g3＝C_g2-C_g1.    (80)

其中，u_c是特征速度。注意(68)-(70)中的速度u和u_i关于u_c是无量纲的。由于存在饱和度的限制：

1＝S₁+S₂+S₃.    (81)

在(68)-(70)中，只有两个方程是独立的。方程u₁和u₂被选作独立方程，同样通过用1-S₁-S₂代替S₃，S₁和S₂也被选作独立变量。

从方程(68)和(69)检查部分流表面是有意义的：

$F_1=\frac{u_1}{u}---\left(82\right)$

$F_2=\frac{u_2}{u}---\left(83\right)$

用典型的流度比率m₂＝0.2和m₃＝10以及二次相对渗透率选择一个示例，

$k_{\mathrm{rl}}=S_l^2.---\left(84\right)$

在图20和21中，部分流表面对于无重力(C_g1＝C_g2＝0)和强重力(C_g1＝5并且C_g2＝10)分别进行绘制。

部分流量曲线对于小的C_g1和C_g2是单调的，而对于大的C_g1和C_g2，它可以变成多值函数。注意依赖于重力和速度的方向，部分流量可以小于零或大于一。如果Newton非线性迭代不考虑这些表面的实际模式，这将引起数值上的困难。

从(68)和(69)，部分流量变为零的边界曲线可被获得：

0＝u-m₂C_g1k_r2-m₃C_g2k_r3，以及    (85)

0＝u+C_g1k_r1-m₃C_g3k_r3.    (86)

对于二次相对渗透率，方程(85)和(86)变为：

$0=u-m_2{C}_{g1}{S}_2^2-m_3{C}_{g2}{(1-S_1-S_2)}^2,---\left(87\right)$

$0=u{+C}_{g1}{S}_1^2-m_3{C}_{g3}{(1-S_1-S_2)}^2.---\left(88\right)$

稳定性分析

如上所述，非线性传输问题的稳定性问题在于流量函数的线性化。数值框架将被检查以判断它们是否能从初始饱和度估计收敛，

$\underset{v\&RightArrow;\&infin;}{\lim }F(S_1^v,S_2^v)=F*,---\left(89\right)$

$\left(\begin{array}{c}{F}_1^v\\ F_2^v\end{array}\right)+\left(\begin{array}{cc}\frac{\&PartialD;F_1}{{\&PartialD;S}_1}& \frac{\&PartialD;F_1}{{\&PartialD;S}_2}\\ \frac{{\&PartialD;F}_2}{\&PartialD;S_1}& \frac{{\&PartialD;F}_2}{\&PartialD;S_2}\end{array}\right)\left(\begin{array}{cc}{S}_1^{v+1}& -S_1^v\\ S_2^{v+1}& -S_2^v\end{array}\right)=\left(\begin{array}{c}{F}_1^{*}\\ F_2^{*}\end{array}\right),---\left(90\right)$

$\left(\begin{array}{c}{S}_1^{v+1}\\ S_2^{v+1}\end{array}\right)=\left(\begin{array}{c}{S}_1^v\\ S_2^v\end{array}\right)+{\left(\begin{array}{cc}\frac{{\&PartialD;F}_1}{{\&PartialD;S}_1}& \frac{{\&PartialD;F}_1}{{\&PartialD;S}_2}\\ \frac{\&PartialD;F_2}{\&PartialD;S_1}& \frac{{\&PartialD;F}_2}{{\&PartialD;S}_2}\end{array}\right)}^{-1}\left(\begin{array}{cc}{F}_1^{*}& -F_1^v\\ F_2^{*}& -F_1^v\end{array}\right),---\left(91\right)$

此外，在每次迭代之后强制0≤S^v+1≤1。

部分流表面

现在将检查部分流表面以理解它们的一般特征。在图22，23和24中，部分流量等值线分别对于无重力(C_g1＝C_g2＝0)，弱重力(C_g1＝0.5并且C_g2＝1)和强重力(C_g1＝5并且C_g2＝10)三种情况进行绘制。同时对该图使用的是m₁＝0.2和m₃＝10。如果不存在重力，每一相的部分流在与总速度相同的方向上移动。在从0.1到0.9以0.1为间隔绘制了等部分流量线的图22中没有负的流量区域。注意，相一和相二的小饱和度的大面积区域具有非常弱的部分流(＜0.1)。如果像在第二中情况中一样，重力很弱，对于相一和相二的负流量范围变得更小，并且部分流量的绝对幅度也变得更小。在图23和24的正部分流中，等部分流量曲线从0到0.8以0.2为间隔进行绘制；但是在负的部分流量区域，对于每幅图片使用不同的等高线值：对于图23A使用(-0.2，-0.4，-0.6，-0.8，-1，-2)，对于图23B使用(-0.05，-0.1，-0.15，-0.2，-0.25，-0.3)，对于图24A使用(-0.02，-0.04，-0.06，-0.08，-1，-2)，而对于图24B使用(-0.05，-0.1，-0.15)。同时绘制的还有u₁＝0(上侧的虚线)和u₂＝0(下侧的虚线)的等高线。由于在负流量范围中的复杂等高线，方程(91)的Jacobian矩阵变成奇异的。Jacobian矩阵的行列式Det变为零的奇异点位置也用实线进行绘制。因为普通的Newton迭代不能穿过这些奇异直线，好的无条件稳定算法应该考虑Jacobian矩阵的奇异点。

零部分流量边界可以对于相1(F₁＝0)和相2(F₂＝0)进行绘制：

$L_1(S_1,S_2):0=1-C_{g1}{m}_2{S}_2^2-C_{g2}{m}_3{(1-S_1-S_2)}^2,$ ，以及(92)

$L_2(S_1,S_2):0=1+C_{g1}{S}_1^2-C_{g3}{m}_3{(1-S_1-S_2)}^2.$

从方程(93)和(94)，交点可用轴S₂＝0来完成。这些点被分别表示为来自(92)的S₁ ^α1和来自(93)的S₁ ^α2。对于将在稍后描述的算法中使用的更好的饱和度的初始估计，在这些零流量曲线上定义了两个点，(S₁ ^β1S₂ ^β1)和(S₁ ^β2 S₂ ^β2)。

$S_1^{\mathrm{\&beta;}1}=2/3S_1^{\mathrm{\&alpha;}1}$ $S_1^{\mathrm{\&beta;}2}=1/3S_1^{\mathrm{\&alpha;}2}---\left(94\right)$

$L_1(S_1^{\mathrm{\&beta;}1},S_2^{\mathrm{\&beta;}1}):0=1-C_{g1}{m}_2{\left(S_2^{\mathrm{\&beta;}1}\right)}^2-C_{g2}{m}_3{(1-S_1^{\mathrm{\&beta;}1}-S_2^{\mathrm{\&beta;}1})}^2,---\left(95\right)$

$L_2(S_1^{\mathrm{\&beta;}2},S_2^{\mathrm{\&beta;}2}):0=1+C_{g1}{\left(S_1^{\mathrm{\&beta;}2}\right)}^2-C_{g3}{m}_3{(1-S_1^{\mathrm{\&beta;}2}-S_2^{\mathrm{\&beta;}2})}^2.---\left(96\right)$

一个点(S₁ ^β3 S₂ ^β3)在 $F_1^{*}<0,$ $F_2^{*}<0$ 并且Det＞0的区域中间被表示(即，图24A中的(0.3，0.3))。

用于在具有重力和毛细管压力的三相流中更新S^v->S^V+1的算法

1.检查总速度和重力的方向；

2.计算C_g1和C_g2；

3.对于相1计算零部分流量的边界(L₁(S₁，S₂)和(L₂(S₁，S₂))；

4.分别计算(L₁(S₁，S₂))和(L₂(S₁，S₂))上的饱和度初始点(S₁ ^β1，S₂ ^β1)和(S^1β2，S₂ ^β2)；

5.检查F^1*和F₂ ^*的符号；

6.从饱和度的初始估计(S₁ ⁰和S₂ ⁰)，计算部分流量估计(F₁ ⁰和F₂ ⁰)和Jacobian矩阵的行列式Det；

7.如果 $F_1^{*}\text{\&GreaterEqual;0}$ 并且 $F_2^{*}\&GreaterEqual;0;$

(a)如果饱和度的初始猜测不能得到正定，将在L₁(S₁ ^β1，S₂ ^β1)的点作为初始猜测；

(b)根据方程(91)迭代饱和度；

(c)如果 $F_1^v<0$ 或 $F_2^v<0$ ，在L₁(S₁ ^β1，S₂ ^β1)设置饱和度；

(d)如果S₁+S₂＞1，规格化饱和度以确保S_i≥0；

(e)如果 $\frac{{\&PartialD;}^2{F}_1^{v+1}}{\&PartialD;S_1^2}\frac{{\&PartialD;}^2{F}_1^v}{\&PartialD;S_1^2}<0$ 或 $\frac{{\&PartialD;}^2{F}_1^{v+1}}{\&PartialD;S_2^2}\frac{{\&PartialD;}^2{F}_1^v}{\&PartialD;S_2^2}<0,$ 低松弛饱和度迭代，即，S^v+1＝(S^v+S^v+1)/2；

8.如果 $F_1^{*}<0,$ $F_2^{*}\&GreaterEqual;0,$ 并且Det≥0，将在L₁(S₁ ^β1，S₂ ^β1)上的点的初始饱和度估计设置为初始猜测；

9.如果 $F_1^{*}<0,$ $F_2^{*}\&GreaterEqual;0$ 并且Det＜0，或者 $F_1^{*}<0,$ $F_2^{*}<0,$ 并且Det＜0，将在L₁(S₁ ^β2，S₂ ^β2)上的点的初始饱和度估计设置为初始猜测；并且

10.如果 $F_1^{*}<0,$ $F_2^{*}<0,$ 并且Det＞0，将在L₁(S₁ ^β3，S₂ ^β3)上的点的初始饱和度估计设置为初始猜测。

对于三相稳定性的可替换算法

当三相流被模拟为两个双相流区域(即，气油和油水区域)的组合时，稳定性算法可进行如下改进。双相流区域的部分流和它们的饱和度(即气油相和油水相)从三相相对渗透率近似中被确定(例如，油的饱和度在两个区域中是相同的，或者油的相对渗透率在两个区域中是相同的)。对于双相流的稳定性分析被应用于Newton-Raphson迭代中的每个区域。如果不稳定的饱和度变化在一个或两个区域中发生，进行低松弛处理以保证Newton-Raphson迭代的稳定收敛。

毛细管压力

在多孔介质内的多相流中，由于它们对岩石的润湿性能和在流体之间的表面张力，流体趋向于具有不同的压力。这种毛细管压力经常用饱和度的简单函数来模拟。在三相流中，我们有两条独立的毛细管压力曲线：

p_c1(S₁，S₂)＝p₂-p₁  以及    (97)

p_c2(S₁，S₂)＝p₃-p₁          (98)

毛细管压力控制小孔中的微小尺寸流，但是它在储藏模拟中扮演很小的角色。不希望毛细管压力改变传输方程的稳定性计算的任何特征。

具有毛细管压的Darcy判据可被写作：

$u_1=-\frac{{\mathrm{kk}}_{r1}}{{\mathrm{\&mu;}}_1}(\&dtri;p_1-{\mathrm{\&rho;}}_{1}g\&dtri;D),---\left(99\right)$

$u_2=-\frac{{\mathrm{kk}}_{r2}}{{\mathrm{\&mu;}}_2}(\&dtri;p_2-{\mathrm{\&rho;}}_{2}g\&dtri;D),---\left(100\right)$

$u_3=-\frac{{\mathrm{kk}}_{r3}}{{\mathrm{\&mu;}}_3}(\&dtri;p_3-{\mathrm{\&rho;}}_{3}g\&dtri;D).---\left(101\right)$

方程(97)-(101)得到：

u₁＝f₁(S₁，S₂)u-f₄(S₁，S₂)(C_g1-C_a1)-f₅(S₁，S₂)(C_g2-C_a2)    (102)

u₂＝f₂(S₁，S₂)u-f₄(S₁，S₂)(C_g1-C_a1)-f₆(S₁，S₂)(C_g3-C_a3)    (103)

$C_{a1}=\frac{k}{{\mathrm{\&mu;}}_1}\&dtri;p_{c1},---\left(104\right)$

$C_{a2}=\frac{k}{{\mathrm{\&mu;}}_1}\&dtri;p_{c2},---\left(105\right)$

C_a3＝C_a2-C_a1.    (106)

注意传输方程不会因为毛细管压力的出现改变其特征。毛细管压力的影响可用重力参数C_g包括。对于总速度和重力具有不同方向的情况(即，C_g＞0的相1)，毛细管压力减小重力的影响。因为毛细管压力包含毛细管压力的空间梯度(

p_c)，并且希望该项相对粘度项μ或重力项C_g保持更小。如果存在任何毛细管压力，第II章的相同算法可通过改进C_g来使用：

${\tilde{C}}_{g1}=C_{g1}-C_{a1}.---\left(107\right)$

Claims (14)

1.一种通过求解非线性S形函数F＝f(S)用储藏模拟器进行建模的方法，该函数表示其中具有流体流的物理地下储层系统中的特性S，所述S形函数具有拐点S^c以及f(S)收敛到的预定值F，所述方法的特征在于包括以下步骤：

a)选择初始值S^v+1作为函数f(S)的可能解；

b)在S^v处对函数f(S)执行Newton迭代(T)以确定下一迭代值S^v+1；

c)判断S^v+1是否相对S^v位于拐点S^c的另一侧；

d)如果S^v+1相对S^v位于f(S)的拐点的另一侧，设置S^v+1＝S¹，一个改进的新估计；

e)判断S^v+1是否在预定的收敛判据内；

f)如果S^v+1不在预定的收敛判据内，设置S^v＝S^v+1；

g)重复步骤b-f，直到S^v+1在预定的收敛判据内；以及

h)选择S^v+1作为S形函数F＝f(S)的令人满意的解S；并且

i)利用所述S形函数F＝f(S)的令人满意的解S，用储藏模拟器对地下储层内的流体流进行建模。

2.权利要求1的方法，其中改进的新估计S^l被设置为S^l＝αS^v+(1-α)S^v+1，其中0＜α＜1。

3.权利要求2的方法，其中0.3＜α＜0.7。

4.权利要求2的方法，其中α＝0.5。

5.权利要求1的方法，其中改进的新估计S^l被设置为函数f(S)的拐点S^c。

6.权利要求1的方法，其中S形函数f(S)是在储藏模拟器中使用以代表地下储层内的流体流的特性的部分流量函数。

7.权利要求1的方法，其中S形函数f(S)在数学上符合下面的数学表达式：

$f\left(S\right)=\frac{\frac{S^2}{{\mathrm{\&mu;}}_1}}{\frac{S^2}{{\mathrm{\&mu;}}_1}+\frac{{(1-S)}^2}{{\mathrm{\&mu;}}_2}}$

其中S＝双相流体流量问题中第一相的饱和度；

f(S)是饱和度S的函数；

μ₁＝第一相中的流体粘度；而

μ₂＝第二相中的流体粘度。

8.权利要求1的方法，其中，如果ε＞|S^v+1-S^v|，则S^v+1在预定的收敛判据内，这里ε是预定的容许极限。

9.权利要求1的方法，其中，如果ε＞|F-f(S^v+1)|，则S^v+1在预定的收敛判据内，这里ε是预定的容许极限。

10.权利要求1的方法，还包括：将S^v+1设置为上或下界S^b，如果在步骤b中S^v+1超过这些界限的话。

11.权利要求1的方法，其中所述S形函数F＝f(S)由查找表和普通解析函数之一来表征。

12.权利要求1、2、3、4、5、6、7、8、9、10或11的方法，其中如果在储藏模拟器的网格单元中部分流量函数f(S)与物理地下储层系统内的具有重力的双相流相关联，则更新Newton迭代，所述更新特征在于包括：

(a)检查网格单元中相速度和重力的方向；

(b)如果相速度是并流，应用常规的Newton稳定性分析；

(i)计算定向的部分流量F_h和F_v，以及它们的二阶导数F_h”和F_v”；

(ii)计算垂直流的C_g，对于水平流，C_g＝0；

(iii)应用饱和度限制0≤S^v+1≤1；

(iv)如果F_h″(S^v+1)F_h″(S^v)＜0或F_v″(S^v+1)F_v″(S^v)＜0，对于饱和度迭代应用低松弛；

(c)如果相速度是逆流；

(i)计算反射饱和度S_r ^＊；

(ii)应用饱和度限制：如果所述单元不邻近边界或者相邻单元并不属于不同饱和度域，那么0≤S^v+1≤S_r ^＊或S_r ^＊＜S^v+1＜1，否则0≤S^v+1≤1；

(iii)如果F_h″(S^v+1)F_h″(S^v)＜0或F_v″(S^v+1)F_v″(S^v)＜0，对于饱和度迭代应用低松弛；并且

(d)利用更新后的Newton迭代，用储藏模拟器对物理地下储层系统的双相流进行建模。

13.权利要求1、2、3、4、5、6、7、8、9、10或11的方法，其中如果在储藏模拟器的网格单元中部分流量函数f(S)与物理地下储层系统的具有重力的三相流相关联，则更新Newton迭代，所述更新特征在于包括：

(a)检查总速度和重力的方向；

(b)计算C_g1和C_g2；

(c)对于相1计算零部分流量的边界L₁(S₁，S₂)和L₂(S₁，S₂)；

(d)分别计算L₁(S₁，S₂)和L₂(S₁，S₂)上的饱和度初始点

和

(e)检查和

的符号；

(f)根据饱和度的初始估计S₁ ⁰和S₂ ⁰，计算部分流量估计F₁ ⁰和F₂ ⁰以及Jacobian矩阵的行列式；

(g)如果

并且

(i)如果饱和度的初始估计没有得到正定，将在的点作为初始饱和度估计；

(ii)根据下式迭代饱和度：

$\left(\begin{array}{c}{S}_1^{v+1}\\ S_2^{v+1}\end{array}\right)=\left(\begin{array}{c}{S}_1^v\\ S_2^v\end{array}\right)+{\left(\begin{array}{cc}\frac{{\&PartialD;F}_1}{{\&PartialD;S}_1}& \frac{{\&PartialD;F}_1}{{\&PartialD;S}_2}\\ \frac{{\&PartialD;F}_2}{{\&PartialD;S}_1}& \frac{{\&PartialD;F}_2}{{\&PartialD;S}_2}\end{array}\right)}^{-1}\left(\begin{array}{cc}{F}_1^{*}& -F_1^v\\ F_2^{*}& -F_1^v\end{array}\right),$

其中在每次迭代后令0≤S^v+1≤1；

(iii)如果或

在设置饱和度；

(iv)如果S₁+S₂＞1，规格化饱和度以确保S_i≥0；

(v)如果

或

低松弛饱和度迭代；

(h)如果

并且Det≥0，设置在

处的点上的初始饱和度估计；

(i)如果

并且Det＜0，或者

并且Det＜0，设置在

处的点上的初始饱和度估计；以及

(j)如果

并且Det＞0，设置在

处的点上的初始饱和度估计；并且

(k)利用更新后的Newton迭代，用储藏模拟器对物理地下储层系统的三相流进行建模。

14.一种通过求解非线性S形函数F＝f(S)用储藏模拟器进行建模的设备，该函数表示其中具有流体流的物理地下储层系统中的特性S，所述S形函数具有拐点S^C以及f(S)收敛到的预定值F，该设备的特征在于包括以下部件：

a)选择初始值S^v作为函数f(S)的可能解的部件；

b)在S^v处对函数f(S)执行Newton迭代(T)以确定下一迭代值S^v+1的部件；

c)判断S^v+1是否相对S^v位于拐点S^c的另一侧的部件；

d)如果S^v+1相对S^v位于f(S)的拐点的另一侧，设置S^v+1＝S¹，一个改进的新估计的部件；

e)判断S^v+1是否在预定的收敛判据内的部件；

f)如果S^v+1不在预定的收敛判据内，设置S^v＝S^v+1的部件；

g)重复部件b-f的动作，直到S^v+1在预定的收敛判据内的部件；

h)选择S^v+1作为S形函数F＝f(S)的令人满意的解的部件；以及

i)利用所述S形函数f(S)的令人满意的解，用储藏模拟器对地下储层内的流体流进行建模的部件。

Images