用于多孔介质中的多相流的多尺度方法

具体实施方式

本文所述的本发明的实施例一般涉及有效确定地下储层中的多相流引起的精细尺度饱和度、尤其是用在储层模拟器中的多尺度方法。如本文将更详细描述的，利用饱和度的输运方程的自适应粗尺度算子克服支配方程的双曲线特征，即，饱和度前沿与基础异构渗透率分布的复杂非线性相互作用。

支配方程和离散化公式

在体积Ω上，像地下地层中的油和水那样的异构域中的二相不可压缩流在数学上可以通过下式表示：

$\Phi \frac{{\partial S}_o}{\partial t}-\frac{\partial }{{\partial x}_i}\left(k\frac{k_{r_o}}{{\mu }_o}\frac{\partial p}{{\partial x}_i}\right)=-q_o$(方程1)

$\Phi \frac{{\partial S}_w}{\partial t}-\frac{\partial }{{\partial x}_i}\left(k\frac{k_{r_w}}{{\mu }_w}\frac{\partial p}{{\partial x}_i}\right)=-q_w$(方程2)

其中，p是压强，S_j是饱和度(下标j代表相；o代表油，w代表水)，以及0≤S_j≤1和S_o+S_w＝1，k是异构渗透率，是相对渗透率(是S_j的函数)，μ_j是粘度，q_j是代表井的源项。符号Φ表示孔隙率，t表示时间。记号S＝S_o将用在下文中。该系统假设毛细压强和重力忽略不计。在Ω上，可以将方程(1)和(2)等效地改写成：

$▿·\lambda ▿p=q_o+q_w$(方程3)

$\Phi \frac{\partial S}{\partial t}+▿·\left(f_{o}u\right)=-q_o$(方程4)

并且总速度变成：

$u=-\lambda ▿p$(方程5)

以及总迁移率和油相分数流(fractional flow)为：

λ＝k(k_o十k_w)    (方程6)

$f_o=\frac{k_o}{k_o+k_w}$(方程7)

这里，对于j∈{o，w}，渗透率异构性是决定天然多孔地层中的流动行为的最主要因素。渗透率k的异构性通常被表示成空间的复杂多尺度函数。此外，因为k趋向于高度不连续，解析空间相关结构和捕获渗透率的可变性需要高度详细的描述。

在边界上，规定通量u·n，其中n是指向外面的边界单位法线向量。方程(3)和(4)是通常通过地下流动储层模拟器有效处理的系统类型的代表性描述。注意，处理不可压缩流的极限情况的能力保证了也可以求解可压缩系统，因为可压缩性使支配方程的系统刚性更低。

对于单元i，方程(3)的离散形式变成：

$\underset{j\in {\aleph }_i}{\Sigma }$(方程8)

这里，表示单元i的相邻单元，λ_ij和t_ij分别表示单元i和j的界面上的迁移率和透过率。离散化输运方程变成：

$\frac{{\mathrm{\Phi V}}_i}{\mathrm{\Delta t}}(S_i^{n+1}-S_i^n)=\underset{j\in {\aleph }_i}{\Sigma }{f}_{o,\mathrm{ij}}{u}_{\mathrm{ij}}=-q_{o,i}$(方程9)

其中，单元i和j的界面上的离散速度通过下式给出：

u_ij＝t_ijλ_ij(p_j-p_i)    (方程10)

限制和延拓运算

原先精细尺度域可以表示成Ω^h，而粗尺度域可以表示成Ω^H，其中H＞＞h。Ω^h和Ω^H不仅分别代表具有网格尺度h和H的空间，而且代表定义在网格上的向量空间。Ω^h中的线性化压强方程和非线性饱和度方程可以写成：

L^hp^h+q^h＝0，p^h∈Ω^h    (方程11)

A^h(S^h)+r^h＝0，S^h∈Ω^h    (方程12)

压强方程是包括线性算子L的p^h的线性方程。由于用方程(12)中的非线性算子A表示的分数流对饱和度的非线性函数相关性，饱和度方程是非线性方程。非线性方程(12)可以使用牛顿法迭代求解。压强的限制和延拓算子可以定义成：

$p^H=R_h^H{p}^h$(方程13)

$p^h=I_H^h{p}^H$(方程14)

然后，可以将压强方程写成：

$R_h^H{L}^h{I}_H^h{p}^H+R_h^H{q}^h=0$(方程15)

当粗尺度网格算子通过下式定义时：

$L^H=R_h^H{L}^h{I}_H^h$(方程16)

可以按粗尺度网格格式将方程(15)写成：

L^Hp^H+q^H＝0，p^H∈Ω^hH    (方程17)

类似地，可以将饱和度方程的限制和延拓算子定义成：

$S^H=R_h^H{S}^h$(方程18)

$S^h=I_H^h{S}^H$(方程19)

然后，可以将方程(12)写成：

$R_h^H{A}^h\left(I_H^h{S}^H\right)+R_h^H{r}^h=0,$S^H∈Ω^H    (方程20)

可以用粗尺度网格算子和由A的非线性引起的τ校正将方程(20)改写成：

$A^H\left(S^H\right)+R_h^H{r}^h-{\tau }_h^H=0,$S^H∈Ω^H    (方程21)

其中，τ校正通过下式定义：

${\tau }_h^H=A^H\left(R_h^H{S}^h\right)-R_h^H{A}^h\left(S^h\right)$(方程22)

多尺度有限体积近似

压强的限制和延拓算子

在通过引用并入本文的美国专利#6823297和#7496488所教导的多尺度有限体积方法中，采用二重粗尺度网格中的基函数(basis function)来构建压强的延拓和限制算子。

