DOI: 10.11817/j.issn.1672-7207.2015.12.027

非线性渗流多尺度有限元数值模拟

张娜¹，姚军²，黄朝琴²

(1. 中国石油大学(华东) 储运与建筑工程学院，山东 青岛，266555；

2. 中国石油大学(华东) 石油工程学院，山东 青岛，266555)

摘要：采用多尺度有限元法对一般的非线性椭圆型方程和与时间相关的非线性方程进行数值分析。阐述多尺度有限元法求解非线性方程的基本原理，推导相应的计算格式并编制计算程序。通过2个不同类型的非线性渗流算例验证该方法的准确性和适用性。最后，以低渗透油藏注水开发为例，对其进行多尺度有限数值模拟。研究结果表明：相对于传统的数值计算方法，多尺度有限元法具有更高的计算精度，在保证计算精度的同时能够减少计算量。

关键词：非线性渗流；多尺度有限元；多尺度；数值模拟

中图分类号：TE319             文献标志码：A         文章编号：1672-7207(2015)12-4584-08

Multiscale finite element method for nonlinear flow in porous media

ZHANG Na¹, YAO Jun², HUANG Zhaoqin²

(1. College of Pipeline and Civil Engineering, China University of Petroleum (Ease China), Qingdao 266555, China;

2. College of Petroleum Engineering China University of Petroleum(East China),Qingdao 266555, China)

Abstract: The multiscale finite element method was used to simulate general nonlinear elliptic equations and time-dependent nonlinear equations. The basic principle of the multiscale finite element method for solving nonlinear equations was described, and the corresponding calculation method was derived. The accuracy and the applicability of the proposed method were verified by two different kinds of nonlinear seepage flow. Finally, based on the study of the two models, the multiscale finite element numerical simulation of the water injection of the low permeability reservoir was carried out. The results indicate that MsFEM has higher accuracy compared to finite element method, and can reduce the computation time.

Key words: nonlinear flow; multiscale finite element method; multiple scale; reservoir numerical simulation

非线性现象广泛存在于自然界和人类社会等众多领域中。随着科学技术的发展，反映现实自然现象的非线性现象逐渐引起人们的极大关注，对非线性问题的研究日渐深入。非线性数学物理方程涉及水文地质、土木工程和石油工程等多个领域，例如非均质多孔介质中非饱和水流问题、 致密油气藏和低渗油藏渗流问题等。近几十年来，对这类非线性方程精确求解是许多国内外专家学者共同关注的问题。实际地层的渗透率分布通常具有很强的随机性, 地层各种参数在空间上具有显著的多尺度特征^[1]。因此，在求解非均质介质中非线性问题的主要困难在于方程系数呈现非线性和空间变异性，而计算精度在很大程度上取决于模型的精度。目前地质建模技术可生成精细地质模型，其中包含大量的地质信息数据，一般的地质模型可达数百万甚至数亿个网格单元，若采用传统的有限差分或有限元法求解^[2-7]，则需对小尺度信息进行分析，导致计算量巨大，很难实现数值模拟。在实际工程计算中，许多研究者采用网格粗化法模拟这类问题^[8]。网格粗化的目的是通过在小尺度上进行计算来预测大尺度上的有效性质，网格粗化后计算量将大大减少，但粗化后的大尺度模型数值模拟不能充分捕捉其小尺度特征。对此，近年来国内外不少研究者提出了一些多尺度数值模拟方法，如多尺度有限元法^[7]、非均质多尺度法。多尺度有限元法(multiscale finite element method，MsFEM)是基于有限元法发展起来的一种多尺度数值方法，不同于网格粗化法，该方法是以具有原分辨率的全局问题为目标并在粗网格上进行求解，因此，该方法在降低计算量的同时还能充分地捕捉到小尺度特征。MsFEM最初源自等^[9]的研究成果；1997年，Hou等^[8]将其推广到一般的二维振荡系数椭圆问题，提出了多尺度有限元法。对于MsFEM的研究主要集中于线性问题，并逐步扩展到多个领域，取得一定研究成果。Efendiev等^[10]将MsFEM进一步推广到非线性问题，对一般非线性椭圆问题进行了理论推导。近年来，MsFEM及其衍生方法主要集中于研究流体在高渗透区域的流动模拟^[11-15]，在低渗透领域的研究还很少。虽然MsFEM的研究时间较短，也未形成有效的通用软件，但其在模拟应用中具有许多独特的自身优势，例如可并行计算、自适应降尺度分析等。因此，本文作者基于Efendiev的理论基础，将MsFEM应用于求解地层中流体非线性渗流问题。

