近年来，页岩油开发成为热点^[1-2]。页岩油储层具有低孔低渗的特点，烃类化合物表现为自生自储，有机质结构与性质显著影响油气的流动^[3]。页岩致密储层经水平井和压裂技术改造，形成尺度较大并与天然裂缝网络沟通的裂缝网络，可提高其产能。人工形成裂缝的尺度与微纳孔隙裂缝相差很大^[4]，导致页岩油在不同尺度裂缝内的流动机理不同，存在表征流体单一、裂缝尺度范围大、划分网格要求精度高和流体参数在裂缝界面处不连续等困难，难以利用简单的解析模型表征油在整个裂缝网络中的流动，导致无法对页岩油藏复杂裂缝中的多相流体流动进行数值模拟表征。目前，广泛应用双重介质模型进行试井解释，但该模型并不适合有大尺度导流裂缝的油藏^[5-6]。早期学者使用单重介质模型模拟导流裂缝，但该模型需要极细密的网格描述裂缝和基质，计算量大、效率低，难以应用到实际的油藏数值模拟中^[7]。后来发展的离散裂缝模型假定裂缝内法向上的物理量保持不变，裂缝节点与基质节点在网格内重合，将裂缝及基质方程叠加后计算，实现了对单相流体流动的模拟，但不适用于多相流体流动。页岩油藏经过水力压裂后，导流裂缝的几何形状较复杂，需要沿裂缝划分非规则网格，但这类离散裂缝模型计算效率低。李宪文等人^[8]利用微地震数据，通过反演得到次生裂缝的分布特征，建立了压裂水平井渗流数学模型，定量表征了渗流场的变化规律。李玉梅等人^[9]利用离散元方法，分析了水力裂缝长度、天然裂缝倾角、内摩擦角及施工净压力对缝网的影响，但未考虑基质与裂缝间的渗流关系。同样，采用数值模拟方法模拟页岩油气藏时，也是利用传统试井方法进行反演，并未考虑多相渗流在页岩网络裂缝的流动^[10-12]。

黎水泉等人^[13-15]对离散裂缝模型进行了改进，提出了基于达西定律的交换流量耦合裂缝及基质方程，利用该方程可以准确地模拟油藏裂缝中的多相流动，但划分网格时要求沿着裂缝的走向划分，计算量大，求解大型矩阵出现奇异，不具有适用性。S. H. Lee等人^[16]把裂缝近似看作基质网格中的井源，提出了嵌入裂缝模型，利用Peaceman公式表征裂缝与基质网格间的流体交换，该模型网格划分效率高，但不能准确表征多相流体在裂缝和基质之间的交换。为此，笔者就多相流体在页岩油藏复杂裂缝中的流动进行了理论模型和数值求解方面的研究，建立了嵌入裂缝多相流体流动的控制模型和网格划分方法，实现了数值模拟的高效计算，为页岩油藏的产能评估提供了新的技术方法。

1.   嵌入裂缝多相流动模型的建立

1.1   嵌入裂缝多相流动模型

1.1.1   控制方程

页岩油藏中裂缝开度很小，其中流体的流动可以被看作一维流动。裂缝和基质中${\rm{\alpha }}$相流体的质量守恒方程分别为：

式中：${w_{\rm{f}}}$为裂缝宽度，m；${\tau _{\rm{f}}}$为裂缝走向；${q_{{\rm{fm,\alpha }}}}$为单位时间内从单位长度裂缝流入基质的${\rm{\alpha }}$（${\rm{\alpha }} = {\rm{w}},{\rm{o}}$）相流体的流量，m³/（m·s）；${v_{{\rm{f,\alpha }}}}$，${v_{{\rm{m,\alpha }}}}$分别为裂缝和基质中${\rm{\alpha }}$相流体的流速，m/s；${\rho _{\rm{\alpha }}}$为${\rm{\alpha }}$相流体的密度，kg/m³；${\mu _{\rm{\alpha }}}$为${\rm{\alpha }}$相流体的黏度，mPa·s；${\phi _{\rm{f}}}$，${\phi _{\rm{m}}}$分别为裂缝和基质的孔隙度；${K_{\rm{f}}}$，${K_{\rm{m}}}$分别为裂缝和基体的绝对渗透率，mD；${K_{{\rm{f,r\alpha }}}}$，${K_{{\rm{m,r\alpha }}}}$分别为裂缝和基质中${\rm{\alpha }}$相流体的相对渗透率；${p_{{\rm{f,\alpha }}}}$，${p_{{\rm{m,\alpha }}}}$分别为裂缝和基质中${\rm{\alpha }}$相流体的压力，MPa；${S_{{\rm{f,\alpha }}}}$，${S_{{\rm{m,\alpha }}}}$分别为裂缝和基质中${\rm{\alpha }}$相流体的饱和度；δ为克罗尼克函数；$\nabla$为拉布拉斯算子；t为时间，s；r为裂缝走向变量。