参照图1，将代表地下储层的区域划分成具有精细尺度单元12的2D精细尺度网格10。具有M个单元22和N个节点24的一个相适应的初级粗尺度网格20(如粗实线所示)构建在原先精细尺度网格上。每个初级粗尺度单元22，即，(i∈{1，...，M})由多个精细尺度单元12组成。也与精细尺度网格10相适应的二重粗尺度网格30(如虚线所示)被构建为使得每个二重粗尺度单元32，即，(j∈{1，...，N})在它的内部恰好包含初级粗尺度网格20的一个节点24。一般说来，每个节点24处在一个二重粗尺度单元32的中心。二重粗尺度网格30也具有M个节点34，即，x_i(i∈{1，...，M})，每一个节点处在初级粗尺度单元22，即，的内部。一般说来，每个二重粗尺度节点34处在一个初级粗尺度单元22的中心。二重粗尺度网格30一般通过连接包含在相邻初级粗尺度单元22内的节点34构成。每个二重粗尺度单元32具有N_c个角36(在二维中四个而在三维中八个)。构建一组二重基函数每个二重粗尺度单元的每个角i一个。

为图1中的粗尺度二重网格30(Ω^D)构建数值基函数。通过求解如下椭圆问题为每个二重粗尺度单元32构建四个二重基函数Θⁱ(i＝1，4)(对于3D为8个基函数)：

$▿·(\lambda ·▿{\Theta }_j^i)=0,$在上    (方程23)

以及边界条件：

$\frac{\partial }{{\partial x}_i}(\lambda ·▿{\Theta }_j^i)=0,$在上    (方程24)

其中，x_i是与的边界相切的坐标。的节点x_k的值通过下式给出：

${\Theta }_j^i\left(x_k\right)={\delta }_{\mathrm{ik}}$(方程25)

其中，δ_ik是克罗内克三角(Kronecker delta)。对于定义一旦计算出基函数，就可以容易地将延拓算子写成：

$p_j^h\left(x\right)=I_H^h\left(p^H\right)=\underset{i}{\Sigma }{\Theta }_j^i{p}_i^H,$$x\in {\Omega }_j^D$(方程26)

其中，精细尺度网格压强而粗尺度网格压强延拓算子是二重网格基函数的线性组合。

图2对于粗尺度网格20上的粗尺度单元示出了如实线所示的八个相邻初级粗尺度或粗尺度单元，具有二重粗尺度单元(j＝A，B，C，D)、如虚线所示的二重粗尺度网格30在图2中通过阴影线示出。粗尺度网格算子可以从每个粗尺度单元的守恒方程中构建。第一步是计算作为粗尺度体积1-9中的粗尺度压强p^H的函数、位于Ω^H内、穿过粗尺度单元界面分段26的通量。这是通过每个相邻粗尺度单元一个地构建八个二重基函数Θⁱ(i＝1，8)实现的。Ω^H内的粗尺度通量可以通过叠加这些基函数，作为粗尺度压强pH的函数获得。从二重基函数Θⁱ_j中提取有效粗尺度透过率。为每个粗尺度单元Ω_i^H构建一组精细尺度基函数Θ_i^k，以便为每个相邻粗尺度单元Ω_k^H构建一个包括Ω_i^H的基函数。精细尺度基函数的局部计算的边界条件从二重基函数中提取。最终，在给出粗尺度解p_i^H的情况下，则穿过界面

${\partial \Omega }_{i_1{i}_2}^H=\partial {\Omega }_{i_1}^H\cap \partial {\Omega }_{i_2}^H$(方程27)

的相通量通过下式近似：

$F_l^{i_1{i}_2}=\underset{i=1}{\overset{M}{\Sigma }}{T}_{l,i}^{i_1{i}_2}{p}_i^H,$l＝o，w    (方程28)

其中，是通过下式给出的相透过率：

$T_{l,i}^{i_1{i}_2}={\int }_{{\partial \Omega }_{i_1{i}_2}^H}\underset{j=1}{\overset{N}{\Sigma }}({\lambda }_l^{\prime }·▿{\Theta }_j^i)·\mathrm{nd\Gamma }$(方程29)

向量n是方向从指向的的单位法线。从粗尺度网格透过率中获取粗尺度网格方程：

L^Hp^H＝Tp^H＝-q^H，(方程30)

其中，T是压强的粗尺度网格透过率算子。

饱和度的限制和延拓算子

为了在粗尺度和精细尺度网格中保持质量守恒，饱和度的限制算子的选择在很大程度上局限于体积平均：

$S_i^H=R_h^H{S}^h=\frac{1}{V_i}\underset{l\in {\Omega }_i^H}{\Sigma }{\upsilon }_l{S}_l^h$(方程31)

其中，υ_l是精细尺度单元l的体积，而V_i是粗尺度单元i的体积。如果给定精细尺度网格10中的相速度，可以将非线性精细尺度网格算子写成：

$A_i^h\left(S^h\right)=\underset{j\in {\aleph }_i}{\Sigma }\frac{k_o\left(S_{o,\mathrm{ij}}^h\right)}{k_o\left(S_{o,\mathrm{ij}}^h\right)+k_w\left(S_{w,\mathrm{ij}}^h\right)}{u}_{\mathrm{ij}}-{\beta }_i(S_i^{h,n+1}-S_i^{h,n})$(方程32)

这里，

${\beta }_i=\frac{{\mathrm{\Phi V}}_i}{\mathrm{\Delta t}}$(方程33)