1  非线性椭圆型方程

传统的有限元(FEM)是20世纪50年代被提出的，随后得到迅速发展和逐步完善，具有成熟的理论基础。FEM进行数值模拟时首先将求解区域进行离散化，形成一系列有限单元，并假定单元内介质是均质；然后，采用多项式插值表示单元内的物理场特征；再利用Galerkin法等方法建立有限元方程，并将局部有限单元方程集合成整体有限元方程；最后，求解整体有限元方程，计算求解区域内物理场特征^[16]。FEM在选取基函数时通常采用线性基、二次基对真实物理场进行近似。当单元内介质为非均质、且物理场的变化是非线性的，这样描述势必产生较大的误差；若进一步精细划分单元网格，则相应的结点数增加，大大增加了计算量。所以，采用FEM模拟多尺度问题存在一定的局限性，主要在于选取的基函数难以真实描述物理问题的实质。

多尺度有限元法(multiscale finite element method，MsFEM) 就是针对多尺度问题提出的解决办法，它克服了FEM的不足，无需在小尺度上精确求解，能够在大尺度上捕获小尺度特征。MsFEM是在宏观尺度进行网格剖分，然后在单元上求解局部微分方程来构造多尺度基函数，最后在粗网格上求解原问题。这样不仅节约大量计算资源，同时保持较高的计算精度，还能够进行并行化。MsFEM与FEM的主要区别在于基函数的选取上。两者基函数对比如图1所示，其中，和分别FEM和MSFEM基函数。通过求解基于控制方程微分算子在单元上的局部方程得到MsFEM的基函数，这样自动将小尺度信息引入到大尺度范围下。若微分算子含有振荡项，则MsFEM基函数也具有相应的振荡特征(图1(c)) 。从图1可以看出：对于非均质单元，MsFEM基函数能够很好地反映单元参数变化，从而MsFEM采用较粗的单元剖分就能刻画出研究区域，从而降低计算量。

图1  FEM与MsFEM基函数对比示意图

Fig. 1  Comparison of basis functions of FEM and MsFEM

根据MsFEM，定义如下空间：

其中：为Soboler空间；u为Soboler空间变量；为p次方函数空间；为d维空间；α=(α₁，

α₂，…，α_d)为d维指数标；。

首先采用MsFEM对一般的非线性椭圆型方程进行分析：

其中：q为源顶；且在处满足p=0；和(，)满足一定条件以确保方程(1)的适定性。令是有界闭区域，n为边界的单位外法向向量，K^c为区域的一个网格剖分，C^c为有限维空间的一集族。MsFEM求解一般的非线性问题，需要建立多尺度映射和相应的多尺度公式。

2  多尺度有限元法

2.1  多尺度映射

不同于线性问题的多尺度有限元求解，非线性问题的多尺度基函数需要通过非线性映射将大尺度函数映射到小尺度上。以下面方程为例说明如何进行多尺度映射：

          (2)

对于，在任一粗网格单元内满足：

          (3)

若单元K为三角形单元且为分段线性函数，则p^f由p^c的结点值决定。对于一维情况，可以得到在区间[x_i-1，x_i]的表达式：

  (4)

其中：x_i-1和x_i为一维剖分的相邻结点。对于二维情 况，p^f可采用多尺度基函数来表示：

               (5)

其中：；x_i为结点；为多尺度基函数，且在粗网格单元K内满足

，        (6)

，              (7)

为分段线性插值函数；为的压力梯度。这样，p^f就可用多尺度基函数来表示。

2.2  多尺度公式

MsFEM求解的目的是得到满足p^f与小尺度解近似相等，因此p^f代入小尺度方程时，需要将得到的方程映射到大尺度上，而这一映射可通过小尺度方程乘以大尺度基函数来实现。故方程(1)相应的多尺度公式为

