﻿ 裂缝性油藏非结构四边形网格数值模拟研究
 西南石油大学学报(自然科学版)  2018, Vol. 40 Issue (4): 105-115

“油气藏地质及开发工程”国家重点实验室·西南石油大学, 四川 成都 610500

Study of Quadrilateral Unstructured Grids in a Numerical Simulation of Fractured Reservoirs
LIU Yong , PENG Xiaolong, DU Zhimin
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China
Abstract: The introduction of unstructured quadrilateral meshes in the numerical simulation of fractured reservoirs is a popular area of research in the field of high-precision numerical simulation. These meshes can be used to create a grid with an unfixed number of connected components under the conditions of complex reservoir boundaries, cross fractures, multiple fractures, multi-directional fractures, complex multilateral wells, etc. Compared to other unstructured meshes, these meshes require fewer grids and have a higher precision. In this study, a mixed percolation model was established by combining the continuous medium model and the discrete fracture model. The numerical solution showed that the introduction of quadrilateral unstructured grids resulted in a significant decrease in the computational efficiency of the numerical simulation. Therefore, three solutions were proposed:morphing of the fracture endpoints with an improved paving algorithm, reordering of two-dimensional and threedimensional mesh numbers, and MPI parallel computing. The computational efficiency was then increased by an average of 70%.
Key words: unstructured quadrilateral mesh     computational efficiency     discrete fracture     mathematical model     reservoir numerical simulation

1 裂缝性油藏基本渗流数学模型

 图1 混合渗流数学模型流动空间示意图 Fig. 1 Schematic diagram of flow space in mixed seepage mathematical model
1.1 基质渗流方程

 $\nabla \cdot \left( {{\rho }_{{\rm o}}}{{v }_{{\rm o}}} \right)+{{q}_{{\rm o}}}-{{\delta }_{{{n}_{{\rm f}}}}}q_{{\rm o, }{{n}_{{\rm f}}}}^{*}=\dfrac{\partial \left( {{\rho }_{{\rm o}}}\phi {{S}_{{\rm o}}} \right)}{\partial t}$ (1)
 $\nabla \cdot \left( {{\rho }_{{\rm w}}}{{v }_{{\rm w}}} \right)+{{q}_{{\rm w}}}-{{\delta }_{n_{\rm f}}}q_{{\rm w, }{\rm n_f}}^{*}=\dfrac{\partial \left( {{\rho }_{{\rm w}}}\phi {{S}_{{\rm w}}} \right)}{\partial t}$ (2)

${{\rho }_{{\rm o}}}$${{\rho }_{{\rm w}}}—油、水的密度，{\rm kg/}{{{\rm m}}^{{\rm 3}}} {{v }_{{\rm o}}}$${{v }_{{\rm w}}}$—油、水的速度，${\rm m/s}$

${{q}_{{\rm o}}}$${{q}_{{\rm w}}}—油、水的质量流量，{\rm kg/(}{{{\rm m}}^{3}}\cdot {\rm s)} {{\delta }_{{{n}_{{\rm f}}}}}—Delta函数，当基质中流体直接向编号为{{n}_{{\rm f}}}的裂缝流动时取1，其他情况取值0； q_{{\rm o}, {{n}_{{\rm f}}}}^{*}$$q_{{\rm w}, {{n}_{{\rm f}}}}^{*}$—单位时间内体积基质块流向裂缝介质中的油、水的质量，${\rm kg/(}{{{\rm m}}^{3}}\cdot{\rm s)}$

$\phi$—孔隙度，%；

${{S}_{{\rm o}}}$${{S}_{{\rm w}}}—油、水的的饱和度，%； 下标{{n}_{{\rm f}}}—裂缝的序号，无因次。 油水两相的运动方程可以分别表示为  {{v }_{{\rm o}}}=-K\dfrac{{{K}_{{\rm ro}}}}{{{\mu }_{{\rm o}}}}\nabla \left( {{p}_{{\rm o}}}-{{\rho }_{{\rm o}}}{\rm g}d \right) (3)  {{v }_{{\rm w}}}=-K\dfrac{{{K}_{{\rm rw}}}}{{{\mu }_{{\rm w}}}}\nabla \left( {{p}_{{\rm w}}}-{{\rho }_{{\rm w}}}{\rm g}d \right) (4) 式中： K—渗透率，{\rm D} {{K}_{{\rm ro}}}$${{K}_{{\rm rw}}}$—油、水的相对渗透率，无因次；