S_ij和u_ij分别表示单元i和j之间的逆风饱和度和总相速度。而且，将粗尺度网格总速度和分数流定义成：

$H_{\mathrm{ij}}^H=\underset{l\in {\partial \Omega }_{\text{ij}}^H}{\Sigma }{u}_l^h$(方程34)

$F_{\mathrm{ij}}^H=\frac{1}{U_{\mathrm{ij}}^H}\underset{l\in {\partial \Omega }_{\mathrm{ij}}^H}{\Sigma }f\left(S_l^h\right)u_l^h$(方程35)

分数流曲线是非线性函数(即，S形)，另外，多相流复杂地与异构渗透率相互作用。粗尺度网格分数流F_ij^H一般说来是不能只通过粗尺度网格饱和度S^H的函数容易地表示的S^h的复杂非线性函数。但是，如果粗尺度网格中的饱和度单调增大或减小，那么在多相流前沿移过网格之后，可以从前一个时间步长或前一次迭代的饱和度变化中估计粗尺度网格分数流曲线。精细尺度网格中的饱和度变化可以通过下式表达成粗尺度网格饱和度变化的分数：

${\xi }_l=\frac{\delta S_l^h}{\delta S_i^H},$$l\in {\Omega }_i^H$(方程36)

可以将分数流曲线线性化成：

$f(S_l^H+{\mathrm{\delta S}}_l^h)\cong f\left(S_l^h\right)+\frac{\partial f}{{\partial S}_l^h}{\xi }_l\delta S_i^H$(方程37)

$F_{\mathrm{ij}}^H(S_i^H+{\mathrm{\delta S}}_i^h)\cong F_{\mathrm{ij}}^H\left(S_i^H\right)+\frac{{\mathrm{\delta S}}_i^H}{U_{\mathrm{ij}}^H}\underset{l\in \partial {\Omega }_{\mathrm{ij}}^H}{\Sigma }\frac{\partial f}{{\partial S}_l^h}{\xi }_l{u}_l^h$(方程38)

借助于粗尺度网格分数流的模型，现在可以将粗尺度网格非线性算子写成：

$A_i^H\left(S^H\right)=U_{\mathrm{ij}}^H{F}_{\mathrm{ij}}^H-{\beta }_i(S_i^{H,n+1}-S_i^{H,n})$(方程39)

现在可以导出将粗尺度网格饱和度内插成精细尺度网格饱和度的延拓算子。

延拓算子I：

对于还未构想可靠延拓算子的区域，可以经由施瓦兹重叠方法直接计算精细尺度网格解。当侵入流体在粗尺度网格中移动时，如巴克利-莱弗里特(Buckley-Leverett)方程中那样的陡峭饱和度前沿经历与基础异构渗透率场的复杂非线性相互作用。施瓦兹重叠延拓方法专门用在求解输运方程的原先多尺度有限体积方法中。虽然这种数值延拓是精确的，但在计算方面成本也很高。本发明的实施例自适应地将这种延拓算子应用于侵入流体首先进入并造成饱和度迅速变化的区域。对区域Ω_i^H求解方程(12)：

$A_i^h\left(S^{h,\upsilon +1}\right)+r^{h,\upsilon +1}=0,$$S^h\in {\Omega }_i^H$(方程40)

其中诺伊曼(Neumann)边界条件为：

$S_{\mathrm{ij},l}^{h,\upsilon +1}=S_{\mathrm{ij},l}^{h,\upsilon },u_{\mathrm{ij},l}^{h,\upsilon +1}=u_{\mathrm{ij},l}^{h,\upsilon },$$l\in {\partial \Omega }_i^H\cap {\partial \Omega }_j^H$(方程41)

上标υ指示迭代等级，方程(41)的近似提供问题局部化，以便可以为每个粗尺度网格求解方程。从方程(31)中可以直接计算出粗尺度网格变量。第一种方法在构建精细尺度网格饱和度时是局部守恒和精确的；但是，它可能要求利用诺伊曼边界条件在初级粗尺度网格中对饱和度进行几次迭代隐性计算，以获得真正全局收敛解。

延拓算子II：

第二延拓算子包含两个步骤：(1)经由粗尺度网格速度的直接内插重构局部守恒精细尺度网格速度；以及(2)显性计算精细尺度网格饱和度。假设粗尺度网格速度和精细尺度网格速度分布从前一个时间步长或前一次迭代(υ)中给出；粗尺度网格i为U_i^H，υ，而精细尺度网格j∈Ω_i^H为u_j^h，υ。从方程(30)的粗尺度网格解中可以获得新的粗尺度网格速度U_i^H，υ+1。在速度变化不大的情况下，可以在新迭代等级上通过下式将粗尺度网格速度变化内插成精细尺度网格速度：

$u_j^{h,\upsilon +1}\left(x\right)=u_j^{h,\upsilon }\left(x\right)+U_i^{H,\upsilon +1}\left(x_0\right)-U_i^{H,\upsilon }\left(x_0\right)$(方程42)

$+\frac{{x-x}_0}{|x_1-x_0|}\otimes (U_i^{H,\upsilon +1}\left(x_1\right)-U_i^{H,\upsilon }\left(x_1\right)-U_i^{H,\upsilon +1}\left(x_0\right)+U_i^{H,\upsilon }\left(x_0\right))$

其中，以及x₀和x₁分别是粗尺度网格的左下角和右上角坐标。如果u_j^h，υ和U_i^H，υ+1分别在精细尺度网格中和在粗尺度网格中是守恒的，则内插精细尺度网格速度u_j^h，υ+1在精细尺度网格中也是守恒的。一旦估计出精细尺度网格速度，就可以通过显性方法廉价地计算饱和度：