，      (8)

其中：

  (9)

2.3  超样本技术

为避免在求解多尺度映射局部问题时由于网格大小与物理小尺度大小相近而引起的共振误差，同时也为进一步提高MsFEM的精度及收敛速度，可采用超样本技术来求解非线性问题的基函数。 超样本技术的思路就是将原单元放大，先在放大的单元上求得临时基函数，然后根据临时基函数得到原单元上的基函数。

以图2中的单元为例，K为原始单元，K_E为放大后的单元，放大后的单元必须与原单元形状相似。首先求出单元K_E上的基函数，令单元K_E各点对应的基函数为ψ_i(i=1，…，4)，单元K各点对应的基函数为 (i=1，…，4)。临时基函数ψ_i满足

图2  超样本单元示意图

Fig. 2  Schematic description of sampling elements

           (10)

且，为Delta函数，当i=j，=1，否则，=0。由于K位于K_E内，所以，基函数可用临时基函数表示，它们满足如下关系式：

；i=1，…，4          (11)

且，由此可求解系数c_j，从而可以通过放大后的超样本单元求得原始单元对应的基函数。

3  与时间有关的偏微分方程

在非饱和多孔介质水流以及油藏多相渗流等问题中都会涉及到与时间t有关的偏微分方程，以在水文力学等领域应用十分广泛的Richards方程为例，讨论MsFEM对此类方程的求解。Richards方程可写为如下形式^[18]：

        (12)

式中：S为含水率；z为向上为正的垂直坐标；k为关于p的渗透率；ρ为流体密度；g为重力加速度。方程(13)的变分问题是求，满足

  (13)

其中：；，满足

，      (14)

，            (15)

为在时刻的近似值。Richards方程涉及到含水量S，每个时间步长上需要重新计算多尺度基函数，这与多尺度有限元法在两相流中的处理相似，相关证明可参考文献[17]。

4  数值算例

4.1  考虑压敏效应的非线性渗流

本文采用不精确牛顿法求解非线性方程。令为C^c的分段线性基函数，其中，N_dof为节点数，那么存在某一数组，MsFEM解可表达为

           (16)

且A与小尺度信息有关，其中，f为组成的向量。因此，求得的A需要满足

F(A)=0                 (17)

其中：F为非线性算子，满足

  (18)

可以发现，A隐式地存在于和p^f中。不精确牛顿法是牛顿迭代法的一种变形，用于求解非线性方程。该方法的求解思路就是给定一初始值A₀，对以下步骤循环迭代直到收敛：

1) ，且 ；

2) 。

其中：为第k步迭代的Jacobian矩阵；为第k次迭代的线搜索步长。可以发现，当β_k=0时，迭代过程为牛顿迭代法。本文选取

β₀=0.001，      (19)

下面对Jacobian矩阵进行研究。令={: z_i为单元K的1个顶点}，Jⁱ={j:z_j为单元的1个顶点}，那么式(19)可写为

  (20)

其中：p^f满足

，      (21)

，        (22)

Z_K为单元K所有顶点集合。F_i(A)并不是完全取决于，所以，存在，满足。令且剖分单元为三角形，根据链式法则，可得到关于的方程：

，

，

对于每个三角形单元K有，那么；若剖分单元为矩形单元，则将1/3换成1/4。

假设p^f已知，那么可根据上述方程得到。综上所述，可得到Jacobian矩阵元素：

，i=j        (23)

，       (24)

若采用超样本技术，则除了求解单元改变外，求解Jacobian矩阵元素步骤与上述过程相似。根据前面所述方程，可得到满足局部齐次方程。

令，则其满足以下局部方程：

，，

，

其中：在相应的单元K上通过计算得到；为超样本单元K_E的结点基函数，采有不精确牛顿法对相应的非线性方程进行求解。

首先考虑压力敏感储层的单相渗流问题，其渗流在区域Ω满足如下椭圆型方程：

            (25)