${{\mu }_{{\rm o}}}$${{\mu }_{{\rm w}}}—油、水黏度，{\rm Pa}\cdot {\rm s} {{p}_{{\rm o}}}$${{p}_{{\rm w}}}$—油、水压力，${\rm Pa}$

${{\rho }_{{\rm o}}}$${{\rho }_{{\rm w}}}—油、水密度，{\rm kg/}{{{\rm m}}^{3}} {\rm g}—重力加速度，g = 9.8 {\rm m/}{{{\rm s}}^{2}} d—到静水液面的深度，{\rm m} 油水两相状态方程可表示为  \left \{ \begin{array}{l} {{C}_{{\rm o}}}=\dfrac{1}{{{\rho }_{{\rm o}}}}\dfrac{\partial {{\rho }_{{\rm o}}}}{\partial {{p}_{{\rm o}}}}\\[5pt] {{C}_{{\rm fo}}}=\dfrac{1}{\phi }\dfrac{\partial \phi }{\partial {{p}_{{\rm o}}}}\\[5pt] {{C}_{{\rm w}}}=\dfrac{1}{{{\rho }_{{\rm w}}}}\dfrac{\partial \phi }{\partial {{p}_{{\rm o}}}} \\[5pt] {{C}_{{\rm fw}}}=\dfrac{1}{\phi }\dfrac{\partial \phi }{\partial {{p}_{{\rm w}}}} \end{array} \right . (5) 式中：{{C}_{{\rm o}}}$${{C}_{{\rm w}}}$—油、水的压缩系数，${\rm MP}{{{\rm a}}^{-1}}$

${{C}_{{\rm fo}}}$${{C}_{{\rm fw}}}—随油、水压力变化的孔隙压缩系数，{\rm MP}{{{\rm a}}^{-1}} 辅助方程  \left \{ \begin{array}{l} {{S}_{{\rm w}}}+{{S}_{{\rm o}}}=1 \\ {{p}_{{\rm c}}}={{p}_{{\rm o}}}-{{p}_{{\rm w}}} \end{array} \right . (6) 式中：{{p}_{{\rm c}}}—毛管压力，{\rm MPa} 联立式(1)\sim式(6)，可得到基质系统完整的渗流模型  \left \{ \begin{array}{l} -\phi {{C}_{{\rm t}}}\dfrac{\partial {{p}_{{\rm o}}}}{\partial t}+\nabla \cdot {{\lambda }_{{\rm w}}}\nabla ({{p}_{{\rm o}}}-{{p}_{{\rm c}}})=\nabla \cdot {{\lambda }_{{\rm w}}}\nabla {{p}_{{\rm c}}}-\\{\kern 40pt} [{{q}_{{\rm wv}}}+{{q}_{{\rm ov}}}-{{\delta }_{{{n}_{{\rm f}}}}}(q_{{\rm o, }{{n}_{{\rm f}}}}^{*}+q_{{\rm w, }{{n}_{{\rm f}}}}^{*})] \\ \phi \dfrac{\partial {{S}_{{\rm o}}}}{\partial t}+{{S}_{{\rm o}}}\phi {{C}_{{\rm to}}}\dfrac{\partial {{p}_{{\rm o}}}}{\partial t}=\nabla \cdot {{\lambda }_{{\rm o}}}\nabla {{p}_{{\rm o}}}+{{q}_{{\rm ov}}}-{{\delta }_{n_{\rm f}}}q_{{\rm o, }{{n}_{{\rm f}}}}^{*} \\ {{S}_{{\rm o}}}+{{S}_{{\rm w}}}=1 \\ {{p}_{{\rm c}}}={{p}_{{\rm o}}}-{{p}_{{\rm w}}} \end{array} \right . (7) 式中：{{C}_{{\rm t}}}—油水两相总压缩系数，{\rm MP}{{{\rm a}}^{-1}} {{\lambda }_{{\rm o}}}$${{\lambda }_{{\rm w}}}$—油，水流度，${\rm D/(mPa}\cdot {\rm s)}$

${{q}_{{\rm wv}}}$${{q}_{{\rm ov}}}—水、油的体积流量，{{{\rm s}}^{-1}} {{C}_{{\rm to}}}—油相总压缩系数，{\rm MP}{{{\rm a}}^{-1}} 1.2 裂缝系统渗流方程 裂缝系统油水两相连续性方程可表示为  \dfrac{\partial {{({{\rho }_{{\rm o}}}{{v}_{{\rm o}}})}_{_{{{n}_{{\rm f}}}}}}}{\partial l}+q_{{\rm o, }{{n}_{{\rm f}}}}^{*}=\dfrac{\partial {{({{\rho }_{{\rm o}}}\phi {{S}_{{\rm o}}})}_{_{{{n}_{{\rm f}}}}}}}{\partial t} (8)  \dfrac{\partial {{({{\rho }_{{\rm w}}}{{v}_{{\rm w}}})}_{_{{{n}_{{\rm f}}}}}}}{\partial l}+q_{{\rm w, }{{n}_{{\rm f}}}}^{*}=\dfrac{\partial {{({{\rho }_{{\rm w}}}\phi {{S}_{{\rm w}}})}_{_{{{n}_{{\rm f}}}}}}}{\partial t} (9) 式中：l—裂缝上发生窜流所在网格边的边长，{\rm m} 在离散介质中，油、水所在的裂缝空间与整体空间相比为小量，故可将油、水密度视为常数处理，则式(8)，式(9)可简化为  \dfrac{\partial {{v}_{{\rm o, }{{n}_{{\rm f}}}}}}{\partial l}+q_{{\rm ov, }{{n}_{{\rm f}}}}^{*}={{\phi }_{{{n}_{{\rm f}}}}}\dfrac{\partial {{S}_{{\rm o, }{{n}_{{\rm f}}}}}}{\partial t} (10)  \dfrac{\partial {{v}_{{\rm w, }{{n}_{{\rm f}}}}}}{\partial l}+q_{{\rm wv, }{{n}_{{\rm f}}}}^{*}={{\phi }_{{{n}_{{\rm f}}}}}\dfrac{\partial {{S}_{{\rm w, }{{n}_{{\rm f}}}}}}{\partial t} (11) 裂缝系统油水两相的运动方程分别为  {{v }_{{\rm o, }{{n}_{{\rm f}}}}}={{\left[-K\dfrac{{{K}_{{\rm ro}}}}{{{\mu }_{{\rm o}}}}\dfrac{\partial \left( {{p}_{{\rm o}}}-{{\rho }_{{\rm o}}}{\rm g}d \right)}{\partial l}\right]}_{{{n}_{{\rm f}}}}} (12)  {{v }_{{\rm w}, {{n}_{{\rm f}}}}}={{\left[-K\dfrac{{{K}_{{\rm rw}}}}{{{\mu }_{{\rm w}}}}\dfrac{\partial \left( {{p}_{{\rm w}}}-{{\rho }_{{\rm w}}}{\rm g}d \right)}{\partial l}\right]}_{{{n}_{{\rm f}}}}} (13) 辅助方程  {{S}_{{\rm o, }{{n}_{{\rm f}}}}}+{{S}_{{\rm w, }{{n}_{{\rm f}}}}}=1 (14) 联立式(10)\sim式(14)，可得到裂缝系统完整的渗流模型  \dfrac{\partial }{\partial l}{{\left[({{\lambda }_{{\rm o}}}+{{\lambda }_{{\rm w}}})\dfrac{\partial p}{\partial l}\right]}_{{{n}_{{\rm f}}}}}+(q_{{\rm wv}, {{n}_{{\rm f}}}}^{*}+q_{{\rm ov, }{{n}_{{\rm f}}}}^{*})=0 (15)  \dfrac{\partial }{\partial l}{{\left({{\lambda }_{{\rm o}}}\dfrac{\partial p}{\partial l}\right)}_{{{n}_{{\rm f}}}}}={{\phi }_{{{n}_{{\rm f}}}}}\dfrac{\partial {{S}_{{\rm o, }{{n}_{{\rm f}}}}}}{\partial t}-q_{{\rm ov, }{{n}_{{\rm f}}}}^{*} (16) 式中：{{p}_{{}}}—压力，{\rm Pa} 1.3 基质-裂缝窜流方程 基质结点向编号为{{n}_{\rm f}}的裂缝结点流动方程为  q_{{\rm ov, }{{n}_{{\rm f}}}}^{*}={{\left[\dfrac{lh{{K}_{{\rm m}}}}{{{\mu }_{{\rm o}}}}\dfrac{({{p}_{{\rm o}}}-p)}{d'}\right]}_{{{n}_{{\rm f}}}}} (17)  q_{{\rm wv, }{{n}_{{\rm f}}}}^{*}={{\left[\dfrac{lh{{K}_{m}}}{{{\mu }_{{\rm w}}}}\dfrac{({{p}_{{\rm w}}}-p)}{d'}\right]}_{{{n}_{{\rm f}}}}} (18) 式中：d'—裂缝上发生窜流所在网格边与基质节点的垂线长度，{\rm m} h—油层的厚度，{\rm m} K_{\rm m}—基质渗透率，D。 \zeta =\dfrac{l\cdot h}{d'}，则式(17)和式(18)化为  q_{{\rm ov, }{{n}_{{\rm f}}}}^{*}={{\left[\zeta {{\lambda }_{{\rm o}}}({{p}_{{\rm o}}}-p)\right]}_{{{n}_{{\rm f}}}}} (19)  q_{{\rm wv, }{{n}_{{\rm f}}}}^{*}={{\left[\zeta {{\lambda }_{{\rm w}}}({{p}_{{\rm w}}}-p)\right]}_{{{n}_{{\rm f}}}}} (20) 建立基本渗流模型后，在数值求解时，采用非结构四边形为网格模型构造系数矩阵。组装的整体矩阵的奇异性与网格划分的质量直接相关，四边形网格离散质量好[6]，整体矩阵带宽相对就窄，矩阵的性质就好，即矩阵方程的奇异性较小，求解所需迭代步数减少，收敛性好。然而，非结构四边形离散时，由于裂缝宽度极小，自适应剖分在裂缝单元端点处产生大量极小化网格[7]，这类极小化网格宽度和裂缝宽度接近，导致构造系数矩阵时矩阵带宽增加，明显加大了系数矩阵的奇异性[8]，不利于数值求解，且可能产生不收敛的状况。消除极小化网格是求解裂缝性油藏渗流方程的关键。 2 改进Paving方法 Paving方法是一种直接的全四边形网格生成方法[9-10]，该方法的输入是区域边界离散所形成的永久边界。在内部网格生成的过程中[11]，边界分为内部边界和外部边界，如图 2所示，正在生成的网格由内部边界和外部边界同时按照固定节点和浮动节点自适应生成。Paving方法生成的网格质量好，效率高[12]，边界灵敏性好，内部不规则网格节点数较少[13]  图2 非结构四边形生成的Paving方法 Fig. 2 Paving method of unstructured quadrilateral generation 将Paving方法引入裂缝性油藏网格离散中，对同一个简单地质模型进行离散时，将直接Paving方法和直角网格适应裂缝宽度的方法作对比研究。方案如表 1所示，离散效果如图 3所示。由图 3a可见，在方案1中，裂缝的上下末端均产生了大量极小化网格；图 3b对方案1的结果直接进行了网格精细1倍处理，产生了大量极小化网格；图 3c直接使用直接网格适应裂缝宽度，产生了和裂缝宽度一样的正方形网格，数量巨大，图中只截取了其中的一部分。 表1 同一地质模型不同离散方法网格质量对比 Table 1 Grid quality comparison of different geological models with different discrete methods  图3 3种离散方法示意图 Fig. 3 Three discrete methods 将Paving算法引入裂缝性油藏中，在有小裂缝存在情况下[14]，需要对其进行改进。提出了自适应裂缝末端形态处理减少极小化网格。基本思路为通过程序自动化识别裂缝后，通过输入参数自动判断产生极小化网格的裂缝，进行裂缝末端自适应形态改进。 裂缝附近离散时，动态计算最短边，逼近裂缝分割出最大面积三角形或四边形以减少极小网格生成[15]，最长边缩放函数为  \left \{ \begin{array}{l} {{g}_{1{\max }}}{{\left( x \right)}}=\arctan x \\ {{g}_{2{\min }}}{{\left( x \right)}}=1-{{\rm e}^{-{{x}^{2}}}} \end{array} \right. (21) 式中：{{g}_{1{\max }}}$${{g}_{2{\min }}}$—最长边缩放函数，无因次；

$x$—到裂缝上发生窜流所在网格边的边长上的距离，$x\in \left( 0, l \right)$${\rm m} 目标优化函数  f\left( E \right)={{a}_{1}}{{g}_{{ j}1}}\left( \dfrac{{{l}^{2}}}{A} \right)+{{a}_{2}}{{g}_{{j}2}}\left( \dfrac{l}{U} \right)+\\{\kern 40pt}{{a}_{3}}{{g}_{{j}3}}\left( \varPsi \left( {{\varphi }_{1}}, {{\varphi }_{2}} \right) \right) (22) 式中：f\left( E \right)—表示剖分边的误差度量函数，无因次； {{a}_{1}}$${{a}_{2}}$${{a}_{3}}—剖分权重参数，无因次； {{g}_{j1}}$${{g}_{j2}}$${{g}_{j3}}—三角形三条边所用的缩放函数，无因次； A—所优化几何图形的面积，{{{\rm m}}^{2}} U—所优化裂缝所在两个几何图形中最短边的边长，{\rm m} \varPsi—所优化裂缝边夹角的最小值函数，{\rm rad} {{\varphi }_{1}}$${{\varphi }_{2}}$—所计算误差度量边两边的夹角，${\rm rad}$

${{a}_{1}}$=${{a}_{{\rm 2}}}$=${{a}_{{\rm 3}}}$=1，$j{\rm =1}$，则有

 $\left \{ \begin{array}{l} f\left( E \right)=\arctan \left( \dfrac{{{l}^{2}}}{A} \right)+\arctan \left( \dfrac{l}{U} \right)+\\{\kern 40pt}\arctan \left( \varPsi \left( {{\varphi }_{1}}, {{\varphi }_{2}} \right) \right) \\ {{l}_{\min }}=\Delta L \\ \Delta L=\dfrac{{{L}_{\max }}}{n} \end{array} \right .$ (23)

${{l}_{\min }}$—最短边目标长度，${\rm m}$

$\Delta L$—自适应最小边长，m；

${{L}_{\max }}$—控制剖分网格最大边长，m；

$n$—控制剖分网格的个数。

 图4 裂缝末端形态变形方法程序流程 Fig. 4 The procedure flow of the shape deformation method at the end of the crack

 图7 裂缝变形计算域的确定 Fig. 7 Determination of the calculation domain of fracture deformation

 图8 油藏模型 Fig. 8 Reservoir model

 图9 生产井的产油量和井底流压 Fig. 9 Oil production and BHP of production wells
 图10 注水井的注水量和生产井的含水率 Fig. 10 Water injection rate of Injection well and water cut of production wells

3 非结构四边形网格编号重新排序提高数值模拟计算效率

 图11 特殊处理三角形索引示意图 Fig. 11 Special processing triangle index schematic

 图12 重排序编程流程图 Fig. 12 Reorder programming flow chart

4 MPI并行计算

MPI是一种消息传递接口，即消息传递函数库的标准[23-24]。由美国并行计算中心发布，是大规模并行计算的环境之一[25]，根据并行MPI的构造方式，通过消息交换处理及计算模块搭建裂缝性油藏数值模拟的并行计算框架。

 图13 裂缝性油藏并行算法示意图 Fig. 13 Schematic diagram of parallel algorithm for fractured reservoirs

5 结语