$S_j^{h,\upsilon +1}=S_j^{h,\upsilon }+\frac{\mathrm{\Delta t}}{v_j}\underset{l}{\Sigma }{u}_l^{h,\upsilon +1}f\left(S_j^{h,\upsilon }\right)$(方程43)

显性饱和度计算的稳定性受CFL数支配，可以将CFL数写成：

$\mathrm{CFL}=u_t\frac{\partial f\left(S_j^h\right)}{{\partial S}_j^h}\frac{\mathrm{\Delta t}}{\mathrm{\Delta x}}$(方程44)

为了显性计算的稳定性，CFL数应该小于1。如果将这种算法应用在的域中(即，在巴克利-莱弗里特前沿后面变稀薄)，那么时间步长大小限制将是次要的。如果存在由稳定性引起的时间步长大小限制，则可以使用多个时间步进或隐性公式：

$S_j^{h,\upsilon +1}=S_j^{h,\upsilon }+\frac{\mathrm{\Delta t}}{v_j}\underset{l}{\Sigma }{u}_l^{h,\upsilon +1}f\left(S_j^{h,\upsilon +1}\right)$(方程45)

延拓算子III

现在设计在粗尺度网格中局部守恒的饱和度的快速内插。如果饱和度分布模式在迭代或时间步长之间变化不大，可以从粗尺度网格饱和度变化中计算精细尺度网格饱和度：

$S_l^{h,\upsilon +1}=S_l^{h,\upsilon }+{\xi }_l\partial S_i^H$(方程46)

这里，假设相对饱和度变化(ξ₁)相对于前一次迭代变化不大。对于在陡峭饱和度前沿之后饱和度变化缓慢的粗尺度网格，这是一种似乎可取的近似。显然，这种内插器的精度依赖于相对于前一次迭代的不变量ξ₁的采用。如下面的数值例子所示，可以识别出可以安全地应用这种简单内插器来达到高数值效率和精度的区域。注意，上面的延拓算子不保证在精细尺度网格中饱和度守恒，但在粗尺度网格中饱和度守恒。因此，这可能在精细尺度网格饱和度中产生非守恒误差。不过，由于这种延拓算子只应用在饱和度缓慢变化，而且粗尺度饱和度总是局部守恒的区域中，所以它们可以保持很小和受限。

守恒速度场的精细尺度压强

原先多尺度有限体积方法观察到，通过粗尺度压强值p_i^H和二重基函数Θⁱ_j构建的精细尺度网格速度在粗尺度二重单元或体积之间的界面附近产生局部质量平衡误差。提出了第二组基函数来构建守恒速度场。计算精细尺度基函数时将从二重基函数中获得的精细尺度通量用作边界条件是保证重构精细尺度速度场守恒的关键。穿过粗尺度单元界面26的精细尺度通量是从二重粗尺度单元的内部提取的(参见图2)。精细尺度基函数(对于3-d为27个，而对于2-d为9个)是利用从每个节点的二重基函数中计算的边界通量求解压强方程构建的。在这种第二基函数手段中存在两个数值难题。计算量因大量第二基函数(对于3-d为27个)而变得相当大，并且基函数方法局限于可以应用叠加规则的线性系统。即使可以应用自适应方法来消除第二基函数的不必要计算，大量的基函数也保证不了这种手段在数值计算上比针对每个时间步长或每次迭代都计算速度和饱和度的直接方法更有效。而且，这种直接手段不要求支配方程严格线性，这可以容易地推广成包括方程的许多非线性效果，即，可压缩性、毛细压强、相对渗透率等。

如图2所示，根据粗尺度网格压强解和二重粗尺度网格基函数，在初级粗尺度网格θΩ^H中构建精细尺度网格压强场。沿着初级粗尺度网格θΩ^H应用诺伊曼边界条件：

$u^{\prime }={u^{\prime }}_o+{u^{\prime }}_w={\lambda }_o·▿p+{\lambda }_w·▿p$(方程47)

并且为精细尺度速度而求解方程(3)。为了获得服从u′的守恒速度场，在粗尺度单元之间的边界上施加如下要求：

u·n＝u′·n    (方程48)

其中，n是界面单元法线向量。局部精细尺度解是具有方程(48)的局部边界条件的方程(3)的解和容易转换成速度的局部压强解。注意，当将延拓算子II或III用于饱和度输运方程时，线性构建或根本不需要计算精细尺度速度。因此，在应用这些延拓算子时可以获得很高计算效率。

顺序完全隐性方案

数值稳定性的隐性多尺度有限体积算法在本领域也是已知的。每个时间步长由一个牛顿循环组成，并且分两个阶段迭代求解多相流问题。首先，在每个牛顿循环中，求解压强方程，并且从压强解中构建速度场。然后，使用所构建精细尺度速度场u在精细尺度网格上求解输运方程。在施瓦兹重叠方法中，利用隐性逆风方案在每个粗尺度单元或体积上局部求解输运问题。将来自前一个迭代等级上的相邻粗尺度体积的饱和度值用于边界条件。一旦输运方程收敛，新的饱和度分布就决定了下一次牛顿迭代的椭圆问题的总迁移率场。已经针对巨大和刚性问题测试了顺序隐性多尺度有限体积方法，并且，甚至对于很大的时间步长，耦合方案也未失败。本发明的实施例采用顺序完全隐性算法；但是，自适应地构建精细尺度特性。通过自适应地构建压强、速度和饱和度，可以不损害数值精度和稳定性地获得高数值效率。

压强的自适应方案