其中：；k₀为原始地层压力下的渗透率，μm²；p₀为原始地层压力，MPa；p为地层压力，MPa；α为储层压力敏感系数，1/MPa；μ为流体黏度，mPa·s。模型长×宽为1 m×1 m，上下边界为不渗透边界，左右边界分别为定压边界，原始地层压力为10 MPa，注入压力和采出压力分别为15 MPa和5 MPa，模型如图3所示。取α=0.05，μ=1 mPa·s。由于该模型无法得到解析解，本文采用FEM求得的精细解作为参考解。令N为各个方向粗网格剖分数目，FEM计算采用的网格剖分取N=512。图4所示为不同剖分情况下，粗网格单元上计算的传统有限元法FEM_c，MsFEM和超样本技术的多尺度有限元法MsFEM_o对应的相对误差。从图4可以看出：FEM_c的计算结果相对误差最大；随网格剖分数增多，MsFEM相对误差相应增大，这是共振效应所致；采用超样本技术的MsFEM_o明显提高了计算精度，而且避免了共振效应; MsFEM_o与参考解吻合较好，从而验证了该算法在处理椭圆型方程的正确性。

经计算发现，MsFEM比FEM节省了约90%的内存；MsFEM_o求解基函数时选取的超样本单元为原来单元的1.5倍，比FEM节省了70%~80%的内存。超样本单元放大倍数不同，得到结果精度不同。因此，本文将超样本单元分别放大为原来单元的1.0倍、1.2倍、1.5倍和2.0倍，对不同超样本单元放大倍数下的计算结果进行了对比，如图5所示。由图5可以发现：超样本单元放大倍数越大，相应的计算误差也越小；超样本单元放大到1.5倍以后相对误差近似相同，说明在放大到一定倍数后主导共振误差已消除。同时，超样本技术增加了计算基函数的计算量，放大倍数越大，计算量也相应增大；超样本单元越接近于原来的粗网格单元，计算耗时才会越短，也越接近MsFEM计算耗时。在对实际非线性问题进行数值模拟时，综合考虑准确性和计算耗时，选取MsFEM_o更适合。

图3  模型示意图

Fig. 3  Schematic of 2D model

图4  不同方法求解的相对误差

Fig. 4  Relative error of solutions using various schemes

4.2  非均质非饱和渗流

采用多尺度有限元法对二维非饱和非线性渗流问题进行模拟，求解以van Genuchten-Mualem模型^[19]为本构关系的Richards方程。van Genuchten-Mualem模型定义如下：

图5  不同放大倍数下MsFEM_o相对误差曲线

Fig. 5  Relative errors of solutions using MsFEM_o with different domain sizes

           (26)

  (27)

               (28)

其中：S_s为饱和流体含量；S_r为束缚流体饱和度；n′，m′和τ为模型参数。

模型参数取值如下：S_s=0.40；S_r=0.04；m′=0.5；n′=2；α=0.15 m^-1；饱和水力传导度K_s呈对数正态分布，均值为1 mm/s，方差为0.025。假设模型区域Ω是长×宽为100 m×100 m的正方形，初始压力水头为-22 m，上表面有水逐渐渗入，区域上边界压力水头为-22 m，下边界压力水头为-2 m，区域左右边界为不渗透边界。图6所示为传导度随机对数正态分布场。求解区域的细网格为250×250，计算中粗网格为25×25。

图7所示为在t=8 h时精细求解与多尺度有限元求解得到的压力水头。从图7可以看出：多尺度有限元与参考解结果吻合较好；多尺度有限元法能够在大尺度计算时，通过多尺度映射和局部更新来捕获小尺度特征。图8所示为在不同时刻多尺度有限元法和精细求解得到的水平方向平均含水率随深度的分布。

图9所示为在不同时刻多尺度有限元法与精细求解的计算时间对比。从图9可知：多尺度求解比精细求解节省了约40%的时间。多尺度有限元法未能在很大程度上减少计算时间，这是由于计算开始时，计算多尺度基函数耗时较长；而随后只需对多尺度基函数进行局部重建，节省了计算时间。从以上对比情况可以看出，多尺度有限元法在保证计算精度的同时节省了计算时间。

图6  随机对数正态传导度场分布图

Fig. 6  Distribution of lognormal random conductivity

