引言

裂缝性油藏数值模拟研究是油气田开发的重要开发技术之一，从地质建模到模型求解过程都非常复杂^[1]。在裂缝形成过程中，由于地质环境不断改变，导致裂缝的分布和裂缝的物性具有很强的非均质性，因而单孔介质模型和双重介质模型均不能很好地对裂缝性油藏进行精细描述。

离散裂缝模型对裂缝的描述可以适应多数复杂的裂缝分布^[2]。离散裂缝模型在三维地质建模中大量应用，但是在数值模拟领域起步较晚。首先，由于裂缝分布复杂，离散后的网格数量相对较多，不利于数值计算；其次，对于倾斜裂缝，一般只能选择非结构网格进行离散，离散后在每一条裂缝末端都产生了较多的极小化网格，增加了构造系数矩阵的复杂度，还会在数值求解时出现异常，导致代不收敛；最后，由于裂缝形状多样，形成的边界复杂，比如，两条裂缝交叉形成新边界、多段裂缝结合形成封闭边界等，这些复杂的边界使得非结构网格离散后，在边界附近产生更多的极小网格，进一步增加了数值模拟的难度。

基于以上问题，建立了由连续介质模型和离散裂缝模型相结合的混合渗流模型，引入非结构四边形网格解决存在复杂裂缝边界的网格离散问题^[3]，研究因非结构四边形网格导致计算效率低下的解决方法^[4]。本文从油藏工程数值模拟角度，基于非结构四边形网格提出了3种计算效率优化方法，并设计了相应的验证模型。

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

针对裂缝性油藏复杂渗流的特点^[5]，采用连续介质模型(基质、微裂缝)和非连续介质(大裂缝)模型对裂缝性油藏进行描述，建立连续双重介质和离散裂缝模型的混合渗流基本数学模型。同时，为了保证离散裂缝网格域和计算域上的无差异性，在非连续介质区域记录所有的裂缝宽度。不同宽度的裂缝构成的非连续介质离散介质区域如图 1所示，图中标红的多条裂缝复杂交叉几何图形为离散介质区域，白色为双重介质区域。

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]。

将Paving方法引入裂缝性油藏网格离散中，对同一个简单地质模型进行离散时，将直接Paving方法和直角网格适应裂缝宽度的方法作对比研究。方案如表 1所示，离散效果如图 3所示。由图 3a可见，在方案1中，裂缝的上下末端均产生了大量极小化网格；图 3b对方案1的结果直接进行了网格精细1倍处理，产生了大量极小化网格；图 3c直接使用直接网格适应裂缝宽度，产生了和裂缝宽度一样的正方形网格，数量巨大，图中只截取了其中的一部分。

将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}$；

下标$j$—缩放函数底标，$j \in \{1, 2\}$。