多相流的隐性粗尺度有限体积算法的最昂贵部分是提供压强限制和延拓算子的二重基函数的重构。因此，为了获得高效率，最好是只有在绝对必要的情况下才重新计算它们。可以利用自适应方案来更新这些基函数。如果对于二重粗尺度单元内的所有精细尺度单元，如下条件：

$\frac{1}{1+{ϵ}_{\lambda }}<\frac{{\lambda }^v}{{\lambda }^o}<1+{ϵ}_{\lambda }$(方程49)

未得到满足，那么，必须重构该控制体积的二重基函数。这里，λ⁰表示最后更新基函数时的总迁移率，λ^v表示当前迭代的总迁移率。ε_λ＞0是用户定义值。注意，如果λ变化了大于1/(1+ε_λ)并小于1+ε_λ的倍数，这个条件(方程49)为真。当然，这些数字密切依赖于用户定义阈值ε_λ。一般说来，更小的阈值触发更多的精细尺度体积，结果是，每个时间步长都要重新计算更多的基函数。对于多种测试情况，将ε_λ取成＜0.2在所得结果中造成边际变化，并且需要重构的基函数的分数(fraction)很小。在本文后面给出的所有数值例子中，将标准ε_λ设置成＜0.2，以便更新基函数。

饱和度的自适应方案

在原先多尺度有限体积方法中，通过施瓦兹重叠方法在精细尺度网格中迭代求解输运方程。通过利用诺伊曼边界条件为每个粗尺度网格求解局部守恒方程，构建局部守恒精细尺度速度场。而且，经由总速度固定的施瓦兹重叠方法迭代求解非线性饱和度。一旦经由原先多尺度有限体积方法大大地优化压强计算，这个过程就变成最耗时部分。

上面导出了饱和度的一个限制算子和三个延拓算子，使得可以不损害数值精度地调适有效算法。在驱替(displacement)过程中，注入流体从注入井移动到生产井。将建立起在饱和度前沿与渗透率的基础异构性之间具有复杂相互作用的类巴克利-莱弗里特状饱和度分布。因此，本发明的实施例通过将模型划分成三个区域而包括自适应算法。区域1被定义成注入流体未到达的地方。区域2被定义成注入前沿已经侵入并且注入流体的饱和度迅速增大的地方。区域3被定义成陡峭前沿移过并且精细尺度单元之间的饱和度分布被大部分建立的地方。饱和度变化像巴克利-莱弗里特方程的稀散波那样缓慢。

在区域1中，由于饱和度没有变化，所以无需计算精细尺度饱和度。区域2要求不需要饱和度发展的任何先前历史的最严格算法，即延拓算子I。在区域3中，可以应用近似线性内插器，即延拓算子II或III，一般说来，算子的选择基于可以容易地从先前所计算的精细尺度饱和度迭代中获得的精细尺度解变化。延拓算子II在计算方面比延拓算子III更昂贵，但在精细尺度上获得局部守恒方案，而延拓算子III只在粗尺度上是局部守恒的。因此，为了保证精细尺度下的非守恒误差是受限的和微小的，在区域3中使用算子III需要良好的标准。

需要一个在区域1，2和3之间建立转变的标准。这种标准可以基于粗尺度网格中的饱和度变化和总速度变化。例如，可以为粗尺度网格i检测从区域1到区域2的转变，并且可以通过下式表示它：

$\left|\right|{\mathrm{\Delta S}}_i^H\left|\right|>{\Delta }_1$(方程50)

其中，Δ₁大于0，一般是大约10^-5到10^-1。从区域2到区域3的转变可以通过下式饱和度和速度两者的变化来识别：

$\left|\right|{\mathrm{\Delta S}}_i^H\left|\right|<{\Delta }_2$和$\left|\right|\frac{{\mathrm{\Delta U}}_i^H}{U_i^H}\left|\right|<{\Delta }_{\upsilon }$(方程51)

其中，Δ₂一般大于大约10^-3，更典型的是大约10^-3到10^-1，Δ_υ一般大于大约10^-3，更典型的是大约10^-3到10^-1。

为了选择区域3中的算子II或III，需要另一条标准。如果方程(46)中的相对饱和度变化ξ₁波动得不厉害，算子III应该产生很小的精细尺度饱和度误差。引入计量每个粗尺度单元的ξ的变化的变量：

${\Xi }_i\equiv \frac{\max \left\{{\xi }_{l,i}\right\}}{\min \left\{{\xi }_{l,i}\right\}},$$l\in {\Omega }_i^H$(方程52)

当在区域3中满足如下条件时：

ΔΞ_i＜Δ_ξ    (方程53)

可以应用算子III。取决于储层模型如何异构，Δ_ξ的典型范围从大约10^-2到10^-5。

由于这些标准极具限制性，输运方程将更多地利用延拓算法I，使其更类似于原先多尺度有限体积算法。虽然计算效率提高了，但数值误差可能随标准变宽松而增大。本发明的其它实施例应用其它区域数量，例如，2个或大于3个(例如，4，5，6个等)来标识使用不同自适应延拓算子构建精细尺度饱和度的输运区域。另外，注意尽管本文所示确定区域转变的标准不是通用的，基于粗尺度特性(即，对于区域1到2，饱和度变化，而对于区域2到3，饱和度和总速度变化)而变化，但在一些实施例中可以应用通用转变标准。

图3是压强和饱和度的延拓和限制算子的示意图。它清楚指明了压强和饱和度计算的顺序算法。该算法可以参照支配方程描述如下：

(1)构建与精细尺度网格相适应的初级粗尺度网格和二重粗尺度网格，使得初级粗尺度单元和二重粗尺度单元是精细尺度单元的聚集体。