图7  t=8 h时不同方法求解的压力水头对比

Fig. 7  Comparison of pressure head distributions obtained by different methods at t=8 h

图8  水平方向平均含水率对比

Fig. 8  Comparison of average water contents in horizontal direction

图9  精细求解和多尺度有限元法计算时间对比

Fig. 9  Comparison of computation time for fine olution and MsFEM solution

4.3  低渗透油藏两相渗流

假设考虑岩石的微可压缩性，忽略毛管压力和流体密度随压力的变化，则低渗透油藏两相流的数学模型为：

           (29)

      (30)

其中：θ表示油相w和水相o；v_θ为θ相流体流速，m/s；q_θ为θ相源汇项，1/s；S_θ为θ相饱和度，且S_o+S_w=1；k_rθ为θ相相对渗透率；λ_θ为θ相黏度，mPa·s；为θ相启动压力梯度，MPa/m。

根据所建立模型，以1/4五点注水的低渗透油藏为例进行研究。藏厚度为10 m，原始地层压力下的渗透率k₀=0.01 μm²，原始地层压力p₀=12 MPa，注入压力和采出压力分别为18 MPa和6 MPa，油的黏度μ_o=5 mPa·s，水的黏度μ_w=1 mPa·s，α=0.04。

油相相对渗透率k_rw=(S^*)²，水相相对渗透率k_ro=(1-S^*)²，其中，残余油饱和度S_or=0.4，束缚水饱和度S_wc=0.25。油相相启动压力梯度分别取0，0.01，0.03和0.05 MPa/m，相应的水相启动压力梯度分别取0，1，3和5 kPa/m。

图10和图11所示分别为不同压力梯度下，MsFEM计算得到生产井含水率和产油量与时间的关系曲线与参考解的对比。从图10和11可以看出：sFEM所得结果与参考解吻合较好，MsFEM在低渗透油藏两相渗流模拟中应用效果较好；启动压力梯度越大，含水率增长幅度越大，采收率也越低，因此，由于低渗透储层存在启动压力梯度导致这类油藏开发难度大。同时，采用MsFEM求解所耗费时间明显缩短，从而，MsFEM在保证计算精度的同时，大节省了计算时间。

图10  含水率随时间变化曲线

Fig. 10  Changes of water-cut with time

图11  产油量随时间变化曲线

Fig. 11  Changes of oil-production with time

5 结论

1) 采用多尺度有限元法对一般椭圆型方程和关于时间的非线性方程进行了研究。阐述了多尺度有限元法的基本理论，并推导了2类方程相应的多尺度数值计算格式。通过2类不同非线性渗流算例验证了该方法的正确性及适用性；并基于这2个算例，对低渗透油藏非线性两相流予以研究。计算结果表明：启动压力梯度对含水率和采收率有显著影响，增大了油藏开发难度。

2) MsFEM的求解结果与精细求解结果基本一致，在粗网格上求解便可捕捉到解的小尺度特征，从而大大节省了计算量又能保证计算精度，适用于求解一般的非线性渗流问题；MsFEM通过多尺度公式将小尺度信息引入到大尺度中，很好地反映了单元内参数变化导致的非线性性，在处理非均质、非线性问题上有显著的优越性。对于非线性渗流问题，MsFEM具有良好的适用性。

本文仅对二维的非线性问题进行了研究，对于三维多相非线性渗流问题有待进一步研究。

参考文献：

[1] HUANG Zhaoqin, YAO Jun, WANG Yueying, et al. Numerical simulation on water flooding development of fractured reservoir with discrete-fracture model[J]. Chinese Journal of Computational Physics, 2011, 28(1): 41-49.

[2] YAN Guoliang, WANG Diansheng, LIU Jinyu, et al. Numerical calculation of permeability for periodic fractured-vuggy porous media based on homogenization method[J]. Journal of the China Coal Society, 2011, 36(9): 1446-1450.

[3] Rutqvist J, Wu Y S, Tsang C F, et al. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2002, 39(4): 429-442.

[4] HUANG Zhaoqin, YAO Jun, LI Yajun, et al. Permeability analysis of fractured vuggy porous media based on homogenization theory[J]. Science China, 2010, 53(3): 839-847.

