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

专利名称：使用改进Newton-Raphson算法求解S形非线性函数的稳定方法和设备的制作方法
技术领域：
本发明一般地涉及使用Newton-R叩hson算法的求解方法，更具 体地，涉及用于求解诸如与储藏模拟器中的传输问题相关的S形函数 的非线性方程的隐式求解算法。
背景技术：
隐式求解算法被广泛应用于解决在诸如储藏模拟的计算机流体 动力(CFD)应用中的非线性时间相关的传输问题。这些隐式求解算 法通常被用来克服显式方法的时步大小限制。这些传输问题在新的时 间层上的解常常4吏用著名的Newton- Raphson方法来获得。
图1简要地说明了用于求解或至少是近似函数F =/(5)的解的 Newton-Raphson方法的一次简单执行。图l描绘了下列函数的图形
(1)
求解/(5)=尸=0确定了一个根，即，函数/(s)与水平轴F-0相 交的地方。S的截断位置用字母A表示。如果该位置，或者到其上的 令人满意的近似可被确定，那么可以认为方程(l)被求解。
Newton-Raphson方法要求对解j故一个初始猜测。该初始猜测在 图1中被记为510。在曲线尸=/^)上的对应点8具有坐标(^， /OS0))。 在点B处画一条与水平轴相交的切线。该切线的斜率是/'(W)。
代数上，在普通点5^上的切线方程是 尸-/(50)-/'(50)05 —圹) (2)
解S是F-0时切线的水平坐标。求解该坐标
lS = lS°-/(tSI°)//'(S°) (3)
该坐标随后被用作/(^ = 0的解的第二猜测或近似，并被记为乂。因此，
5"'-W—/(^)//'0S1。) (4)
对方程(1)的根的更好的估计可通过再次获得在坐标0SV(5"1)), 即，点C上的切线，并确定它与水平轴的交点来确定。这可是用下式 来实现
S"S1-/(S'),/'(S'〉 (5)
可以看到f相比y是对方程(l)的解的更好近似。通过获得产生
在^处的截距的在点D， (SV(^))上的切线可以得到更好的估计。总 的来说，如果对F =/(5)的解的初始近似或猜测是51",下一近似5"+1
可用下式确定
S"+'= s"-f(S")/f(S"〉 (6)
方程(6)提供了 Newton-Raphson方法的基本公式。通过对函数 f(S)重复使用Newton-Raphson方程(6),对于方程(l)的越来越好的近 似可被确定。然而，这假设用来找寻/(5)的解的Newton-Raphson方 法的使用是无条件稳定的。事实并不总是如此。
在储藏模拟中，待求解函数通常是传输问题，它在解的最优估计 附近被线性化，并在每次迭代中被隐式求解。在一个典型的储藏模拟 器中，该传输问题在每次Newton-Raphson迭代中对于储藏才莫型中的 每个单元进行求解。 一旦在多次Newton-Raphson迭代之后获得了 一 个时步的收敛解，该模拟随后继续进行下一时步。如果流量函数如图 1所示是凸的(冲击)，凹的(稀疏波)，或线性的(接触间断)， 这种方法简单并且无条件稳定，这是大部分计算机流体动力(CFD) 问题的情况。然而，对于在大部分储藏模拟中遇到的S形流量函数， 如果所选时步过大，Newton-Raphson方法不收敛。在储藏模拟中， 这个问题通常用经验性时步控制技术来克服。然而，得到的时步大小 常常是保守的或过大的，导致了多余的计算或时步缩减。
因此，需要一种用于求解具有S形流量函数的传输问题的隐式求 解算法，它不管时步大小无条件稳定。该发明提供了对于这个需求的 解决方案。

发明内容
提供了一种用于求解非线性S形函数F-Z(^)的方法，该函数表
示物理系统中的特性S。
S形函数具有拐点Se和/(5)收敛到的预定值
尺初始猜测^被选作函数/(5)的可能解。Newton迭代(T )在 的^处被执行以确定下一迭代值^"。接着判断S^是否相对Sv位于 拐点se的另一侧。如果SV"相对SV的位于拐点的另一侧，5"+1被设置 成^'， 一个改进的新的估计。然后判断^"是否位于预定的收敛判据 内。如果不是，那么5"被设置成5"+1,并重复上述步骤直到5"+1位于 预定的收敛判据内。随后，5"+1的收敛值被选作S形函数F-y(^)的令 人满意的解。改进的新的估计^最好被设置成拐点Se，或者^和5"+1 之间的平均值，即，^ = 0.5(5"+ 5"+1)。可替换地，除了0.5之外的松 弛因子可被使用，例如从0.3至0.7中选择的值。
此外，求解算法对于具有重力和毛细管压力的双相和三相流进行
描述o
本发明的目标是提供一种使用无条件稳定的改进
Newton-Raphson框架求解S形非线性函数以得到基于物理系统的问 题的近似解的方法和设备。


本发明的这些及其它目标，特点和优点将参考下面的说明书，待 审批的权利要求和附图变得更好理解，其中
图1清晰地描绘了使用传统Newton-Raphson方法的函数F =/(5) 的迭代求解方法；
图2是对于传统储藏模拟器中一个时步的流程图，其中算子T 表示使用Newton迭代求解线性化的传输函数；
图3对于包含注入和生产井位置的2D测试情况显示了 30x30网 格的渗透场A:;
图4A, 4B和4C显示了在几个初始时步之后，来自图3的渗透场A的饱和度分布；
图5A和5B显示了凸和凹的函数/(5)各自的收敛图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形函数的收敛图IOA，IOB,IOC，IOD,IOE,IOF，IOG，IOH和101分别显示了对于 三种粘度比/V/^和三个时步大小，在0.5个注入的孔隙体积倍数 (PVI)之后的饱和度场。
图11显示了使用第一和第二改进Newton框架以及使用具有常 数松弛因子的Newton框架在0.5个PVI并且/z,/^-l之后的收敛历
史；
图12是对于一个时步具有非常数速度场的流程图，其中求解压 力方程用算子P表示，而求解线性化的传输方程用算子T表示；
图13显示了对于三种模拟的一个时步的收敛历史，两种使用改 进Newton-Raphson框架执行，第三种在所有迭代步骤中使用常数松 弛因子；
图14是在本发明的优选实施例中执行的用以求解对于诸如饱 和度的物理特性S的S形函数F-y(A中的步骤的流程图15清晰地显示了使用根据本发明的改进Newton-Raphson框 架的8形函数/^=_/(5)的迭代求解过程；
图16A和16B显示了相1的部分流量曲线(A)m-l而(B) m
=0.2;
图17A和17B显示了相2的部分流量曲线(A) m = 1而(B) m
=0.2;图18A和18B描绘了在单元/和w之间的部分流量在上侧图 中，单元/是逆风的，而在下侧图中，单元w是逆风的，(A)无逆流； (B)逆流在单元/中；(C)逆流在单元w中；以及(D)逆流在单元/和w 中；
图19A和19B显示了具有逆流和并流的油相(较轻的流体)的 部分流量曲线；
图20A和20B显示了对于(A)相1和(B)相2: m产0.2, m3=10 并且Cgl= Cg2=0的没有重力影响的部分流表面；
图21A和21B显示了对于(A)相1和(B)相2: m尸0.2, m3=10. Cgl= 5， Cg2=10的有重力影响的部分流表面；
图22A和22B描绘了对于(A)相1和(B)相2: m产0.2， m3=10. Cgl=0， Cg2=0的没有重力影响的部分流表面；
图23A和23B描绘了对于(A)相1和(B)相2: m产0.2， m3=10. Cgl=0.5并且Cg2=1.0的具有小的重力影响的部分流表面；
图24A和24B描绘了对于(A)相1和(B)相2: m尸0.2， m3=10. Cgl= 5并且Cs2=10的具有大的重力影响的部分流表面；
具体实施例方式
I.储藏模拟的理论背景和数值研究 考虑下面的双曲线方程
- 0 on n (7)
方程(7)是表示在双相表层流问题中单相传输的模型方程。因变 量S是该相的饱和度(OS^l) ， /是^的函数，0$/^l，而u是无散 度总速度场。为了简单而不失一般性，u被认为是常数，并且假设
公式(7)的非线性传输问题可用如图2的流程图所示的传统 Newton- Raphson方法迭代地求解。该过程得到新的时间层上的解 5*n+1。 S形函数的解的初始猜测y被给出。该猜测可以是来自先前时步的饱和度值^。在每次迭代V中，线性化传输方程最好对于储藏模
型或子模型中的每个单元进行求解(用算子T表示)，得到一个更新 的饱和度场5"+1。 一旦从一次迭代到下一迭代的最大绝对饱和度变化 没有超出一个适当定义的门限￡，这个非线性循环被认为收敛，并且 犷+1构成了在新的时间层上的解5"+1。 研究方程(7)的2D离散形式
v+lw A, f "+l "+l 、 A,广"+l"+l 、 (8、
3 、 乂 、
这是基于隐式Euler和一阶展开。为了简单起见，方程(8)对于速 度向量u-^"/的两个分量均大于等于零，并假设正交的均匀间隔的
网格的情况进行书写。每个网格单元可用索引对,V确定，并且在两个 坐标方向上的网格间距被定义为Ax和Ay。上标《和v分别表示旧的时 间和迭代级别，而A/是时步大小。为了求解方程(8)，方程(8)中的流量 /"+1被线性化近似
W+l V+l V
、_^
77
V+I V
S —S
(9)
最后，的近似5"+1通过求解线性化系统被获得
w十",x+'《o+s":w -x-—-"《- a》
它是通过将方程(9)代入方程(8)得到的。在每次求解方程(10)时强
制0SSV+1￡1是很重要的。
对于在本说明书中提到的数值研究，30x30的网格单元被用于建 模。图3显示了用于通过求解压力/;的椭圆方程计算总速度『-;iAp的 渗透场A:和注入和生产井位置
▽.(AVp)-—鋒D (11)
这里，人是总流度义^& + &]，其中，^是相y的相对渗透率，
、M 〃2 JA是该相的粘度。源项《对生产井是1，对注入井是-1，并且SeO是
饱和度的初始条件。对于简单二次方程的情况，典型的s形部分流量 函数可被表示为
匕
= A (12)
--J--
其中，M和^是两相的粘度(注入相用下标l标记)。为了说明
上面参考图2解释的隐式求解算法对于过大的时步是不稳定的，给出 了三个数值实验的结果，对于这些实验，//,/^被分别选取为1， IO和
0.1。图4显示了在少许初始时步之后的饱和度分布。尽管时步相对注 入的孔隙体积倍数(PVI)较小，该结果显示了以饱和度前端的棋盘 形图案为特征的强烈的振荡行为。图2的求解算法的非线性循环不收 敛，并且迭代次数被限制为200。现在研究这种不稳定性或者振荡行 为的起因，并描述可得到对于Newton-Raphson方法的使用无条件稳 定的框架的变型。
理解非线性传输问题的稳定性问题的关键在于流量函数的方程 (9)的线性化。首先，传输问题f(S)的非线性被分离。然后判断对于这 些在0SS、1范围内的初始值，迭代框架是否满足收敛条件
/0SV) + /'CS^)0 v +1 — =尸. (14)
在y周围的f(S)。 Newton迭代乂>式如下获得
"+ 1="尸-/(。 (15) /'(Sv)
此外，在每次迭代之后强制0SSS1。现在研究条件(13)满足的S° 和F的组合。如果函数f是S的线性，凸的或凹的函数，条件(13)总 是满足。由于线性情况很简单，只给出了凸的和凹的情况的例子，分别用f(S"s2和f(S)-l-(l-S)2表示。
图5A和5B显示了这些凸函数和凹函数的收敛图。收敛图通过 选择初始值SG和流量函数f的目标值F被创建。然后迭代次数被确定， 它要求收敛到距目标值F预定的容许偏差的范围内。水平轴表示起始 值S、而垂直轴表示目标值F。该处理在目标值F处对于大量起始值 O^S^l进行重复。接着该处理对于大量(^F2的目标值进行重复。 随后，对于每组值(S， F)达到收敛所需的迭代次数的确定被画出以 产生收敛图。
快速收敛被表示为较暗的阴影，而白色显示了参数空间中迭代方 程(15)不收敛的区域。在图5A和5B中不存在白色区域，因此该迭代 框架对于凸函数和凹函数f(S)是无条件稳定的。
然而，如果流量函数/是S形函数，情况就不同了。 S型函数的 特征是具有在拐点Sc处相连的凸部分和凹部分，如图6A和56B所示。 拐点Se的特征是函数f(S)的二阶导数f，(S)在其周围改变符号的点。
图6A和6B中，上侧图分别显示了对于//1//^2 = 1和//1/^2=0.1由方
程(12)定义的/。对应的收敛图4皮显示在图6C和6D中。大的白色区 域的存在清楚地显示了应用于方程(12)的流量函数/的迭代框架对于 这些^///2值不是无条件稳定的。如果SQ的初始猜测很接近最终解， 即，如果|/(^-。l不太大，则能够达到收敛。如果求解传输问题的话，
这与限制时步大小或最大饱和度变化是一致的。
图6C和6D中的收敛图显示了如果S、 Se,迭代方程(15)对于所 有可能目标值F均收敛，其中Se是/的拐点。这个发现表明如果fCr) 和"5"+1)的二阶导数具有相反的符号，设置5*v+1= S、如果这在Newton 迭代的最后被完成，改进的Newton-Raphson框架4皮产生，它对于S 形流量函数是无条件稳定的。
图7A和7B显示了对于图6A和6B的S形流量函数的收敛图。
不管怎样，如果cr)和(5"。的二阶导数具有相反的符号，设置y"-^r，
这就提供了第一改进Newton-Raphson框架。这可通过计算下式进行 校验
13f'( Sv)f'( Sv+1)<0 (16)
如果二阶导数的乘积是负的，那么5"和5"+1在拐点Se的两侧。
水平轴表示起始值S、而垂直轴是目标值F。暗色区域表示快速 收敛，而不收敛对于白色区域中的参数对被获得。注意，在图7A和 7B中不存在白色区域，因此第一改进Newton-Raphson框架是无条件 稳定的。将该改进应用于传输问题的求解算法上是直观的，并且允许 任意大的时步。
确定拐点Se的精确位置是很难的。如果求解部分流量函数依赖 于多个变量的传输方程的系统，这就更加真实。在某些迭代中，查找 表可被用来表示部分流量函数f(S)，并且拐点不能被简单地解析地计
算。此外，找到复部分流量函数的拐点sQ是计算复杂的。
因此，可选择地，并且最好是提出在现实设置中更容易实现的第
二改进框架。第二改进框架在每次Newton迭代结尾时S"1的改进方
式上与第一框架不同，其中y和^"位于曲线/的拐点^的两侧。代
替将^+1设置为5"，可使用一种有选择性的低松弛法。例如，如果二
阶导数，即/(5) = 32//迅2的符号在新旧迭代级别上不相同，Sv+1可用
0.5(5"+犷+1)代替。其它级别的松弛也可被使用。也就是说，低松弛可
由下式推广
S"o^+(1-a『1 (17)
其中，S^S的新的估计；而01=低松弛系数。
第二增强框架的流程图在图8中描绘。由于确定/的二阶导数相 对容易，与当5"和5""在y的两侧时将sv"设置为^的第一改进框 架相比，第二可替换框架可以在现有的储藏模拟器中被容易地实现。
与图6C, 6D和图7A， 7B中显示的收敛图相似，使用可选低松 弛的第二可替换改进对应的收敛图被显示在图9A和9B中。水平轴表 示起始值5^，而垂直轴表示目标值F。暗色区域表示快速收敛，而不 收敛对于白色区域中的参数对被获得。与图7A和7B中的收敛图 一样， 没有白色区域显示了第二可替换改进得到了一个无条件稳定的框架。
如果被应用于非线性传输方程的求解算法，第一和第二改进均用作预测。第一和第二改进框架几乎是等效的。此外，如果时步大小在
未改进或传统Newton-Raphson系统的稳定范围内，由本发明的改进 引起的效率上的损害没有在数值研究中被发现。另一方面，如果通过 在每次Newton迭代结尾的任何地方均使用低松弛来获得稳定性，而 不是如上所述在改进框架中可选地使用，那么效率会被严重损害。对 于使用由方程(12)定义的S形部分流量函数的上述测试情况，数值结 果被显示在图10A-I中。所有模拟用第一和第二改进框架来执行，并 且证实任意大的时步可被使用。图10A-I分别显示了在时间-0.5PVI 之后当V^"(顶行——图IOA， 10B和10C) ， //,//^"0 (中间行
-图IOD， 10E和10F)，以及//,///,0.1 (底行-图10G, 10H
和101)时获得的饱和度场。在这些图左侧列中的结果使用100倍小 的时步，或A,-o.oo5PVi来计算，对此未改进的Newton-Raphson框架 证明也是稳定的。注意这些时步是用来获得如图4所示的不收敛结果
的大小的 一半。图的中间和右侧列分别显示了 Ak 0.1PVI和A, = 0.5PVI时
计算的结果。
图11显示了使用上述第一和第二改进Newton-Raphson框架的 结果，其中只有当二阶导数/"『)和/" 的符号在迭代步骤中改变 时，5""被分别可选地设置为^和0.5(5" + 5*v+1)。进一步，研究在其 中低松弛被应用于每次迭代的结尾的第三框架，即单纯低松弛框架。 垂直轴显示了在迭代步长中的最大饱和度变化，而垂直轴显示了迭代 的次数。图11对于时步大小是0.5PVI，并且粘度比率a/^是一的情
况进行绘制。对于执行的所有模拟，第一和第二改进框架均被测试， 并且在所需迭代方面没能观察到明显的差异。虽然前两个框架的收敛 率基本相同，使用简单低松弛的第三框架被显示更为低效。
如果流度因子人是饱和度的函数，改进Newton-Raphson框架也 可被用作求解算法的积木。这对于石油行业特别重要。对于二次相对 渗透的情况，流度因子X被定义为
、 乂
15其中k是渗透率； S=第一相的饱和度
Hi=第一相的粘度；而 第二相的粘度。
作为k是饱和度的函数的结果，速度场u不再是常数，并且如果 X变化，需要通过求解压力方程(18)进行更新。该算法通过图12中的 流程图进行概述。注意图8的流程图再次出现在更新整个速度场u所 需循环中的虚线框内。这个收敛解与通过完全隐式的框架获得的解是 一样的，其中压力和传输方程被同时求解。
图13显示了对于两种模拟的一个时步的收敛历史， 一个由使用 可选低松弛的改进Newton-Raphson求解算法执行，另一个由使用简 单或单纯低松弛的方法执行。测试情况与前面的研究相同，时步大小 是0.5PVI，而粘度比/Va"。使用改进算法，求解线性化传输方程 需要大约60次，而椭圆压力方程只需要九次(收敛图中的峰的个数)。 简单或单纯低松弛在另 一方面被证明是更低效的。
用于求解对于物理特性的S形函数的典型步骤
本发明包含用于稳定求解非线性S形曲线F=f(S)，以确定物理系 统的特性S的表示的方法和设备。图14显示了可被用来依照上面章 节描述的原理实现本发明的步骤的优选示范性流程图。图15清晰地 描述了改进Newton框架。
在改进Newton-Raphson算法中，在S形函数的二阶导数/"0T)和 /"(f )的符号在相继迭代之间改变了的情况下，下一迭代值^+1被设 置成新的估计或稳定性限制S1。改进的新的估计，或稳定性限制5"1 被理想地选择以使得改进Newton- Raphson方法保证稳定收敛。
第一步骤110获得表示物理系统中的特性S的S形函数f(S)。例 如，特别感兴趣的物理系统是地下储藏中的流体流量。S形函数f(S) 可用数学函数表示，例如上面的方程(12),或使用查找表的方式获取 S形函数的要素。像储藏模拟领域的技术人员所知的那样，诸如部分 流量函数的S形函数可根据相对渗透率度量从实验室测试中得到。求解对于物理系统的非线性问题领域的技术人员将会认识到使
用改进Newton-Raphson框架的本发明可被用来求解表示除储藏模拟 之外的物理系统的特性的特征的S形函数。
S形函数f(S)通过包含至少一个在拐点Se处相连的凸部分和凹部 分来表征。函数f(S)被求解以确定满足条件F-f(S)的物理特性S，其 中F是目标值。例如，在储藏模拟中，当求解在线性化传输问题的网 格单元中的饱和度值S的S形部分流量函数f(S)时，来自先前时间层 的流量值F被理想地用作目标值F。可替换地，目标值F也可以通过 数值模拟中的非线性迭代来估计。
初始值或猜测5"在步骤120中被选作满足条件F = f(S)的物理特 性S的可能解。^的初始值最好用初始条件或来自先前时步的结果进 行估计。尽管次选的，^"的初始值也可用内插方法或随机化方法得 到。
在图8中用算子(T)表示的Newton迭代在步骤130中使用当前迭 代的值5"在函数f(S)上被执行以确定下一迭代值Sv+1。该Newton迭 代对应于将方程(15)应用于S形函数f(S)以确定下一迭代值Sv+1 。在图 15中，初始猜测是S、
可选的检测可在步骤140中进行以查看5""是否在可接受的范围 内。例如，如果在诸如与方程(12)中类似的函数的S形函数的分母接 近零，对于5"+1的估计值可能相当大。如果^"超过某一边界限制， 例如饱和度S小于零或大于一，那么可在步骤150中将5""设置为诸 如0.0或1.0的值以保持物理特性S在现实范围内。在图15的例子中， 可以看到点(S。，f(sO))的正切与表示流量值F的水平线的截距没有出现 在饱和度上限Sb-l.O内。因此，5"+1在图15中第一迭代的情况下被 设置为上限S、1.0。
然后在步骤160中判断^T和5""是否位于S形函数f(S)的拐点 Se的两側。进行这种判断最好的方式是查看f(S)二阶导数/(S)在Sv 和Sv+1处的符号是否彼此相反。这可用公式(16)计算/("和/(f)的 乘积来检验。如果二阶导数的乘积是负的，那么^和5""位于拐点的两侧。 尽管是次优的，其它进行这种判断的方式是明确计算拐点S^的位置，
然后查看5"v和^"是否在的Se的同一侧。拐点&可通过求解/'(。 = 0
来计算。这在f(S)足够简单时可解析地实现，而在f(S)比较复杂时在 数值上实现。在使用查找表定义S形函数f(S)的情况下，f(S), /'0S)和 /(。的初始值表可被产生。然后查找表可被检查以确定拐点Sc的位 置。在图15中，S。和Sb的值位于拐点Se的两侧。
如果5"v和5""位于拐点Se的两侧，那么为了保证在求解S时无 条件稳定，在步骤170中，^+1被设置成新的估计或稳定性限制511。 该改进估计^被选择以保证Newton-Raphson迭代的收敛。在第一优 选实施例中，这个改进估计^被设置成等于拐点Sc。在如图15所示 的第二或更好的方式中，低松弛方法被用来计算S、理想地，^根据 下面的方程进行计算
=必"+(卜牵("+|) (〗9) 其中0<a<l。
优选地，a的范围是0.3<a<0.7，更优选地，0.4<a<0.6,而最优 选地，a=0.5。使用a=0.5, ^+1或^被设置成0.5(^和^+1)。
如果判断5"和^v+1位于拐点Se的同 一侧，那么Sv+1最好保持不 变而不是被设置成改进估计S1，即y或sL0.5(5"+5^1)。
在步骤180中，判断当前迭代值^"是否达到预定的容许条件。 例如，可以要求在5"和的连续迭代值之间的差足够小。
数学上，这可被表示为 ￡>1^+1-叫 (20)
其中，￡是预定的容许极限。
可替换地，为了考虑到解5"v+1可以足够接近真实解，也可以要求
e>n+')i (2i)
如果^"不被认为是一个令人满意的解，那么在步骤190中^ 被设置成Sv+1，并且继续步骤130-180直到^"的值收敛到预定的容 许极限内。图15举例说明了对于饱和度S1， 82和83等的连续估计如
18何收敛到饱和度s的真实值。如果达到充分的收敛，物理值s被认为
是令人满意的解。在储藏模拟器的情况下，当对于在其上执行改进
Newton迭代的每个网格单元达到收敛之后，对于下一时间层Sn+1在 网格单元中的饱和度s随后被设置成5*v+1。
本发明还包括设备或机器可读的确实包含指令程序的程序存储 设备。这些指令被机器执行以实现稳定求解表示物理系统中的特性S 的非线性S形函数F = f(S)的方法步骤。该S形函数具有拐点Se以及 f(S)收敛到的预定值F。所述方法包括下列步骤。选择F-f(S)的可能 解的初始值。在函数f(S)上在5"处执行Newton迭代以确定下一迭代 值5"+1。然后判断5"+1是否相对5"位于拐点Se的另一侧。如果^和 ^"位于相对侧，那么5""被设置成等于改进的新估计5"。接着进行 检查以通过判断5""是否在预定收敛判据内来判断5"+1是否是特性S 的令人满意的解。如果5""没有令人满意地收敛，那么5"被设置成等 于5"+1。重复上述步骤直到5"+1在预定的收敛判据内。当达到收敛时， Sv+1被选作问题F = f(S)的令人满意的解。
II.对于具有重力和毛细管压力的双相流的算法
重力和毛细管压力在Newton-Raphson方法的数值稳定性上的 影响将对于双相流的非线性传输方法进行测试。重力和/或毛细管压力 在储藏模拟中的影响可导致储藏单元中的逆流。这种逆流在传统的储 藏模拟器中会造成严重的数值不稳定。通过扩展上面在第I章中描述 的方法，将得到对于具有重力和毛细管压力的双向流的非线性传输方 程无条件稳定的算法。
首先，将在部分流量曲线上研究重力的影响。然后用特征参数描 述部分流量曲线的特征。接着描述对于具有重力的双向流的无条件稳 定算法。改进也可被提供给无条件稳定算法以包含毛细管压力的影 响。
重力影响
对于具有重力的双向流的Darcy判据可被写作这里，A是渗透率，"是相速度，/^是相对渗透率，；/是粘度，D 表示深度而p是密度。下标表示相特性。
从方程(22)和(23)，可以容易地推导出
"2=/2(S)"+/3(S)Cg (25)
其中S表示相1 (即，水)的饱和度。总速度被定义为
w - w, + w2 (26)
,=F^V (27)
,2 ='"、 ■ (28)
+《2
,=W^,r2 (29) 3 、+怖、2
这里，m是粘度比率//,///2。
重力引起的速度可用变量G表示，
c ,(A-A)VD (3Q)
其中，"c是特征速度而D是深度。注意，方程(24)和(25)中的速
度《，"/和"2关于Wc也是无量纲的。
从方程(24)和(25)分析部分流量曲线是有意义的。在图16和17 中，部分流量分别对于相1和相2进行绘制，
(31)
F2-&, (32)
具有二次相对渗透率， 、《 (33)如果(XwC^1，部分流量曲线是单调的，反之，它对于wCV〉1 或附Q o可能变成多值函数。注意部分流量可以小于零或大于一， 这取决于重力和速度的方向。如果Newton非线性迭代不能考虑这些 曲线的模式，这会导致数值上的困难。注意来自负Cg和正"的部分 流量曲线等价于来自正Cg和负《的情况。在实际的具有固定坐标的 模拟中，G可保持恒正，但是平均速度可正可负。注意不失一般性，
的具有负G正平均速度的情况。
如果油(较轻的流体)方程;故选作初级方程，从方程(25)，油相 的部分流变成零S(/或一 S/的状态
以及
<formula>formula see original document page 21</formula>(35)
<formula>formula see original document page 21</formula> (36)
如果对于 i ，部分曲线变成单调曲线；反之如果M〈。，
或对于 w 〉 0 ， /|"| 〉 1 ， 部分流量曲线包括部分曲线变成多值(非单
调)的区域。当部分流量曲线具有最小值或最大值点时，那么这里存
在在非线性迭代中饱和度并不穿过的拐点(参见图19)。如杲相对渗 透率类似方程(33)和(34)中的二次方，该拐点可被解析地计算。
三次方程可被解析求解，并且总能获得至少一个实数根<formula>formula see original document page 21</formula>(38) 其中<formula>formula see original document page 21</formula>—; (43)
饱和度可被分割成两个范围(1) S〉5^和S〈S二总的来说， 饱和度不穿过这个反射点，但如果该单元在无流边界附近，或者挨着 具有其它范围的饱和度的单元，该状态应该被允许穿越。由于存在允 许流体在没有逆流的情况下进行累加的边界，部分流量曲线并不适 用。
有限差分近似
图16和17中的部分流量曲线表明双相流可以是并流和逆流，取 决于特征变量G，粘度比率附和饱和度51。如图18所示，当单元/ 和邻近单元n中的部分流没有引起逆流时，普通的有限展开差分近似 可被明确地使用。然而，如果逆流在单元/或单元n中发生，如图18 所示，上行流是由相决定的。
作为结果，从在压力解决方案中计算的总速度分离的相带有在单 元/或n之间的并流或逆流的详细信息。在图18中，描述了能够在单 元/和它的邻近单元n之间发生的所有可能情况。假设来自压力解决 方案中的总速度"/是在单元/和n的分界面上的速度，并且在单元 和n内的总速度是一致的，称之为"。
在单元Z中，并向和对向的相速度分别用"/'和"/'表示其总速度 。这些相速度可从方程(24)和(25)中获得。该方程可被表示为 < ，C,), (45)
类似地，单元n中的相速度可通过下式获得 ";- Cg，w)， (47)
22如果单元/是上行流，
如果单元n是上行流，
(49)
(50)
注意假设某一相速度可以是并向或对向的，这依赖于流和重力方 向。此外，如果0〈wC^1，逆流并不出现。
在顺序隐式算法中，饱和度方程使用固定的总速度来迭代求解， 以保持质量守恒(即，发散自由度)。该算法可用总速度和部分流量 公式直接实现。然而，对于逆流的情况，依赖于单个饱和度的原始部 分流必须被改进以保持质量守恒。如果相速度",从上行流《,单元中确 定，该单元a,的部分流可被表示为
a/
(51)
〃'
那么总速度可被写为
(52)
々 々 '
(53)
〃,，
(54)々w.
(55)
稳定性检查
非线性传输问题的稳定性问题在于流量函数的线性化。注意饱和
值框架关于该框架能否对于所有单元从初始饱和度估计05S°S1收敛 到满足如下收敛判据的最终饱和度Xv进行检查lim，，) = F0, (56)
其中F。是右矢量。在方程(56)中，饱和度用矢量符号表示。随着 F(S)在5"附近线性化，
F0r)+f芸).(S^-SO-F。. (57)
、跟人
可以获得迭代方程
f+ 尸V(芸y. (58) 此外，在每次迭代之后，对于所有7'强制0SS}V+1S1。
在有限差分模拟的三维模型中，重力参数G依赖于流的方向
(即，对于水平流Fh， Cg=0，而对于垂直流Fv， Cg邦)。作为结果， 必须对于水平流Fh和垂直流Fv考虑两条部分流量曲线。对于附Cg〉1， 垂直流Fv可能变成非单调曲线。较轻的流体饱和度(5*。)的方程被求解， 而另一相由饱和度限制(Sw-l-51。)进行计算。 重力的考量
重力的其它考量的影响有可能改变部分流量曲线F-yr》，方程
(31)和(32)，以使得这些曲线可能包含反射点&*。如上所述，Cg，方 程(30)被计算，并在最终估计该部分流量曲线/Y"时在方程(24)和(25) 中使用。
在没有逆流的情况中，部分流量曲线将仍像第I章所述没有任何 反射点S/。然而，如果存在逆流，那么诸如图19A-B所示的曲线可 能包含一个反射点&*。
区域I和II可被定义为由反射点5"/分开的区域，区域I是在拐 点之下0<S <5/的区域，而区域II是在S <1.0之上的区域。
在正在估计^的网格单元靠近无流边界的情况下，上述第I章 的方法被使用，而不考虑任何反射点S/。这是因为逆流不可能在这 种条件下发生。可替换地，如果在相邻单元i和n中的y'位于不同区 域I和II，那么上述第I章的方法被使用以允许饱和度^"穿过反射 点&*。在网格单元i是(l)不靠近无流边界，或者(2)单元i和n位于相同 区域的情况下，那么^"对于5"所在的特定区域被确定。例如，如果 5"v在区域I中,那么^T"将继续位于区域I中。如果/("/ 0， 那么无需低松弛来确保稳定性。尽管如此，在5""位于区域I外侧的 情况下，5^将被设置成0或者&*。如果/("/'(f )<0，即，在该区 域中拐点se的两侧，那么低松弛将被使用以确保稳定性。再次，如果 ^"位于5"所在的区域I的外侧，在计算5""的更新或低松弛值，即 0.5(5"+^")之前，5"v"将被设置成为下限或上限，0或S/中的一个。
类似地，在区域II， 5"和^+1被检查以查看5"和5""是否位于 区域II中拐点^的同一侧。如果是的话，无需松弛来保持稳性定—— 只需要进行检查以确保5"+1在区域11的范围内。如果^和^+1位于 该区域中拐点Se的两側，即/'(^)/，+')<0，那么低松弛被使用，即 5"+1=0.5(^+^+1)。再次，在进行后续的低松弛计算之前5""必须落在 区域II内。
毛细管压力可以很容易的并入到上述步骤中。《被简单地设置成 C「C。，就像下面在方程(65)中描述的那样。
在该点上，已对于这个特定的网格单元i确定了本次迭代v+l的 5"+1。如第I章所述，收敛检查在模型或子模型的每个网格单元上进 行以查看收敛是否在所有网格单元中均已达到以使得一个时步可被 认为已完成，并且在该时步无需进一步的Newton迭代。
在网格单元中更新^->^+1的算法
下面的算法将用于具有重力的双相流。
1. 检查相速度和重力的方向。
2. 如果相速度是并向的，使用常规的稳定性分析；
*计算定向的部分流量及其二阶导数(/V和); * 计算对于垂直流的Cg (对于水平流，Cg = 0);并且 *使用饱和度限制0ST+1￡1 。
參 &口果f;'Gr+')《'(^)〈o或《(sv+')《cr)<o ，对于饱和度迭4戈4吏用 低松弛，即5"+1=0.5(5"+5"+1)。
25如果相速度是对向的；
* 计算反射饱和度5V、并且
*使用饱和度限制如果该单元不靠近边界或者相邻单元并不 属于不同的饱和度区域，那么0￡5"+1^/或5;*<^+1<1，否则0SS"+Sl。
* 如果《or+')《or)〈o或《or+')《or;KO ，对于饱和度迭^U吏用
低松弛，即^+1=0.5(5"+5"+1)。 毛细管压力
在多孔介质内的多相流中，由于对岩石的润湿性能和在流体之间 的表面张力，流体趋向于具有不同的压力。这种毛细管压力经常用饱 和度的简单函数来模拟
(59)
假设相l是润湿的，而相2是非润湿的。通常，毛细管压力控制 小孔中微小尺寸流，但是它在储藏模拟尺寸中扮演很小的角色。因此， 不希望毛细管压力改变传输方程的稳定性计算的任何特征。
具有毛细管压力的Darcy判据可被写作
M- —(\7Pl - p,gVZ)) (60)
M ~(Vp ^gV2 ) (61) 从方程(59)-(61)，<formula>formula see original document page 26</formula> (63) 其中<formula>formula see original document page 26</formula>
传输方程不会因为毛细管压力的出现改变其特性。毛细管压力的 影响可通过改进重力参数G被包含。对于总速度和重力具有不同方向的情况(即，cpo的相i)，毛细管压力减小重力的影响。毛细管
压力包含毛细管压力的空间梯度()，并且希望该项相对粘度项//
或重力项cv保持更小。如果存在任何毛细管压力，相同的稳定性算 法可通过改进G来使用
^-C广C。. (65)
III.对于具有重力和毛细管压力的三相流的算法 现在将描述对于具有重力和毛细管压力的三相流的无条件稳定 算法。在储藏模拟中，三相流区域在模型中非常有限，并且通常模型 中出现的三相区域是由两个双相流的结合引起的。结果，并不清楚在 实验室的三相相对渗透率测量可被直接应用于储藏模拟。此外，三相 相对渗透率的测量在技术上非常困难，M丄Oak, L.E. Barker, and D.C. Thomas, T/rree—/;/r"5^ re/ //ve/ e/7we 6///(v o/Bre" saw^s/Wwe, J. Pet. Tech., pages 1057-1061， 1990。对于在实验测量中给出的不确定性 以及在储藏模型中网格块的相对渗透率的尺度问题，对于稳定性分析 的相对渗透率的简单幂律模型被使用。对于三相流的部分流量曲线的 一般特征将被检查。根据部分流量曲线的特征，无条件稳定的算法对 于具有重力和毛细管压力的三相流进行设计。 三相流
由于物理环境的复杂性以及在测量三相相对渗透率中的实际困 难，多孔介质中的三相流并没有被很好地理解。然而，在储藏模拟中， 网格单元中的三相流通常被发现主要通过双相流的合并而出现。尽管 如此，这使得储藏模拟中的三相相对渗透率可以是水油或油气系统的 两个相对渗透率的平均值。
改进的Stone II方法常常在许多模拟器中被使用，其中油的相对 渗透率用复杂的函数形式模拟
<formula>formula see original document page 27</formula>
这里，A^。， A 和/^g分别是油，水和气的相对渗透率。A^。ew是在 自然的水的饱和度下的油的相对渗透率，&。w是油水系统中油的相对渗透率，而A^。g是油气系统中油的相对渗透率。水和气的相对渗透率 与双相系统给出的相同。改进Stone II方法的油的相对渗透率被发现 对于一些系统变成负的。不足的实验证据不允许对于三相相对渗透率 使用诸如Stone II的复杂非线性函数形式。A. Heiba， Three- phase relative permeability, CYXF7 C 7>cAw/c / Afe附wawJw附，1987， 提出了 一种A^。w和A^。g的直线内插。对应实验测量中给出的不确定性和储藏 模拟中网格块的相对渗透率的尺度问题，对于稳定性分析的简单幂律
模型被使用。该分析可对于任意合理的三相相对渗透率的模型(无负 值或非单调曲线)被进一步扩展。
具有重力的三相流将首先被检查。与上面在双相流稳定性分析中 的一样，毛细管压力将通过改进重力数值在稍后被并入。 对于三相流的Darcy判据可被写作
",-^dA，). (67)
这里，A是渗透率，《是相速度，&是相对渗透率，a是粘度，D 表示深度，而/>是密度。下标/表示相特性(即，i-l，2，3分别表示水 相，油相和气相)。从方程(67)，可以容易地推出
",-肌S2，&)"-肌S2，S》Cg,-肌S2,S3)C《2, (68)
w2-/2(S1，S2，JS> + /4(1S1，152，S3)Cgl-/6(>51，S2,S3)C^，and (69)
w3 -^(Sp^S^+^S'^AX^ —/6(S,，S2，S3)Cs3. (70)
总速度被定义为
w - w〗+ w2 + w3
而部分参数由下式给出
_^2、2_ (73)
28
(71)
(72)/3
W人3
(74)
/4
fcrl + m2ftr2 +附人3
(75)
/5
(76)
一 ，3^f3 人--
(77)
其中，附i是关于相1的粘度的粘度比率a/m。由重力引起的速度 用变量Cgi表示
肌
(79)
Cg3 = Cg2 — Cg，
(80)
其中，"c是特征速度。注意(68)-(70)中的速度"和"i关于"c是 无量纲的。由于存在饱和度的限制
l-5*,+& + &, (81)
在(68)-(70)中，只有两个方程是独立的。方程"！和"2被选作独 立方程，同样通过用l-5V&代替513， ^和5^也被选作独立变量。 从方程(68)和(69)检查部分流表面是有意义的
/T — "I 户l 一 —
(82)
(83)用典型的流度比率w2 = 0.2和m3 = 10以及二次相对渗透率选择一个
示例，
& =《. W
在图20和21中，部分流表面对于无重力(Cw= Cg产0)和强重 力(Cg尸5并且Cg产lO )分别进行绘制。
部分流量曲线对于小的Q,和Cg2是单调的，而对于大的Cg7和
Cg2，它可以变成多值函数。注意依赖于重力和速度的方向，部分流量 可以小于零或大于一。如果Newton非线性迭代不考虑这些表面的实
际模式，这将引起数值上的困难。
从(68)和(69)，部分流量变为零的边界曲线可净皮获得 0 = w-w2C^r2-/ 3CSA3:以及 (85)
0 = " + C^,,-w3Cf^3. (86)
对于二次相对渗透率，方程(85)和(86)变为
0 - w - w 2 -《2(1 -《-(87)
0 = " - 附 - - S2)2. (88) 稳定性分析
如上所述，非线性传输问题的稳定性问题在于流量函数的线性 化。数值框架将被检查以判断它们是否能从初始饱和度估计收敛，
30<formula>formula see original document page 31</formula>此外，在每次迭代之后强制0^SV+、1。
部分流表面
现在将检查部分流表面以理解它们的一般特征。在图22， 23和 24中，部分流量等值线分别对于无重力(Cg尸Cg产0)，弱重力(Cg尸 0.5并且)和强重力(Q尸5并且Cg尸lO )三种情况进行绘制。
同时对该图使用的是^ = 0.2和m3 = 10 。 如果不存在重力，每一相的部分
流在与总速度相同的方向上移动。在从0.1到0.9以0.1为间隔绘制了 等部分流量线的图22中没有负的流量区域。注意，相一和相二的小 饱和度的大面积区域具有非常弱的部分流(<0.1)。如果像在第二中 情况中一样，重力4艮弱，对于相一和相二的负流量范围变得更小，并 且部分流量的绝对幅度也变得更小。在图23和24的正部分流中，等 部分流量曲线从0到0.8以0.2为间隔进行绘制；但是在负的部分流 量区域，对于每幅图片使用不同的等高线值对于图23A使用(-0.2, -0.4， -0.6， -0.8，國l， -2),对于图23B寸吏用(-0.05， -0.1， -0.15， -0.2， -0.25，國0，3)，对于图24A^f吏用(-0.02， -0.04， -0.06， -0.08, -1， -2 )， 而对于图24B使用(-0.05， -0.1, -0.15)。同时绘制的还有w产0 (上 侧的虚线)和《2=0 (下侧的虚线)的等高线。由于在负流量范围中的 复杂等高线，方程(91)的Jacobian矩阵变成奇异的。Jacobian矩阵的 行列式"W变为零的奇异点位置也用实线进行绘制。因为普通的Newton迭代不能穿过这些奇异直线，好的无条件稳定算法应该考虑 Jacobian矩阵的奇异点。
零部分流量边界可以对于相l(F产O)和相2(F尸0)进行绘制
^队&): 0 -1—C,2《-Cg2 w3 (1—S, - S2 )2 ， 以及 (92)
￡2): 0 -1 +- c-3 3 (", - &)2. (93)
从方程(93)和(94),交点可用轴5"2=0来完成。这些点被分别表示 为来自(92)的和来自(93)的5T2。对于将在稍后描述的算法中使用 的更好的饱和度的初始估计，在这些零流量曲线上定义了两个点，OS,
Sf =2/3^' Sf2-1/3S,。2 (94) A (S尸，Sf): 0 = 1 - Cslm2 (Sf )2 - Cf 2m3(1 - 1 - )2, (95) Z^S/^Sf ):0-l + Cg'(S/")2 —CVt73(卜S尸—Sf )2. (96)
一个点(&P3 51,)在《<0,《<0并且/)">0的区域中间净皮表示 (即，图24A中的(0.3，0.3))。
用于在具有重力和毛细管压力的三相流中更新^->^+1的算法
1. 检查总速度和重力的方向；
2. 计算Cg,和C^;
3. 对于相1计算零部分流量的边界(丄,(SA))和(A(SA));
4. 分别计算仏(S, A))和(丄2 (S A))上的饱和度初始点(V, W)和
5. 检查《和《的符号；
6. 从饱和度的初始估计(5/和&"，计算部分流量估计(iV" 和尸2°)和Jacobian矩阵的4亍列式Det;
7. 如果《>。并且《》0;
(a) 如果饱和度的初始猜测不能得到正定，将在a(《'，《)的点作 为4刀始猜测；
(b) 根据方程(91)迭代饱和度；(c) 如果《、o或《、0，在W,sf)设置饱和度;
(d) 如果5^+5^1，规格化饱和度以确保&^0;
(e) 如果，1 s， <o或<o低松弛饱和度迭代，即，
、,as, 《as22 《
5"v+H5"+,1)/2.，
8. 如果《<0，《20，并且Z)论O，将在5f')上的点的初始
饱和度估计设置为初始猜测；
9. 如果《<0， F2*20，并且1)^<0，或者《<0，《<0，并且/)^<0，
将在A(《2，《2)上的点的初始饱和度估计设置为初始猜测；并且
10. 如果《<0，《<o，并且/)">0，将在^os/"，sf)上的点的初始
饱和度估计设置为初始猜测。 对于三相稳定性的可替换算法
当三相流被模拟为两个双相流区域(即，气油和油水区域)的组 合时，稳定性算法可进行如下改进。双相流区域的部分流和它们的饱 和度(即气油相和油水相)从三相相对渗透率近似中被确定(例如， 油的饱和度在两个区域中是相同的，或者油的相对渗透率在两个区域
中是相同的)。对于双相流的稳定性分析被应用于Newton-Raphson 迭代中的每个区域。如果不稳定的饱和度变化在一个或两个区域中发 生，进行低松弛处理以保证Newton-Raphson迭代的稳定收敛。 毛细管压力
在多孔介质内的多相流中，由于它们对岩石的润湿性能和在流体 之间的表面张力，流体趋向于具有不同的压力。这种毛细管压力经常 用饱和度的简单函数来模拟。在三相流中，我们有两条独立的毛细管 压力曲线
A!(SA)-P广A以及 (97) 威,S2)"广/v (98)
毛细管压力控制小孔中的微小尺寸流，但是它在储藏模拟中扮演 很小的角色。不希望毛细管压力改变传输方程的稳定性计算的任何特 征。
具有毛细管压力的Darcy判据可被写作
33<formula>formula see original document page 34</formula>
方程(97)-(101)得到
<formula>formula see original document page 34</formula>
注意传输方程不会因为毛细管压力的出现改变其特征。毛细管压 力的影响可用重力参数&包括。对于总速度和重力具有不同方向的 情况(即，C,0的相1)，毛细管压力减小重力的影响。因为毛细管 压力包含毛细管压力的空间梯度(VPJ ，并且希望该项相对粘度项p
或重力项Cg保持更小。如果存在任何毛细管压力，第II章的相同算 法可通过改进Cg来^f吏用
<formula>formula see original document page 34</formula>
权利要求
1.一种求解非线性S形函数F＝f(S)的方法，该函数表示物理系统中的特性S，所述S形函数具有拐点Sc以及f(S)收敛到的预定值F，所述方法包括以下步骤a)选择初始值Sv+1作为函数f(S)的可能解；b)在Sv处对函数f(S)执行Newton迭代(T)以确定下一迭代值Sv+1；c)判断Sv+1是否相对Sv位于拐点Sc的另一侧；d)如果Sv+1相对Sv位于f(S)的拐点的另一侧，设置Sv+1＝Sl，一个改进的新估计；e)判断Sv+1是否在预定的收敛判据内；f)如果Sv+1不在预定的收敛判据内，设置Sv＝Sv+1；g)重复步骤b-f，直到Sv+1在预定的收敛判据内；以及h)选择Sv+1作为S形函数F＝f(S)的令人满意的解S。
2. 权利要求1的方法，其中改进的新估计5"被设置为^ = 0^+(1國01)5"+1，其中0<a<l。
3. 权利要求2的方法，其中a=0.5。其中0.3<a<0.7。其中稳定性极限^被设置为函数f(S)的
4. 权利要求2的方法，
5. 权利要求l的方法, 拐点Se。
6. 权利要求l的方法, 用的部分流量函数。
7. 权利要求l的方法， 数学表达式其中S形函数f(S)是在储藏模拟器中使 其中S形函数f(S)在数学上符合下面的<formula>formula see original document page 2</formula>其中S=双流体流量问题中第一相的饱和度；f(S)是饱和度S的函数；第一相中的流体粘度；而 ^2=第二相中的流体粘度。
8，权利要求l的方法，其中，如果o^+'-^|,则5""在预定的收敛判据内，这里^是预定的容许极限。
9. 权利要求l的方法，其中，如果o^-/(圹')|，则5"+1在预定的收敛判据内，这里s是预定的容许极限。
10. 权利要求l的方法，还包括将^"设置为上或下界Sb,如 果在步骤b中5"+1超过这些界限的话。
11. 权利要求1的方法，其中所述S形函数F = f(S)由查找表和 普通解析函数之一来表征。
12. —种机器可读的确实包含指令程序的程序存储设备，所述指 令可被所述机器执行以执行求解非线性S形函数F = f(S)的方法步骤， 该函数表示物理系统中的特性S，所述S形函数具有拐点SC以及f(S) 收敛到的预定值F，该方法包括以下步骤a) 选择初始值^作为函数f(S)的可能解；b) 在5"处对函数f(S)执行Newton迭代(T)以确定下一迭代值c) 判断sv"是否相对5"位于拐点y的另一侧；d) 如果^"相对y位于f(S)的拐点的另一侧，设置^+1=5*'， 一 个改进的新估计；e) 判断^+1是否在预定的收敛判据内；f) 如果5"+1不在预定的收敛判据内，设置5" = 5"+1;g) 重复步骤b-f，直到5""在预定的收敛判据内；以及h) 选择^+1作为S形函数F-f(S)的令人满意的解。
13. —种在储藏模拟器的网格单元中在与具有重力的双相流相 关联的部分流量函数/(A上更新Newton迭^的方法，该方法包括(a)检查网格单元中相速度和重力的方向(i)如果相速度是并流，应用常规的Newton稳定性分析；(ii) 计算定向的部分流量Fh和Fv，以及它们的二阶导数Fh"和/V';(iii) 计算垂直流的Cg,对于水平流，Cg = 0;(b) 应用饱和度限制0SSV+1￡1;(i) 3口果《'0^+1)《'(^)<0或《(^+1)《0^)<0 ，对于饱和度迭4戈应用低松弛；(c) 如果相速度是逆流； 计算反射饱和度S/;(i) 应用饱和度限制如果所述单元不邻近边界或者相邻单元 并不属于不同饱和度域，那么0SS""sy/或5V、5""〈1，否则0ST+、1; 并且(ii) :i口果《or+')/Vor) <o或《'or+"《or)<0 ，对于饱和度迭4戈应 用低松弛。
14. 一种在储藏模拟器的网格单元中在与具有重力的三相流相 关联的部分流量函数/(A上更新Newton迭代的方法，该方法包括(a) 检查总速度和重力的方向；(b) 计算C^/和C^;(c) 对于相1计算零部分流量的边界1^(S!， S。和L2(S^ S2);(d) 分别计算L"S^ S2)和"(Sh S2)上的饱和度初始点(W,《)和(e) 检查《和《的符号；(f) 根据饱和度的初始估计S/和5^，计算部分流量估计F/和 /V以及Jacobian矩阵的行列式；(g) 如果《^o并且F;20;(h) 如果饱和度的初始猜测没有得到正定，将在a(Y',W)的点作 为初始猜测；(i) 根据方程(91)迭代饱和度；(j)如果^<0或/^<0，在L,c,o设置饱和度；(k)规格化饱和度以确保5^0;(1)如果^V^〈o或^V^〈0，低松弛饱和度迭代； (m)如果《、0， F;》0，并且Z)论O，将在A(《，Sf)上的点上的初始饱和度估计设置为初始猜测；(n)如果《<0， f"o，并且""<0，或者《<0， f;<o，并且Z)"〈0，将在^(^2,sf)上的点上的初始饱和度估计设置为初始猜测；以及(0)如果《<0， F2*<0，并且/>^>0，将在L,(《3，Sf)上的点上的初 始饱和度估计设置为初始猜测。
全文摘要
提供了一种用于求解非线性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^I，一个改进的新的估计。这个改进的新的估计S^I最好被设置成拐点S^C，或者S^V和S^V+1之间的平均值，即，S^I＝0.5(S^V+S^V+1)。重复上述步骤直到S^V+1在预定的收敛判据内。此外，描述了对于具有重力和毛细管压力的双相和三相流的求解算法。
文档编号G06G7/48GK101657824SQ200680014944
公开日2010年2月24日 申请日期2006年3月15日 优先权日2005年3月15日
发明者哈姆蒂·A·切莱皮, 星·H·李, 詹妮·帕特里克 申请人:切夫里昂美国公司