(2)从方程(23)中导出基函数(Θ)，并且从方程(24)和方程(25)中导出局部边界条件。

(3)使用基函数(Θ)从方程(29)中计算粗尺度透过率。

(4)使用方程(28)从通量中构建粗尺度网格守恒方程，并且从方程(30)中计算粗尺度网格压强p^H。

(5)从方程(34)中计算粗尺度网格速度，并且从方程(39)中计算粗尺度网格饱和度。

(6)根据粗尺度网格速度和饱和度变化，使用方程(50)和方程(51)的标准识别区域1，2和3。并且，应用方程(53)的标准来决定区域3中的内插方法。

(7)对于区域1，注入流体的饱和度变化被认为是零。

(8)对于区域2，与基函数(Θ)一起，从粗尺度压强场p^H中构建方程(26)中的精细尺度压强p^h。从精细尺度压强解p^h(x)中构建初级粗尺度网格的诺伊曼边界条件。这个压强解将产生守恒精细尺度相速度场u_o，u_w，u_g。然后，将相速度用于使用方程(40)为精细尺度饱和度求解。

(9)对于ΔΞ_i≥Δ_ξ的区域3，使用方程(42)直接从粗尺度速度内插精细尺度速度。使用精细尺度速度直接计算精细尺度饱和度。这里，可以使用显性(方程43)或隐性(方程45)离散化方案来导出精细尺度饱和度。

(10)对于ΔΞ_i＜Δ_ξ的区域3，其中粗尺度网格中的饱和度变化处在线性渐近域中，如果有必要，使用方程(46)从粗尺度网格饱和度线性内插精细尺度网格饱和度。

(11)使用新饱和度分布更新迁移率场λ，并且在需要的情况下重新计算基函数(包括更新有效粗尺度透过率)。这里，应用自适应方案。

(12)如果应用隐性解决方法，通过重复步骤2-9继续进行下一次牛顿拉夫森(Newton Raphson)迭代，直到收敛为止。

(13)通过重复步骤2-9执行下一个时间步长。

图4-6例示性地将上面的多尺度算法浓缩成流图。方法(100)包括为地下储层提供模拟模型，该模拟模型包括定义多个精细尺度单元的精细尺度网格、和定义多个作为精细尺度单元的聚集体的粗尺度单元的粗尺度网格(步骤110)。响应于预定饱和特性定义与每个粗尺度单元相对应的饱和区(步骤120)。例如，粗尺度网格速度和饱和度变化可以用于通过使用方程(50)和方程(51)识别饱和区。为每个饱和区确定精细尺度饱和度(步骤130)。例如，对于一些饱和区，可以将精细尺度饱和度指定成零值，或对于一些饱和区，可以例如通过使用方程(40)、(43)、(45)和(46)的计算来计算精细尺度饱和度。汇集饱和区以便获取多个精细尺度单元上的精细尺度饱和度分布(步骤140)，饱和区的汇集要求在粗尺度网格上的每个粗尺度单元中组合精细尺度饱和度。这个过程可以一直迭代进行到精细尺度饱和度场收敛到特定容限。

在一些实施例中，像在图5中所强调的那样，响应于预定饱和特性为粗尺度单元定义三种饱和区(步骤120′)。第一区被定义成注入流体未侵入粗尺度单元的地方。第二区被定义成注入流体已侵入粗尺度单元并且注入流体的饱和度和总速度大于或等于预定数量的地方。第三区被定义成注入流体已侵入粗尺度单元并且注入流体的饱和度和总速度小于预定数量的地方。另外，可以响应于精细尺度饱和度分布输出视觉显示(步骤150)。视觉显示的例子可以包括模拟压强分布、饱和度分布、饱和度前沿的特性、和累计采油和油分数的表示。

图6强调了为三个识别的饱和区(步骤120″)确定精细尺度饱和度(步骤130′)的本发明一个实施例。为了识别三个饱和区，首先例如通过预定数量确定注入流体是否在粗尺度上侵入(步骤122)。如果注入流体未侵入饱和区，则将注入流体的精细尺度饱和度指定成零(步骤132)。如果注入流体侵入了饱和区，则确定与最后时间步长或迭代相比，注入流体的饱和度变化和总速度变化是否小于预定数量，这在数学上可以描述成和(步骤124)。如果注入流体的饱和度变化和/或总速度变化不小于预定数量，则通过使用诺伊曼边界条件，经由施瓦兹重叠方法求出注入流体的精细尺度饱和度(步骤134)，以便导出精细尺度相速度，然后可以应用精细尺度相速度来求解输运方程。对于这在数学上可以表示成如果注入流体侵入了饱和区，但饱和度变化和总速度变化小于预定数量，则使用延拓算子内插注入流体的精细尺度饱和度。首先确定注入流体的相对饱和度变化是否小于预定数量(即，ΔΞ_i＜Δ_ξ)(步骤126)。如果相对饱和度变化大于预定数量，则从粗尺度速度内插精细尺度速度，然后可以直接计算精细尺度饱和度(步骤136)。这里，可以应用显性或隐性离散化方案。如果相对饱和度变化小于预定数量，则可以直接从粗尺度饱和度线性内插精细尺度饱和度，这在数学上可以表示成(步骤138)。一旦构建了所有饱和度区，就将饱和度汇集在一起，以便获取多个精细尺度单元上的精细尺度饱和度分布(步骤140)。并且，这个过程可以一直迭代进行到精细尺度饱和度场收敛到特定容限。可以响应于模拟流体流动输出视觉显示。例如，视觉显示可以包括模拟压强分布、饱和度分布、饱和度前沿的特性、和累计采油和油分数量，但不局限于这些。