页岩油藏裂缝的开度很小，其中基质中的压力场在跨越裂缝时保持连续，多相流体在裂缝–基质中存在不同的流动特征，导致界面处相饱和度发生间断和基质的相饱和度场在裂缝处存在不连续，如图1所示。

图  1  油藏中垂直于裂缝方向上的油相饱和度分布曲线

Figure  1.  Saturation distribution curve of oil phase perpendicular to the fracture direction in oil reservoirs

下载: 全尺寸图片 幻灯片

裂缝与两侧基质的物质流动交换量并不相等，基质与裂缝的交换量可以表示为：

裂缝流入两侧基质的${\rm{\alpha }}$相流体流量$q_{{\rm{fm,\alpha }}}^ \pm$满足：

式中：$n_{\rm{f}}^ \pm$为裂缝两侧边界的外法线方向，其中$\pm$表示裂缝两侧。

1.1.2   方程离散

嵌入基质的裂缝被规则的基质网格分割为${N_{\rm{f}}}$段一维裂缝网格。页岩油藏裂缝宽度为${w_{\rm{f}}}$，取第$k$段裂缝网格为研究对象，裂缝内流体的流动被认为是一维的，采用一维离散格式把裂缝中${\rm{\alpha }}$相质量守恒方程全隐式离散为：

式中：$\Delta t$为时间步长，s；${{n}}$为$t = {{n}}\Delta t$时刻的变量；$\Delta {L_{{k}}}$为第$k$段裂缝网格的长度，m；${{\phi _{{{\rm{f}},k}}}}$为第$k$段裂缝网格的孔隙度。

$\left( {{v_{{\rm{f,\alpha }}}}} \right)_{{{k}} + \frac{1}{2}}^{{{n}} + 1}$和$\left( {{v_{{\rm{f,\alpha }}}}} \right)_{{{k}} - \frac{1}{2}}^{{{n}} + 1}$为通过裂缝网格2个界面的流速，可以表示为：

式中：各组分流体的相对渗透率$K_{{\rm{f,r\alpha }}}$使用上限值。

为实现裂缝与基质之间的流量交换，在裂缝两侧建立了位于裂缝网格垂直平分线上的参考点${M_ + }$和${M_ - }$，其与裂缝的距离分别为${d_ + }$和${d_ - }$（见图2）。由于压力场在裂缝附近连续分布，且其导数在单侧基质中也保持连续，所以可以通过裂缝及基质网格中的压力单侧插值确定参考点${M_ + }$和${M_ - }$上的压力，裂缝与两侧基质的流量交换$q_{{\rm{fm,\alpha }}}^ +$和$q_{{\rm{fm,\alpha }}}^ -$可以表示为：

图  2  裂缝网格两侧的参考点

Figure  2.  Reference points on both sides of fracture meshes

下载: 全尺寸图片 幻灯片

假设基质网格中的空间步长为$\Delta x = \Delta y = h$，以空间坐标为$(i,j)$的网格为研究对象，使用二维九点离散格式，基质中$\alpha$相质量守恒方程可以全隐式离散为：

式（12）中，各组分流体的相对渗透率${K_{{\rm{m,r\alpha }}}}$使用上限值，而系数$T$满足：

由于裂缝穿过基质网格，导致流体的物理量在网格中发生间断，该网格不参与实际计算，而是为其他网格的计算提供边界条件。裂缝两侧基质中的物理量分别保持连续，通过单侧插值的方法利用基质网格中的压力${p_{{{{\rm{m}},i - 1,j}}}}$和${p_{{{{\rm{m}},i,j + 1}}}}$及裂缝压力${p_{{{{\rm{f}},k}}}}$确定网格$(i,j)$中${S_1}$部分的压力场，得到网格$(i,j)$进入网格$(i - 1,j)$和$(i,j + 1)$的流量${q_{i - \frac{1}{2},j}}$和${q_{i,j + \frac{1}{2}}}$（见图3，图中的虚线界面可以看作整个油藏基质的内边界，且满足第二类边界条件），为网格$(i - 1,j)$和$(i,j + 1)$的流量计算提供边界条件。同理，该类网格与相邻常规网格的流量交换都可以使用上述方法求出。