若${{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所示：程序先载入所有裂缝和网格的剖分参数，通过读取控制极小化网格参数的阈值，将此阈值作为程序判断极小化网格的标准，然后计算各个裂缝的自适应变形量，再将此变形量与周围动态生成的网格进行衔接得到最优化的比例，然后判断新生成的网格是否满足当前所设置的极小化网格的阈值，如果满足，则继续进行下一步，如果不满足，则重新设定变形量与衔接四边形的边长大小，循环执行直到全部满足极小化网格阈值的设定条件。

裂缝末端形态变形方法示意图见图 5，将正方形中的一条裂缝标记为矩形${ ABDC}$，在裂缝宽度方向的${ AB}$和${ CD}$两个末端进行形态变形，通过上述算法计算的自适应宽度将矩形变形为${ E}{{{ B}}^{'}}{{{ D}}^{'}}{F}{{{ C}}^{'}}{{{ A}}^{'}}$。

3种对极小化网格处理方式设置及结果如表 2所示，处理效果如图 6所示。图 6a为对裂缝上下末端都未进行处理的效果图，图上明显可见密集的极小化网格；图 6b对裂缝上端进行了处理，通过程序自适应处理后，在裂缝的上端已经没有极小化网格存在；图 6c对裂缝上下端均进行了处理，处理后整个区域极小化网格已经全部消失，证明了算法的优越性。

从处理前后后数据对比可以看出，处理后极小化网格明显减少，裂缝上下端均处理后极小化网减少到0，同时网格总数量也减少了35.71%。

这种变形处理后油藏模型实际网格域和计算域发生了改变^[16]，网格域和计算域如图 7所示，图 7a中根据前述自适应变形法则红色区域则为变形改变的几何区域${EB}{{{ B}}^{'}}$、${EA}{{{A}}^{'}}$、${FD}{{{ D}}^{'}}$、${FC}{{{C}}^{'}}$，图 7b所示为数值模拟计算程序对裂缝真实的计算区域矩形${ABDC}$。为了保证计算的精确性，离散裂缝模型的计算域与原模型保持一致，数模程序直接调用前述裂缝建模所记录的裂缝宽度。对于基质网格域的改变下面以实例进行论述。

设某油藏有压裂井：Well-1井、Well-2井和Well-3，油藏模型如图 8所示。

井和所在压裂缝基础模拟数据如表 3所示，网格域厚度离散设定为0.1 m，垂直渗透率为1 D，计算结果如表 4所示。

生产井和注入井生产数据在裂缝末端处理对比图见图 9，图 10。图 9所示为模型定产量开发，图 10所示为模型定注水量开发。

裂缝末端变形处理和未处理的井底流压曲线接近，指标相同(图 9)；裂缝末端变形处理与未处理的注水前期曲线和指标基本相同，处理后的结果在生产中期含水率稍有降低(图 10)。这是由于减少了极小化网格，迭代步数大大降低，累计误差减小，更接近真实的含水率，后期两种方法的含水率回归基本一致，可见裂缝末端变形处理保留了油藏含水率上升规律。因此，在基质-微裂缝的双重介质模型网格域的改变对计算结果的影响甚微，而整体计算效率却提高了近90%。

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

在裂缝性油藏数值模拟中，引入非结构四边形网格，优势明显，但是其存储和索引机制使得寻址效率极大降低。结构化网格采用数组顺序存放，通过数组下标直接访问，存储空间和访问时间消耗极少，而是非结构网格的编号具有随机性，存储方式复杂，除了需要存储网格编号外，还需要存储与相邻网格关系。以相同节点数的结构网格和非结构网格为例，数值离散后，网格组装系数矩阵求解通量项的个数比为3：7，因此，非结构网格对内存的消耗和对处理器的消耗是结构网的两倍多，再加上寻址效率问题，耗时更长。

非结构网格序号的无序性使得计算效率低下^[17]，对内存的需求量大，甚至会导致计算迭代不收敛。而非结构四边形网格在生成过程中，会同时存在四边形和三角形，编号更加复杂，若不重新排序，计算花费更大。采用非结构四边形网格的编号重新排序方法，可实现邻近网格的编号顺序最优化^[18]，从而降低数值模拟时计算机内存和处理器的消耗，提高裂缝性油藏数值模拟的效率。

非结构四边形网格相比其他非结构网格特点是平面产生的网格中包含四边形和三角形，对应生成三维模型时有六面体和五面体，在处理网格文件时采用一种新的记录方法使得三角形和四边形记录格式的一致性，方便程序统一排序。具体操作为：统一二维平面编号定为4个点，四边形按照节点编号逆时针排序，由于三角形只有3个点，因此将最后一个点复制作为重复标记的第四个点，同理，将三维立体编号统一定为8个点，将三角形五面体编号最后一位如法复制作为重复标记，如图 11所示。

排序步骤为：先二维排序然后三维排序，排序原理为依次循环取得各个网格编号后，计算编号邻差，取其平均值重新编号，同时将旧编号赋值给替换的新编号，直到各个网格编号邻差小于预先设置的排序阈值后停止重新排序。二维网格编号排序结束后进行三维网格编号重新排序(图 12)。

以本文算例为例，模拟油藏生产20 a，排序前计算步数为205，计算时间120 min，排序后计算步数为105，计算时间45 min，经过二维三维网格编号重新排序计算效率提高了62.5%。

4 MPI并行计算

裂缝性油藏网格数量庞大^[19]，许多复杂的裂缝性油藏模型在单机上无法顺利进行模拟计算^[20]，针对非结构四边形网格的特点，根据油田实际情况，多层油藏达到上万网格时，即便采用了裂缝末端变形处理和网格编号重新排序后，计算时间仍达到3 d以上，距离油藏数值模拟计算夕发朝至(16 h)的目标仍然较远，难以及时进行参数调节和优化，使得工程师无法快速对油气藏开发进行决策^[21]，而数百万的非结构四边形网格更难以计算。本文提出利用成熟的MPI并行计算方法以提高非结构四边形数值模拟计算效率^[22]。

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

油水两相裂缝性油藏数值模拟常规单机计算过程为：采用IMPES算法，同步交叉求解压力和饱和度^[26]，其中求解压力占据70%计算时间步。而在裂缝性油藏多层求解时，不同层由注水井和采油井分别联通，每一层油层的压力可以独立计算^[27-28]，从而可以将每一层的数值计算并行分发到各个计算节点进行求解，求解完各层压力即可求得饱和度。

根据此思想裂缝性油藏MPI并行算法流程图如图 13所示。

完成程序后，利用西南石油大学国家重点实验室大型刀片式并行机对10层算例模拟10 a计算，计算耗时如表 5所示

由表 5可见，应用并行计算方法后，程序运行时间由9.01 h降低到了2.86 h，加速比高达3.15。但计算效率仍然不高，主要原因是模型各层的求解效率依旧低下，裂缝性油藏非结构四边形网格并行计算方法的研究空间仍然巨大。

5 结语

引入了非结构四边形网格到裂缝性油藏油水两相数值模拟中。针对计算效率低下问题，提出了3种处理方法：改进Paving算法的裂缝末端形态变形法、非结构四边形网格编号重新排序方法和MPI并行计算方法。实例验证结果表明，对于单层网格，前二种方法计算效率平均提高76%以上，对于多层计算，综合运用3种方法平均计算效率提高了70%左右。本文提供的方法可有效的应用到裂缝性油藏非结构四边形网格的数值模拟计算中，以期探寻一体化解决复杂裂缝性油藏数值模拟的新方法。