[5] HUANG Zhaoqin, YAO Jun, WANG Yueying, et al. Numerical study on two-phase flow through fractured porous media[J]. Science China, 2011, 54(9): 2412-2420.

[6] ZHI Qingong, ZHAO Yizhong, CHENG Yuanfang, et al. Finite element simulation on percolation of low permeability reservoir with permeability anisotropy[J]. Global Geology, 2012, 15(4): 265-269.

[7] YANG Renfeng, JIANG Ruizhong, LIU Shihua, et al. Numerical simulation of nonlinear seepage in ultra-low permeability reservoirs[J]. Acta Petrolei Sinica, 2011, 32(2): 299-306.

[8] Hou T Y, Wu X H. A multiscale finite element method for elliptic problems in composite materials and porous media[J]. Journal of Computational Physics, 1997, 134(1): 169-189.

[9] I, Osborn E. Generalized finite element methods: Their performance and their relation to mixed methods[J]. SIAM Journal on Numerical Analysis, 1983, 20(3): 510-536.

[10] Efendiev Y, Hou T Y. Multiscale finite element methods: Theory and Applications[M]. Heidelberg: Springer, 2009: 40-95.

[11] Gulbransen A F, Hauge V L, Lie K A. A multiscale mixed finite element method for vuggy and naturally fractured reservoirs[J]. SPE Journal, 2010, 15(2): 395-403.

[12] Efendiev Y, Galvis J, Wu X H. Multiscale finite element methods for high-contrast problems using local spectral basis functions[J]. Journal of Computational Physics, 2011, 230(4): 937-955.

[13] Chu C C, Hou T Y. A new multiscale finite element method for high-contrast elliptic interface problems[J]. Mathematics of Computation, 2010, 79(272): 1915-1955.

[14] Iliev O, Lazarov R, Willems J. Variational multiscale finite element method for flows in highly porous media[J]. Multiscale Modeling & Simulation, 2011, 9(4): 1350-1372.

[15] Efendiev Y, Galvis J, Hou T Y. Generalized multiscale finite element methods (GMsFEM)[J]. Journal of Computational Physics, 2013, 251(23): 116-135.

[16] 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003: 441-461.

WANG Xucheng. The finite element method[M]. Beijing: Tsinghua University Press, 2003: 441-461.

[17] 姚军, 张娜, 黄朝琴, 等. 非均质油藏多尺度混合有限元数值模拟方法[J]. 石油学报, 2012, 33(3): 442-447.

YAO Jun, ZHANG Na, HUANG Zhaoqin, et al. Multiscale mixed finite element numerical simulation for heterogeneous reservoir[M]. Acta Petrolei Sinica, 2012, 33(3): 442-447.

[18] 陈曦, 于玉贞, 程勇刚. 非饱和渗流Richards方程数值求解的欠松弛方法[J]. 岩土力学, 2012, 33(S1): 237-243.

CHEN Xi, YU Yuzhen, CHENG Yonggang. Under-relaxation methods for numerical solution of Richards equation of variably saturated flow[J]. Rock and Soil Mech, 2012, 33(S1): 237-243.

[19] VanGenuchten M T. A closed-form equation for predicting the hydraulic conductivity of un-saturated soils[J]. Journal of the Soil Science Society of America, 1980, 44(5): 892-898.

(编辑  赵俊)

收稿日期：2015-01-13；修回日期：2015-04-20

基金项目(Foundation item)：国家重点基础研究发展规划(973计划)项目(2011CB201004)；国家自然科学基金资助项目(51234007)；中国石油大学(华东)优秀博士学位论文培育计划项目(LW120201)；中国博士后科学基金资助项目(2015M582169)(Project (2011CB201004) supported by the National Basic Research Development Program (973 Program) of China; Project (51234007) supported by the National Natural Science Foundation of China; Project (LW120201) supported by the Outstanding Doctoral Dissertation Training Program of China University of Petroleum; Project (2015M582169) supported by the China Postdoctoral Science Foundation)

通信作者：姚军，博士，教授，从事复杂介质油气渗流理论和油藏数值模拟研究；E-mail：Rcogfr-upc@126.com