图  3  内边界替代覆盖裂缝的基质网格

Figure  3.  Replacement of matrix meshes covering fractures by the inner boundary

下载: 全尺寸图片 幻灯片

求解裂缝控制方程时，若裂缝位于油藏内部，则在两端取绝流边界条件；若裂缝延伸至油藏边界，其两端的边界条件由油藏的外边界条件确定。

1.2   多裂缝交叉网络多相流动模型

页岩油藏中压裂裂缝网络中存在多裂缝交叉现象（见图4），在交叉点处多条裂缝相互连通，但由于交叉点体积仅为e²量级，远小于网格尺度，导致无法模拟多相流体在多裂缝交叉网络中的流动。设置交叉点体积仅为e²量级，可以避免在计算中出现奇异值。

图  4  交叉于一点的3条裂缝

Figure  4.  Three fractures intersecting at one point

下载: 全尺寸图片 幻灯片

为避免小体积交叉点而导致求解矩阵出现奇异值，利用裂缝网格间的界面传导率计算得到交叉处的交换流量表征裂缝网格交叉点的物理量。单相流动中，交叉点O处流体满足质量守恒方程：

式中：${v_{iO}}$为流体通过$i$条裂缝的流速，m/s；$l$为汇集到交叉点的裂缝网格数。

由于交叉点$O$的长度极小，即${L_O} \ll {L_i}$，通过裂缝$i$流入交叉点$O$的流速可以表示为：

将式（14）代入式（15），得到交叉点O处的压力为：

根据式（14）和式（15），可得单相流动中裂缝网格之间的流速为：

式中：${T_{ij}}$为裂缝网格$i$与$j$之间的界面传导率。

多相流动时，裂缝网格间的流速为：

2.   计算实例

参考页岩油藏开发中的一口实例井，采用停泵压降分析软件反演某页岩油藏一口多级压裂水平井的裂缝参数。某页岩油藏区域为2 000 m×2 000 m×100 m，水平井的水平段长1 000 m，水力压裂生成了20条裂缝；井储系数为0.2 m³/MPa。地层裂缝区域初始压力为40 MPa, 渗透率为4.93 mD，孔隙度为8.0%，地层温度为80 ℃；油藏日产液量40 m³，生产时间为365 d。页岩油藏相对渗透率曲线如图5所示。

图  5  页岩油藏相对渗透率曲线

Figure  5.  Relative permeability curve of the shale oil reservoir

下载: 全尺寸图片 幻灯片

地层中存在天然裂缝，但难以确定具体位置和尺寸。为此，假设水平井附近有单条裂缝、相互交叉2条裂缝和相交于一点的3条裂缝。水平井多级水力压裂产生的裂缝与天然裂缝如图6所示，其中水力压裂产生的裂缝半长最长为110 m，最短为50 m，平均有效裂缝半长为91.75 m。图6设置了与人工裂缝相交及不相交的单条裂缝、与人工裂缝相交的交叉裂缝、与人工裂缝不相交的交于一点的3条裂缝等3种情况，用于表征天然裂缝可能与人工裂缝构成的裂缝网络形式。

图  6  多级压裂水平井裂缝示意

Figure  6.  Fractures of a multistage fractured horizontal well

下载: 全尺寸图片 幻灯片

为验证建立的模型，计算多相流体在多裂缝交叉网络中流动时，需要从宏观物理量上反映有天然裂缝和无天然裂缝的区别，设以定产液量40 m³/d生产，计算其井底压力数据变化，验证人为设置裂缝对井底压力的影响。含有单条裂缝、相互交叉2条裂缝和相交于一点的3条裂缝的多级压裂页岩油水平井的网格划分如图7所示。人工裂缝和天然裂缝附近网格较密，不含有裂缝的区域网格稀疏。可以看出，水力压裂裂缝区域网格加密，可以实现水力压裂裂缝与天然裂缝的相交；多条裂缝相交于一点的裂缝也可以进行网格划分，并实现与人工裂缝的交叉。

图  7  含有不同类型裂缝的多级压裂水平井离散网格

Figure  7.  Discrete meshes of a multistage fractured horizontal well with different types of fractures

下载: 全尺寸图片 幻灯片