例子

采用带有各种边界条件(源/宿或狄利克雷(Dirichlet)边界条件)的储层模型来广泛地测试用于压强和饱和度计算的本发明多尺度方法。在同构储层、具有短各向同性渗透率相关长度的异构储层模型、和具有长各向异性渗透率相关长度的高度异构介质中研究二相流。前两个例子在模型的两个对角上包括一个注入井和一个生产井，最后一个例子牵涉到在入口和出口边界上具有恒压条件的线性驱替过程。

假设流体是不可压缩的，并应用二次相对渗透率模型(和)。水与油之间的粘度比是1∶5(不利驱替)。压强的非线性收敛容限(ε_p)被设置成1psi，饱和度的非线性收敛容限(ε_s)被设置成10^-4。在自适应输运算法中，为从区域1(在注入流体侵入之前)到区域2(由注入流体的初始侵入引起的陡峭饱和度变化)的转变选择转变标准Δ₁＝10^-5。从区域2到区域3使用转变标准Δ₂＝10^-2和Δ_υ＝0.1(饱和度前沿移过区域之后的缓慢饱和度变化)。在不同储层模型中使用相对饱和度变化(Δ_ξ)的稍微不同的阈值。已经发现，高度异构模型需要Δ_ξ的限制性容限，以便保持数值精度。在本文提供的例子中，对于前两个例子，应用Δ_ξ＝10^-3，而对于第三个例子，应用Δ_ξ＝10^-4。

精细尺度解是参考解，压强和饱和度误差的L₂范数通过下式定义：

$e_p=\frac{{\left|\right|p^{m,s}-p^f\left|\right|}_2}{\left|\right|p^{\mathrm{init}}\left|\right|}$(方程52)

e_s＝||S^ms-S^f||₂    (方程53)

在例3的线性驱替过程中，可以将入口和出口边缘之间的压强差作为特征压强用于归一化e_p。

例1：同构介质

考虑同构渗透率为k＝100md的700ft×700ft二维储层模型。尽管该模型是二维的，但为了便于描述工作条件起见，假设该模型在第三方向具有1ft的单位厚度。将精细尺度网格70×70均匀粗化成10×10粗尺度网格。由于每个粗尺度块包含7×7个精细尺度单元，所以放大因子是49。储层原来是充满油的，将水注入以替代油。水从左上角注入，而产品从右下角出现。初始储层压强是2000psi。水的注入速率在储层条件下是恒定的(50桶/天)，储层流体以相同速率产出。在粗尺度网格中均衡分配注入和生产率(例如，在左上角粗尺度网格中注入，而在右下角粗尺度网格中产出)。

图7和8分别例示了参考解决方案(A)、未使用输运计算自适应的原先多尺度有限体积解决方案(B)、和使用自适应输运计算的多尺度手段解决方案(C)在t＝0.8PVI(注入的孔隙体积)上的压强和饱和度分布的比较。在这些图中，FM表示精细尺度方法，MSOM表示在输运计算上不具有自适应的多尺度施瓦兹重叠方法，MSAT表示使用自适应输运计算的多尺度手段。通过图7和8中的三种方法计算的压强和饱和度的差别几乎是察觉不出的。

图9A-C显示了针对三个时间步长，通过MSAT自适应计算的饱和度。在图9中，白色正方形指示在该时间步长上在压强和饱和度计算的迭代中至少一次应用延拓算子I(施瓦兹重叠)的区域2。由于输运方程对于饱和度来说是非线性的，所以该算法需要多次牛顿迭代来求解它。甚至在同一时间步长中，如果饱和度变化小于转变标准，在迭代期间也可以将延拓算子从I切换成II。图9中的黑色正方形指示在该时间步长上至少一次使用延拓算子II的区域3。注意，延拓算子I主要应用在陡峭饱和度前沿周围的区域中。而且，如果饱和度和总速度变化很小，则应用延拓算子II。在以后的时间里，由于在模型的大部分中已经良好地建立起饱和度分布，所以将延拓算子III广泛应用在饱和度计算中。

表1

在表1中，为各种时间步长列出了相对于参考解决方案(FM)的MSOM和MSAT的L₂范数，并且还包括了压强和饱和度计算的自适应比。e_p和e_s分别是压强和饱和度的L₂误差范数。f_p表示更新基函数的分数(fraction)，f_I表示需要通过延拓算子I进行精细尺度输运计算的粗尺度块的分数，f_II表示通过延拓算子II采用精细尺度输运计算的粗尺度块的分数。所有自适应统计通过从第一时间步长到当前时间步长计算的平均分数来度量。例如，在计算从t＝0到t＝0.8PVI的饱和度中应用延拓算子I的粗尺度网格的分数是6.01％。

图10例示了这个例子的生产过程中的累计采油和油分数曲线。在图10中，累计采油和油分数的数值误差几乎是注意不到的。

从这个数值例子中，首先注意到，来自自适应输运计算(MSAT)的压强与来自未使用输运计算自适应的原先多尺度有限体积(MSOM)的压强一样精确。(MSAT)方法在饱和度计算方面产生比MSOM稍高的数值误差，但仍然相当小。其次，用于压强计算的基函数更新随着前沿从注入井移动到生产井不断减少。由于总迁移率变化一般小于饱和度变化，驱替过程中的压强修正相当适中。其结果是，需要更新小百分比的基函数(例如，在1.5PVI下，1.95％)。类似地，总速度也极缓慢地变化，在1.5PVI模拟期间，只要求3.44％的速度更新。相比之下，饱和度前沿在从注入井移动到生产井时经历了宽阔的转变区。需要原先精细尺度输运计算的粗尺度网格模型的分数在4.55％到26.00％之间变化。通过放宽转变标准(Δ1和Δ2)可以增大自适应性，但是，这可能产生较不精确的结果。