本算例中给出了定产量和初始压力，但页岩油藏地层压力随着生产时间增长而发生变化，存在生产时间较短、压力难以波及到天然裂缝和无法有效反映存在多条裂缝流体的流动引起的压力变化的问题。因此，假设图6中的水平段仅有人工裂缝，给定产液40 m³/d和初始压力40 MPa, 计算出初始时刻的地层压力（见图8），可以看出，多级压裂水平井裂缝的压力为基本上相同；压力波及范围区域如图9所示，可以看出，无论多级压裂水平井中的裂缝长短，缝内压力在40 MPa下的作用下压力影响区域边界大致相同。闷井10 d以后，缝内附近的地层压力下降，其压力波及范围快临近天然裂缝（见图10）。选此时压力35 MPa为计算压力，对不含有天然裂缝和含有天然裂缝的页岩油藏裂缝网络进行数值模拟，比较其压力变化。

图  8  不含有天然裂缝的多级水力压裂地层压力分布

Figure  8.  Formation pressure distribution of a multistage hydraulically fractured horizontal well without natural fractures

下载: 全尺寸图片 幻灯片

图  9  不含有天然裂缝的多级水力压裂水平井在初始压力40 MPa下压裂改造后的压力波及区域

Figure  9.  Pressure swept zone of a multistage hydraulically fractured horizontal well without natural fractures with an initial pressure of 40 MPa

下载: 全尺寸图片 幻灯片

图  10  不含有天然裂缝的多级水力压裂水平井在闷井10 d后的地层压力分布

Figure  10.  Formation pressure distribution of a multistage hydraulically fractured horizontal well without natural fractures after 10 d of shut-in

下载: 全尺寸图片 幻灯片

不含有天然裂缝和含有天然裂缝的多级水力压裂水平井的井底压力相对对数曲线如图11何图12所示（Ⅰ为井储段，Ⅱ为表皮段，Ⅲ为线性流段，Ⅳ为过渡段或线性流段，Ⅴ为线性流段，Ⅵ为边界流段）。

图  11  不含有天然裂缝的多级水力压裂水平井井底压力相对对数曲线

Figure  11.  Relative logarithmic curves of the bottomhole pressure of a multistage hydraulically fractured horizontal well without natural fractures

下载: 全尺寸图片 幻灯片

图  12  含有天然裂缝的多级水力压裂水平井井底压力相对对数曲线

Figure  12.  Relative logarithmic curves of the bottomhole pressure of a multistage hydraulically fractured horizontal well with natural fractures

下载: 全尺寸图片 幻灯片

从图11和图12可以看出，不含有天然裂缝和含有天然裂缝的多级水力压裂水平井井底压力在井储段存在明显差异，不含有天然裂缝的井底压力变化较快，而含有天然裂缝的井底压力变化持续时间较长，说明致密页岩油藏中存在裂缝沟通的情况，这与图6中存在天然裂缝与水力压裂裂缝沟通相一致。表皮段（Ⅱ）、线性流段（Ⅲ）、过渡段或线性流段（Ⅳ）和线性流段（Ⅴ）的压力变化趋势基本一致，其中过渡段或线性流段（Ⅳ）既呈现出过渡段的特征，又具有线性流的特征，在含有天然裂缝地层中的时间较短，这与天然裂缝与水力压裂裂缝的沟通有关。边界流段（Ⅵ）的压力变化趋势差异较大，不含有天然裂缝的多级水力压裂水平井井底压力相对对数曲线中的压力导数斜率较高，而含有天然裂缝的井底压力导数斜率较低，这主要是由缝内压力变化波及到天然裂缝区域，天然裂缝内的多相流体流向水力压裂裂缝。

3.   结　论

1）嵌入裂缝多相流动模型和多裂缝交叉网络多相流动模型，可以对页岩油藏复杂裂缝中的多相流体流动进行数值模拟表征，并根据地层压力变化规律得到页岩油藏随生产时间的变化情况。

2）利用裂缝与两侧基质的流量交换，实现复杂裂缝嵌入的网格划分，计算多相流体的流动。但油水在致密页岩油藏中的渗流速度较低，导致天然裂缝附近的压力变化与人工裂缝附近的压力变化相比较小。

3）利用井底压力变化特征可以表征天然裂缝与水力压裂裂缝的沟通，根据井储段、过渡段和边界流段的差异判断沟通水力压裂裂缝的天然裂缝规模及天然裂缝与水力压裂裂缝的距离，为评价页岩油藏储层提供科学依据。