例2：具有对数正态渗透率的异构介质

图11描绘了用在第二个例子中的渗透率场。储层模型具有相关长度适中的异构渗透率场。除了渗透率场按以毫达西为单位的对数渗透率的均值为4和方差为2的对数正态分布之外，应用与前一种情况相同的储层模块，并且给出空间相关长度为0.2。渗透率通过顺序高斯模拟方法生成。

图12和13分别比较了三种不同方法FM、MSOM和MSAT在t＝0.8PVI上的压强和饱和度分布。在存在高度变化和相关渗透率场的情况下，水饱和度分布在图13中呈现与前例的图8中的简单对称饱和度分布形成对比的复杂结构。来自MSOM和MSAT的多尺度方法的数值结果的精度较高。最大饱和度误差处在渗透率极低的区域中。

图14A-C表明饱和度计算的自适应性与同构情况下的自适应性类似。延拓算子I主要应用在饱和度前沿附近，然后使用延拓算子II，直到建立起饱和度分布。本发明的自适应饱和度计算(MSAT)得出与未用输运自适应计算的原先多尺度方法(MSOM)同样精确的数值结果。在表2中，为各种时间步长列出了相对于参考解决方案(FM)的MSOM和MSAT的L₂范数，并且还给出了压强和饱和度计算的自适应比。

表2

图15画出了生产过程中的累计采油和油分数。多尺度方法(MSOM和MSAT)下的生产历史与精细尺度参考方法(FM)下的生产历史很难区分。这种具有对数正态渗透率的数值例子清楚地证明了自适应输运算法在再现原来多尺度结果以及参考精细尺度网格模拟结果时是高度有效的和精确的。

例3：利用狄利克雷边界条件的异构模型

图16示出了第三个例子中的对数渗透率场。该渗透率场与220×54精细尺度网格和产生11×9放大因子的20×6粗尺度网格相关联。初始压强是4000psi。左边界保持在注水的4000psi恒压上，而右边界保持在生产储层流体的1000psi恒压上。由于生产和注入速率随时间连续变化，所以特征时间可以通过下式表示：

${\tau }_o=\frac{\phi L_x^2\bar{\mu }}{\bar{k}|p^{\mathrm{left}}-p^{\mathrm{right}}|}$(方程54)

其中，和分别表示特征粘度和渗透率，L_x是x方向的模型维度。将水和油粘度的算术平均表示成而将渗透率的体积平均表示成

图17A-C和18A-C分别比较了参考解决方案(FM)和多尺寸方法(MSOM和MSAT)在t＝4τ₀上的压强和饱和度分布。尽管该模型包含高度异构渗透率场，但由于狄利克雷边界条件和沿着流动方向大得多的几何尺寸，所以图17-C中的压强场呈现沿着流动方向几乎一维的轮廓。相反，如图18A-C所示，饱和度分布清楚地指示导致极复杂、宽阔饱和度前沿的复杂异构分布。

图19A-C示出了例3中三个时间步长的自适应饱和度计算。同样，白色和黑色正方形表示在这个时间步长上在迭代饱和度计算中分别至少一次应用延拓算子I和II的区域。延拓算子I只用在饱和度前沿周围或总速度显著变化的几个区域中。在表3中更清楚地表明了这些统计。在所有列出的时间步长中应用算子I的分数小于5％；而在最后时间里，在整个区域中几乎都应用算子III。

表3

图20画出了生产过程中的累计采油和油分数。尽管高度异构和相关长度很大，但精细尺度结果和多尺度结果极其一致。多尺度自适应输运计算和原先多尺度方法的结果基本相同。由于高异构性，在本例中可以看到精细尺度结果与多尺度结果之间的细小差异，但多尺度模拟中的误差在工程应用上是可接受的。在表3中，也为各种时间步长列出了相对于参考解决方案(FM)的MSOM和MSAT的L₂范数，并且还给出了压强和饱和度计算的自适应统计。即使在这种高度异构问题中，压强和饱和度的自适应比也与前面的例子相差不大。由于新饱和度自适应只有2-5％的区域需要延拓算子I(昂贵的施瓦兹重叠方法)，所以包括输运自适应计算的方法(MSAT)比原先多尺度有限体积算法(MSOM)高效得多。

这些例子证明了使用自适应输运计算的本发明多尺度结果与精细尺度解决方案极其一致。而且，流动和输运方程的自适应得出了在计算方面比以前的方法有效得多的算法。在本发明的范围内可以想到将粗尺度网格分段成数量不是三个，例如，2个和4个或大于4个(即4，5，6等)的区域，以标识可以将不同自适应延拓算子用于构建精细尺度饱和度的输运区域。但是，如上所述，发现三个输运区域足够了，并被用在这些例子中。类似地，在本发明的范围内想到可以利用非均匀网格尺寸、歪斜网格取向、交替网格单元形状、和多于一个二重粗尺度网格。但是，发现使用均匀初级粗尺度网格和二重粗尺度网格足够了。

虽然在本说明书中结合本发明的某些优选实施例对本发明作了描述，并且为了例示的目的阐述了许多细节，但对于本领域的普通技术人员来说，显而易见，可以不偏离本发明的基本原理地容易对本发明加以变更，并且可以相当大地改变本文所述的某些其它细节。