引用本文
周彤, 王海波, 李凤霞, 李远照, 邹雨时, 张驰. 层理发育的页岩气储集层压裂裂缝扩展模拟[J]. 石油勘探与开发, 2020,47(5): 1039-1051.
ZHOU Tong, WANG Haibo, LI Fengxia, LI Yuanzhao, ZOU Yushi, ZHANG Chi. Numerical simulation of hydraulic fracture propagation in laminated shale reservoirs[J]. Petroleum Exploration & Development, 2020,47(5): 1039-1051.  
Doi: 10.11698/PED.2020.05.18
Permissions
Copyright©2020, 《石油勘探与开发》编辑部
《石油勘探与开发》编辑部 所有
层理发育的页岩气储集层压裂裂缝扩展模拟
周彤1, 王海波1, 李凤霞1, 李远照2, 邹雨时3, 张驰2
1.中国石化石油勘探开发研究院,北京 100083
2.中石化重庆涪陵页岩气勘探开发有限公司,重庆 408014
3.中国石油大学（北京）,北京 102249

联系作者简介：邹雨时（1985-）,男,博士,辽宁沈阳人,中国石油大学（北京）非常规油气科学技术研究院副研究员,主要从事岩石力学与水力裂缝扩展理论研究工作。地址：北京市昌平区府学路18号,中国石油大学（北京）289信箱,邮政编码：102249。E-mail：yushizou@126.com

第一作者简介：周彤（1986-）,男,博士,山东济宁人,中国石化石油勘探开发研究院高级工程师,主要从事非常规油气储集层压裂改造的相关研究工作。地址：北京市海淀区学院路31号,中国石化石油勘探开发研究院采油工程研究所,邮政编码：100083。E-mail：zhout1986@126.com

收稿日期: 2019-11-28
基金项目: 国家科技重大专项“页岩气水平井压后裂缝参数评价”（2016ZX05060001-032）; 中石化科技项目“涪陵页岩气水平井压后效果评价技术研究”（P18052-5）
摘要

以涪陵页岩气田焦石坝背斜主体区为研究对象,通过室内岩石力学实验和直剪实验明确了层理发育的岩石力学各向异性特征及其参数;结合室内实验评价结果,基于离散元方法建立了考虑力学各向异性、层理弱面和纵向应力差异的目标储集层三维压裂裂缝扩展模型,分析了不同层理弱面发育密度、层理强度以及压裂工程参数（射孔簇数、排量和压裂液黏度）条件下的水力裂缝展布规律。研究表明：综合考虑层理弱面和纵向应力差异的影响,研究区3~4 MPa的隔层应力差异将缝高控制在应力遮挡层内,缝高低于40 m;若不考虑层理弱面影响,缝高扩展预测结果明显较高。高密度层理缝的开启增加了水力裂缝复杂性,但显著限制了水力裂缝缝高的延伸。通过降低簇数、提高排量、增加前置液阶段高黏度压裂液用量及比例,可以减少层理弱面与纵向应力差异对缝高扩展的限制,促使裂缝纵向延伸。利用裂缝扩展模型对涪陵页岩气田焦页A井压裂施工段进行模拟,模拟结果与微地震裂缝监测结果相符。图13表4参34

关键词: 页岩; 层理; 水力压裂; 裂缝扩展规律; 纵向应力差异; 焦石坝背斜; 页岩气储集层
中图分类号:TE371 文献标志码:A 文章编号:1000-0747(2020)05-1039-13
Numerical simulation of hydraulic fracture propagation in laminated shale reservoirs
ZHOU Tong1, WANG Haibo1, LI Fengxia1, LI Yuanzhao2, ZOU Yushi3, ZHANG Chi2
1. Research Institute of Petroleum Exploration and Development, Sinopec, Beijing 100083, China
2. Shale Gas Exploration and Development Co., Ltd., Sinopec, Chongqing 408014, China
3. China University of Petroleum, Beijing 102249, China
Abstract

The main area of the Jiaoshiba anticline of the Fuling shale gas field was taken as the research object, laboratory rock mechanical experiments and direct shear experiments were conducted to clarify the mechanical anisotropy characteristics and parameters of rock samples with rich beddings. Based on the experimental results, a 3D fracture propagation model of the target reservoir taking mechanical anisotropy, weak bedding plane and vertical stress difference into account was established by the discrete element method to analyze distribution patterns of hydraulic fractures under different bedding densities, mechanical properties, and fracturing engineering parameters (including perforation clusters, injection rates and fracturing fluid viscosity). The research results show that considering the influence of the weak bedding plane and longitudinal stress difference, the interlayer stress difference 3-4 MPa in the study area can control the fracture height within the zone of stress barrier, and the fracture height is less than 40 m. If the influence of the weak bedding plane is not considered, the simulation result of fracture height is obviously higher. Although the opening of high-density bedding fractures increases the complexity of hydraulic fractures, it significantly limited the propagation of fracture height. By reducing the number of clusters, increasing the injection rate, and increasing the volume and proportion of high-viscosity fracturing fluid in the pad stage, the restriction on fracture height due to the bedding plane and vertical stress difference can be reduced, and the longitudinal propagation of fractures can be promoted. The fracture propagation model was used to simulate one stage of Well A in Fuling shale gas field, and the simulation results were consistent with the micro-seismic monitoring results.

Keyword: shale; lamina; hydraulic fracturing; fracture propagation law; longitudinal stress difference; Jiaoshiba anticline; shale gas reservoir
0 引言

页岩储集层具有低孔隙度、超低渗透率等特点, 水平井多级多簇压裂形成大规模复杂裂缝网络是有效开发页岩油气的关键。不同于常规储集层, 页岩作为层状沉积岩石, 具有显著的非均质性。同时, 页岩储集层中发育大量局部间断如断层、层理面和天然裂缝等, 使其具有显著的结构各向异性, 对水力压裂裂缝扩展产生一定影响。

目前, 研究页岩中复杂裂缝网络的扩展机制, 已成为非常规油气开发领域的重要课题。近年来, 国内外学者针对页岩储集层水力裂缝扩展进行了大量的理论研究, 建立能够模拟多裂缝扩展的压裂模型, 包括线网模型[1]、离散裂缝网络模型[2, 3]和非常规裂缝模型[4], 以及其他基于有限元法[5, 6]、边界元法[7, 8, 9, 10]、扩展有限元法[11, 12]、离散元法[13, 14, 15, 16]和相场法[17, 18]的模型。上述模型把岩石简化成各向同性材料, 重点研究高角度天然裂缝影响下的裂缝展布形态, 没有考虑层理、力学各向异性等页岩的典型特征对裂缝扩展形态的影响。室内压裂模拟与现场裂缝监测结果表明, 层理缝对水力裂缝垂向延伸及其扩展形态存在影响[19, 20, 21, 22, 23]。水力裂缝遇到层理面后可能产生贯穿、转向、终止或阶梯式延伸等结果, 从而导致最终裂缝网络形态的不确定性。因此, 多数数值模型无法体现出层理发育的页岩储集层与裂缝性储集层水力裂缝扩展的本质区别。

现场井下取心资料显示, 页岩储集层普遍发育高度密集的层理裂缝（每米数条至数百条）, 这种强各向异性特征给裂缝扩展模拟带来巨大挑战。利用有限元的黏聚单元分析天然裂缝弱面对压裂裂缝扩展的影响, 计算精度高[24], 但当天然裂缝密度预设过高时, 会导致计算稳定性、收敛性变差。对于弱面发育的页岩岩体, 其岩体介质更趋近于不连续介质的离散体。离散元法作为一种不连续介质数值模拟方法, 在处理岩土裂缝等不连续介质大变形问题（如高度复杂的裂缝网络）时具有较大优势[25]。Zhao等[13]应用二维颗粒离散元法模拟水力裂缝与天然裂缝的作用行为, 模型通过弱化颗粒间粘结强度来模拟天然裂缝。Zangeneh等[14]应用二维离散元法模拟水力裂缝网络, 模型中地层被多组节理所分割, 节理之间为可变形的岩块。Nagel等[15]应用三维离散元法研究了裂缝性储集层中水力裂缝类型及影响因素。但上述模型并没有考虑层理缝。Zou等[26, 27]建立了三维离散元裂缝网络模型, 探讨了工程尺度下层理缝对页岩储集层水力裂缝扩展的影响, 但是没有考虑储集层纵向应力非均匀性对水力裂缝扩展形态的影响。

为此, 本文以涪陵页岩气田焦石坝背斜主体区龙马溪组页岩储集层为研究对象, 开展了岩石力学各向异性特征实验测试, 基于离散元方法建立并求解了考虑页岩层理、岩石各向异性等因素的三维页岩储集层复杂裂缝扩展模型, 结合室内实验测试结果, 探讨了纵向应力差异与层理影响下的页岩储集层水力裂缝扩展规律。

1 区域概况

以焦石坝背斜主体区为研究目标, 其上奥陶统五峰组— 下志留统龙马溪组富有机质含气泥页岩段集中分布在五峰组— 龙马溪组底部, 底板为上奥陶统涧草沟组浅灰色瘤状灰岩。根据储集层物性特征, 将龙马溪组、五峰组页岩层从下到上共划分为9个小层。五峰组— 龙马溪组页岩段岩性主要为灰黑色炭质、硅质泥页岩, 主力含气层段层理发育, 自下而上层理发育程度逐渐降低[28]。其中, ①小层层理缝极发育, 每米具有上百条层理缝; 顶部⑧— ⑨小层层理缝发育程度较低, 取心观察结果见表1。储集层高角度天然裂缝总体发育规模小, 且相互之间不连通。

表1 不同小层层理发育特征

基于室内Kaiser地应力测试修正测井曲线, 解释得到焦页A井纵向应力剖面（见图1）。纵向上最小水平主应力存在差异, ①— ④小层平均地应力大小为49 MPa, 且层内应力差异小（仅为1~2 MPa）; ⑤小层地应力明显增大, 且与①— ④小层存在3~4 MPa的应力遮挡（达到53 MPa）; ⑥小层应力相对于①— ⑤小层进一步增大, 最高值接近56 MPa。另外, 目的层底板涧草沟组为较好的下部应力遮挡层, 隔层应力差大于12 MPa。

图1 焦页A井目的层段综合评价图

2 页岩各向异性特征

以涪陵页岩气田为研究对象, 利用高温高压动态岩石三轴测试系统, 对同一口井不同小层、不同取心方向的页岩试样进行三轴岩石力学测试, 结果如图2所示。测试结果表明：平行层理方向（0° ）取心时, 测得的弹性模量Eh值普遍高于垂直层理方向（90° ）取心的弹性模量Ev值。由层理缝极发育的①小层到层理缝极不发育的⑧小层, Eh/Ev值由1.26降至1.06（见图2a）。当加载应力垂直作用于层理面时, 外力压实作用会导致层理微裂隙闭合, 轴向应变较大, 获得的弹性模量与变形模量偏小, 且层理裂隙越发育, 这个现象越明显。所以随着围压的增加, 层间裂隙压实作用增强, 不同方向测得的弹性模量差异程度降低（见图2b）。另外, 受到层理裂隙的影响, 不同层理方向取心试样测试后破裂形态不同, 0° 取心页岩产生纵向的张拉劈裂裂缝, 而90° 取心试样则产生张性与剪切裂缝共存的共轭剪切破坏。

图2 不同取心方向页岩力学测试结果

为研究层理的力学性质及其影响下的页岩抗剪强度各向异性, 利用直剪试验机对平行层理方向与垂直层理方向的标准试样进行抗剪强度测试, 实验结果通过线性回归可以近似求得内摩擦角与黏聚力（见图3）。直剪实验获取的摩尔-库伦强度包络线表明：相同法向正应力条件下, 当剪应力方向平行层理方向时岩体发生剪切破坏所需要的剪应力低于垂直层理方向。同时, 层理发育程度越高, 平行层理方向剪切强度越低。层理极发育的①小层平行层理方向剪切时测得黏聚力为5.24 MPa, 层理发育程度较弱的⑧小层平行层理方向剪切时黏聚力为15.18 MPa。另外, 受层理发育程度影响, 平行层理方向剪切时裂缝破裂面断面平整、光滑; 垂直层理方向剪切时产生了基质张拉劈裂和层理面剪切滑移等复合破坏模式, 断面整体不平整, 伴随局部劈裂、剥离产生散体剥落碎片。

图3 不同剪切方向时的摩尔-库伦强度包络线

综上, 由于叠层间粘接弱面或层理的存在, 页岩岩石弹性参数和抗剪切强度等力学参数均具有各向异性。因此, 针对层理发育的页岩气储集层, 如何准确刻画层理弱面发育程度及其影响下的各向异性特征, 是体现页岩储集层与裂缝性储集层中水力裂缝扩展差异的关键。

3 三维裂缝扩展模型

模型控制方程主要由缝内压裂液流动方程、岩体变形方程和断裂破坏准则构成, 利用有限元和离散元的混合方法求解[26, 27]。根据离散元方法[29], 将地层模型求解域（求解域为Ω f）离散成若干个块体单元（三棱柱单元）, 块体单元之间通过虚拟弹簧链接, 传递相互作用力, 弹簧的断裂代表岩石的破裂。所有接触块体单元之间存在的节理单元构成压裂液流动的连通裂缝网络, 利用有限元法计算其内流体压力分布, 将此压力作为外部载荷作用在裂缝面上（即块体的接触面）, 然后计算块体变形和弹簧受力状态, 弹簧破裂（即裂缝扩展）由最大拉应力准则和摩尔-库伦准则决定。同时, 为了研究页岩岩石力学各向异性对裂缝扩展形态的影响, 在模型中利用横观各向同性线弹性材料本构方程来代替各向同性本构方程[26, 27]。

3.1 流动方程

缝内流体流动视为不可压缩牛顿流体的平面板流, 不考虑重力的作用, 满足连续性方程[28]：

$\begin{matrix} {{q}_{x}}=-\frac{{{w}^{3}}}{12\mu }\frac{\partial p}{\partial x}\\{{q}_{y}}=-\frac{{{w}^{3}}}{12\mu }\frac{\partial p}{\partial y}\\\end{matrix}$ （1）

$\frac{\partial w}{\partial t}+\frac{\partial {{q}_{x}}}{\partial x}+\frac{\partial {{q}_{y}}}{\partial y}+{{q}_{\text{l}}}=0$ （2）

由于页岩基质渗透率极低, 因此将基质块体视为不可渗滤体, 忽略压裂液的滤失, 即ql=0。那么, 在裂缝网络中全局的质量平衡方程如下：

$\int_{{{\Omega }_{\text{f}}}}{w\text{d}s}=tQ$ （3）

当N条裂缝起裂时, 进入各裂缝的流量之和为Q, 各裂缝的流量大小取决于各条裂缝在扩展延伸过程中的宽度和压力。模型假设缝内完全被流体填满, 且水力裂缝端部没有流动。

3.2 岩体变形方程

线弹性动态平衡方程如下[30]：

${{\sigma }_{ij, j}}+{{b}_{i}}-\rho {{a}_{i, t}}-\alpha {{u}_{i, t}}=0$ （4）

模型块体边界为固定边界, ui=0; 在块体接触面上施加接触力, 水力裂缝产生时在裂缝壁面上施加流体压力pi, ${{p}_{i}}={{\sigma }_{ij}}\times {{n}_{j}}$。由于应力应变关系符合线弹性本构方程[30], 那么：

${{\sigma }_{ij}}={{D}_{ijkl}}{{\varepsilon }_{kl}}$ （5）

页岩作为层叠沉积岩, 可视为横观各向同性材料, 即弹性特征在层理面内相同, 而垂直层理面方向不同[31]。用5个弹性常数表征横观各向同性岩石的线弹性特征, 包括Eh、Ev、υ h、υ v、Gv。如果页岩储集层地层是水平的, 则其柔度系数矩阵为：

$A=\left( \begin{matrix} \frac{1}{{{E}_{\text{h}}}} & -\frac{{{\nu }_{\text{h}}}}{{{E}_{\text{h}}}} & -\frac{{{\nu }_{\text{v}}}}{{{E}_{\text{v}}}} & 0 & 0 & 0 \\ -\frac{{{\nu }_{\text{h}}}}{{{E}_{\text{h}}}} & \frac{1}{{{E}_{\text{h}}}} & -\frac{{{\nu }_{\text{v}}}}{{{E}_{\text{v}}}} & 0 & 0 & 0 \\ -\frac{{{\nu }_{\text{v}}}}{{{E}_{\text{v}}}} & -\frac{{{\nu }_{\text{v}}}}{{{E}_{\text{v}}}} & \frac{1}{{{E}_{\text{v}}}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{2\left( 1+{{\nu }_{\text{h}}} \right)}{{{E}_{\text{h}}}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{{{G}_{\text{v}}}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{{{G}_{\text{v}}}} \\\end{matrix} \right)$ （6）

其中, 剪切模量Gv为[32]：

${{G}_{\text{v}}}=\frac{{{E}_{\text{h}}}{{E}_{\text{v}}}}{{{E}_{\text{h}}}\left( 1+{{\nu }_{\text{v}}} \right)+{{E}_{\text{v}}}}$

相邻的块体与块体之间存在法向与切向弹簧, 弹簧会由于拉伸或剪切滑移从而产生断裂, 即张性破坏或剪切破坏。弹簧的断裂根据最大拉应力准则与摩尔-库伦准则进行判断。当-Fn(n)< Fn, max（拉应力为负）或$|F^{(n)}_s|$< Fs, max时, 相邻接触点之间未发生破坏, 第n步的正应力与切应力根据（7）式和（8）式进行求解。

$F_{n}^{(n)}=F_{\text{n}}^{(n-1)}-{{k}_{\text{n}}}\text{ }\!\!\Delta\!\!\text{ }u_{\text{n}}^{(n)}$ （7）

$F_{\text{s}}^{(n)}=F_{\text{s}}^{(n-1)}-{{k}_{\text{s}}}\text{ }\!\!\Delta\!\!\text{ }u_{\text{s}}^{(n)}$ （8）

当Fn(n)≥ Fn, max=AT0时, 发生张性破坏, 弹簧断裂, 相邻块体发生分离（Δ un> 0）, 相互作用力消失, 即：

$\left\{ \begin{matrix} F_{n}^{(n)}=0 \\ F_{\text{s}}^{(n)}=0 \\\end{matrix} \right.$ （9）

当$F^{(n)_s}≥ {{F}_{\text{s, max}}}=A{{S}_{0}}+{{F}_{n}}^{\left( n \right)}\text{tan}\varphi$时, 发生剪切破坏, 相邻的块体与块体之间产生切向位移（Δ un≤ 0, Δ us> 0）。此时, 不同于张性破坏, 块体间还存在压应力与摩擦阻力, 即：

$\left\{ \begin{align} & F_{\text{n}}^{(n)}=F_{\text{n}}^{(n-1)}-{{k}_{\text{n}}}\text{ }\!\!\Delta\!\!\text{ }u_{\text{n}}^{(n)} \\ & F_{\text{s}}^{(n)}=F_{\text{s, re}}^{{}}=F_{\text{n}}^{(n)}\tan \varphi \\ \end{align} \right.$ （10）

3.3 方程迭代耦合

在缝内流体压力的作用下, 基质块体发生变形后改变了裂缝的宽度; 随着缝宽的变化, 缝内流量也相应发生改变, 继而影响缝内流体压力的大小, 即缝内流体压力与裂缝缝宽相互影响。采用弱耦合方法来实现裂缝内压裂液流动和固体变形的迭代过程, 分别离散流动方程和岩体变形方程, 并顺序循环求解[26, 27]。基质块体在发生张性或剪切破坏前, 初始节理单元的渗透率与基质块体渗透率相等, 因此将基质块体渗透率等效为裂缝的初始宽度, 即w0=（12K0）1/2, 为流体最初流动设计原始缝宽[33]。然后, 在当前时间步长内, 需要选取合适的实验解pm和wm, 求解下一时间步的压力pm+1, 进而由（11）式求取wm+1。通常取0< β ≤ 0.5, 时间步长足够小有利于方程迭代收敛[34]。

$\left\{ \begin{matrix} {{p}_{m+\tfrac{1}{2}}}=\frac{1}{\gamma \left( {{w}_{m}} \right)}\left( \frac{\Delta {{w}_{m}}}{\Delta {{t}_{m}}}-f \right) \\ {{p}_{m+1}}=(1-\beta ){{p}_{m}}+\beta {{p}_{m+\tfrac{1}{2}}} \\\end{matrix} \right.$ （11）

3.4 模型验证

利用模型分别模拟了简单双翼对称垂直裂缝和水平径向缝的扩展形态, 并与经典压裂模型解析解进行对比, 主要输入参数包括E=35 GPa, υ =0.2, μ =5 mPa· s, Q=5 m3/min, 结果如图4所示。由于经典PKN模型和径向模型假设裂缝扩展时的能量消耗主要用于裂缝内流体流动, 不考虑岩石断裂韧性、抗拉强度等力学破裂参数的影响, 因此数值模型中抗拉强度设定为零。数值模型计算结果与PKN和径向模型计算结果基本吻合, 差异分别为3.6%和4.9%, 从而验证了模型的可靠性。

图4 数值模型与经典模型模拟结果对比

4 模拟结果与分析

基于涪陵页岩气田焦石坝背斜主体区储集层参数, 建立层状地层模型, 考虑地层为水平产状, 网格单元（三棱柱）边界沿着预设层间界面迹线（层理弱面）。前期开展①— ⑨小层裂缝扩展模拟时, 发现裂缝不受上部⑦— ⑨小层影响。由于模型计算效率的问题, 下文建模时未考虑⑦— ⑨小层的影响。模型在X、Y、Z方向分别宽200 m, 长500 m, 厚50 m。根据室内实验与测井地应力解释结果, 将9个小层分为3段纵向上存在应力差异的层位（见表2）。其中, 底部①— ④小层总厚度为28 m; 中部⑤小层厚度为10 m, 与⑥小层应力差为3 MPa; 上部⑥小层厚12 m。水平井穿过①— ④小层中部, 网格模型如图5所示。为模拟单段3簇压裂, 假设每个射孔簇仅产生一条水力裂缝。基础压裂模拟参数为：一段3簇, 簇间距25 m, 排量14 m3/min, 压裂液黏度2.5 mPa· s。

表2 地应力剖面设置参数

图5 最小水平主应力网络模型

研究区层理缝整体发育, 尤其是底部五峰组— 龙马溪组①— ④小层, 每米达上百条层理缝。在裂缝扩展模型模拟过程中, 需要采用等效方式将层理缝近似成压前“ 显式” 的层理弱面, 层理弱面间隔为2~10 m。同时, 为了对比研究天然裂缝对裂缝扩展形态的影响, 在模型中设置了离散的高角度天然裂缝, 缝长20 m, 缝高10 m, 随机离散分布。由于储集层高角度天然裂缝总体发育规模小, 且相互之间不连通, 因此天然裂缝线密度设定为0.05条/m。参考室内岩石力学测试结果, 确定模拟输入的页岩基质、天然裂缝与层理弱面的抗拉强度、黏聚力和内摩擦角等力学参数（见表3）。

表3 压裂模拟主要输入参数
4.1 不考虑层理的影响

为说明非均质层段中层理弱面对人工裂缝展布的影响, 首先以均质储集层段（不考虑层理和天然裂缝）水平井单段压裂裂缝扩展模拟形态作为对比（见图6a）。在均质层段中应用上述压裂施工参数进行模拟, 压后形成3条单翼主裂缝, 缝长450~480 m。外侧簇裂缝在应力阴影影响下, 向外侧偏转; 中间簇裂缝在外侧裂缝挤压下, 沿缝高方向的延伸情况优于外侧裂缝。近井筒缝高穿过⑤、⑥小层界面, 缝高达50 m。随着至井筒距离增大, 人工裂缝高度逐渐减小, 由于受到隔层应力差异的限制, 可见在缝尖端处的缝高限于④小层顶部, 无法进入⑤小层, 高度仅为28 m。考虑天然裂缝影响（不考虑层理）, 随着天然裂缝的开启, 水力裂缝复杂度提高, 主裂缝缝长减小, 约330 m。由于局部穿层高角度天然裂缝的存在, 局部缝高增大（见图6b）。整体来看, 不考虑层理时近井筒裂缝均贯穿两个应力遮挡层。

图6 不考虑层理影响时水平井单段压裂裂缝扩展模拟形态（模拟注入15 min）

4.2 考虑层理的影响

4.2.1 层理弱面密度的影响

室内压裂物理模拟发现, 页岩压后会形成层理与水力裂缝交织而成的“ 栅栏形” 复杂缝网[16], 层理弱面密度会显著影响人工裂缝扩展形态。基础压裂施工模拟参数条件下（Nc=3, Q=14 m3/min, μ =2.5 mPa· s）, 当层理弱面为中等强度（TBP=4 MPa, SBP=15 MPa, φ BP=25° ）、弱面间距为2 m时, 不同注入时刻裂缝扩展形态如图7所示。扩展前期（模拟注入5 min）, 水力裂缝在①— ④小层内延伸, 同时伴随层理开启（见图7a）; 注入15 min, 受隔层应力差的影响, 缝高扩展受限从而稳定在28 m, 水力裂缝沿缝长方向延伸, 伴随着层理的持续开启（见图7b）; 注入45 min, 中间簇水力裂缝突破④小层并扩展至⑤小层, 最终被层理裂缝截止, 缝高稳定在38 m（见图7c）。

图7 考虑层理影响时水平井单段压裂不同模拟时刻裂缝扩展形态（弱面间距为2 m）

图8为不同层理弱面密度时的人工裂缝扩展数值模拟结果。层理作为大面积、连续的弱面, 可以增加水力压裂改造的裂缝密集度, 从而提高储集层改造的充分性与改造效果。层理密度越大, 单位改造体积内水力裂缝密度越大, 但过多层理的开启会严重制约水力裂缝在长度与高度方向上的扩展, 储集层改造体积大幅度减小。不同层理密度条件下, 人工裂缝长度、高度与施工时间的关系如图9所示。当dBP为8 m时（弱发育）, 缝长与缝高快速扩展, 15 min时近井缝高穿过⑤、⑥小层界面, 之后整体人工裂缝高度稳定在40 m, 最终模拟缝长达400 m; 当dBP为2 m时缝长小于300 m, 水力裂缝截止于⑤、⑥小层界面内, 近井筒处裂缝最大高度为37 m。总体上, 近井裂缝纵向扩展主要截止于⑤、⑥小层界面, 缝高约为38 m, 明显小于不考虑层理时的情况, 而远端人工裂缝纵向截止于④、⑤小层界面, 缝高约为28 m。

图8 不同层理密度时裂缝扩展形态（模拟注入45 min）

图9 不同层理密度时的裂缝扩展动态

4.2.2 层理弱面强度的影响

层理强度是对层理容易开启与否的描述, 在数值模型中利用层理的抗张强度和剪切强度（包括黏聚力、内摩擦角）来表征。为了研究层理强度对裂缝扩展形态的影响, 实验设置4组对比方案, 通过输入不同的TBP、SBP和φ BP值来模拟强度不同的层理弱面, 强度分级见表4所示。

表4 层理弱面强度力学参数表征

层理间距为2 m, 裂缝扩展动态参数与模拟结束时的裂缝形态如图10和图11所示。低强度层理弱面页岩储集层人工裂缝被近井层理裂缝截断, 仅沿着层理水平扩展, 缝高小于5 m（见图11a）。相比之下, 层理力学强度中等— 高时缝高扩展顺利, 当水力裂缝穿过④小层后, 受⑤小层隔层应力差异影响缝高停止扩展一段时间, 压裂模拟后期（40 min）中间簇裂缝在中高— 高层理强度条件下穿出纵向应力遮挡层（见图10）。因此, 随着层理强度的提高, 层理开启难度增加, 压裂液造主裂缝效率更好, 缝长增加（见图11d）。层理强度较弱时, 层理容易被水力裂缝开启, 裂缝沿层理转向、分流, 缝高扩展受限。

图10 不同层理强度条件下裂缝扩展动态

图11 不同层理强度条件下的裂缝扩展模拟结果（模拟注入45 min）

4.2.3 压裂工程参数的影响

为考察压裂施工参数与裂缝扩展形态的相关性, 数值模拟过程中恒定地质参数, 考虑中等强度层理（TBP=4 MPa, SBP=15 MPa, φ BP=25° ）, 层理间距为2 m; 改变簇数、排量与压裂液黏度, 模拟结果如图12所示。

图12 不同工程参数条件下的裂缝扩展模拟结果（模拟注入45 min）

在相同注入参数条件下, 单段簇数越少, 越有利于提高缝内流体压力, 促进人工裂缝缝高的扩展。在簇数为2, 簇间距为25 m时, 形成2条横切裂缝, 缝长为400 m, 近井人工裂缝纵向扩展至接近⑥小层顶部, 缝高达48 m, 而远井人工裂缝纵向截止于⑤、⑥小层界面, 缝高为38 m（见图12a）。随着簇数的增加, 缝高减小。簇数为3（见图7c）和5（见图12b）时, 近井人工裂缝纵向截止于⑤、⑥小层界面, 缝高38~40 m, 远井人工裂缝纵向截止于④、⑤小层界面, 缝高28 m, 缝长约为300 m。与3簇压裂相比, 5簇压裂时沿缝长方向的缝高大小的维持能力变弱, 远离井筒裂缝的缝高很快降低至28 m, 缝高整体变小。

注入排量与压裂液黏度也是页岩储集层体积压裂改造需要考虑的重要工程参数。当排量由12 m3/min增加到16 m3/min, 人工裂缝穿层扩展趋势明显, ⑤小层改造范围增加（见图12c和图12d）; 提高压裂液黏度, 可以显著降低层理和纵向隔层应力差异对缝高扩展的限制作用。当黏度为1.0 mPa· s时, 水力裂缝充分开启了近注入点层理, 人工裂缝纵向截止于④、⑤小层界面, 缝高28 m; 同时由于层理裂缝分流作用显著, 主裂缝缝长仅168 m, 大幅度降低单井缝控储量, 不利于长期稳产（见图12e）。当压裂液黏度提高至2.5 mPa· s时, 缝高到达⑤小层, 不同位置的层理均得到不同程度的开启。当压裂液黏度为25 mPa· s时, 水力裂缝沿最大主应力方向扩展（垂直裂缝）, 缝高突破至⑥小层（见图12f）; 但提高压裂液黏度大幅度降低层理缝的开启程度, 降低裂缝复杂性。

页岩储集层水力压裂以提高裂缝复杂性、增大储集层改造体积为目标, 从而获取产能最大化。根据模拟结果来看, 高度发育的层理缝提高了水力压裂裂缝的复杂性, 但同时制约了裂缝缝高的延伸, 减小储集层改造体积。因此, 对于层理发育、开启容易的压裂层段, 需要降低簇数、提高排量、增加前置液阶段高黏度压裂液用量及比例, 从而降低层理对缝高扩展的限制; 而对于层理不发育、开启难的压裂层段（层理不发育的上部气层或上覆应力大的深层页岩气储集层）, 则是以提高裂缝复杂性、增加泄气裂缝面积为压裂目标, 需要增加簇数、提高排量、降低压裂液黏度（采用低黏度滑溜水）。

5 与微地震监测结果的对比

以涪陵页岩气田焦页A井为例, 该井水平段主要穿行①— ③小层, 压裂施工共计28段（2簇压裂16段、3簇压裂12段）, 簇间距17~25 m, 平均单段用液量1 965 m3, 支撑剂用量50.7 m3。压裂过程中对28段进行微地震监测, 各段微地震事件数5~135个。监测结果显示裂缝长度为229~483 m, 缝网带宽34~204 m, 裂缝高度16~49 m（其中缝高30~40 m的裂缝占52%）。利用三维裂缝扩展数值模型对各段进行模拟, 与微地震裂缝监测结果对比（见图13, 取监测的第3段进行展示）。模拟得到的各段裂缝缝长为250~430 m, 缝高25~50 m, 带宽76~167 m, 整体情况较为相符。

图13 模拟结果与微地震监测对比示意图（红点表示微地震信号）

6 结论

基于离散元方法建立了层理发育的页岩储集层压裂裂缝扩展模型, 考虑了页岩储集层层理弱面和纵向应力差异的影响。采用涪陵页岩气田焦石坝背斜主体区龙马溪组页岩气水平井实际参数进行数值模拟分析, 研究发现高密度层理开启可以增加改造体积的裂缝复杂性, 但缝高、缝长明显受抑制; 低强度层理开启, 人工裂缝转向水平扩展, 井底压力较低, 缝高扩展受限。不考虑层理影响时, 存在7 MPa隔层应力差异, 中间簇局部水力裂缝缝高可以穿过应力遮挡, 缝高接近50 m; 考虑层理作用时, 存在3~4 MPa隔层应力差异, 缝高受控于应力遮挡层内, 缝高低于38 m。因此, 预测页岩储集层裂缝扩展时需要充分考虑层理弱面与纵向应力差异的影响, 否则会使缝高扩展预测结果与实际存在差异。对于层理异常发育的压裂层段, 需要降低簇数、提高排量、增加前置液阶段高黏度压裂液用量及比例, 降低层理对缝高扩展的限制, 从而提高裂缝纵向延伸高度、增大储集层改造体积。

符号注释：

ai, t— — 加速度, m/s2; A— — 接触面积, m2; A— — 柔度系数, Pa-1; bi— — 单位体积的力, N/m3; dBP— — 弱面间距, m; Dijkl— — 弹性张量, Pa; E— — 各向同性均质储集层弹性模量, Pa; Eh, Ev— — 平行层理和垂直层理方向的弹性模量, Pa; f— — 右端向量, m/s; Fn— — 法向力, N; Fn, max— — 块体单元间弹簧发生张性断裂时所需要的法向力, Pa; Fs— — 切向力, N; Fs, max— — 块体单元间弹簧发生剪切断裂时所需要的切向力, Pa; Fs, re— — 弹簧断裂后的残余剪切阻力, N; Gv— — 剪切模量, Pa; GR— — 自然伽马, API; i, j, k, l— — 张量指标, 无因次; kn, ks— — 法向和切向弹簧刚度, N/m; K0— — 基质渗透率, m2; m, m+1— — 当前迭代步和下一迭代步; n-1, n— — 当前时间步和下一时间步; nj— — 裂缝面法向单位向量, 无因次; N— — 裂缝条数, 条; Nc— — 簇数; p— — 压力, Pa; pi— — 裂缝壁面施加的流体压力, Pa; pm— — 当前迭代步的压力实验解, Pa; pm+1/2— — 当前迭代步到下一迭代步过程中的压力实验解, Pa; p0— — 孔隙压力, Pa; ql— — 滤失速率, m/s; qx, qy— — x与y方向上的流速, m2/s; Q— — 总注入量, m3/s; RLLD— — 深侧向电阻率, Ω · m; RLLS— — 浅侧向电阻率, Ω · m; s— — 裂缝面积, m2; SBP— — 层理缝剪切强度, MPa; S0— — 基质剪切强度, Pa; t— — 时间, s; TBP— — 层理缝抗拉强度, MPa; T0— — 基质抗拉强度, Pa; w— — 动态裂缝宽度, m; w0— — 初始缝宽, m; wm— — 当前迭代步的缝宽实验解, m; x, y— — 坐标系两个方向, m; ui— — 位移, m; ui, t— — 速度, m/s; α — — 单位体积内的阻尼, （kg· s）/m3; β — — 迭代修正系数; γ — — 实验解求解系数, m/（Pa· s）; ε kl— — 应变张量, 无因次; μ — — 流体黏度, Pa· s; υ — — 泊松比, 无因次; υ h, υ v— — 平行层理和垂直层理方向的泊松比, 无因次; ρ — — 岩石密度, kg/m3; σ h— — 水平最小主应力, Pa; σ H— — 水平最大主应力, Pa; σ ij— — Cauchy张量, Pa; σ ij, j— — Cauchy张量导数, N/m3; σ v— — 垂向应力, Pa; φ — — 内摩擦角, （° ）; φ BP— — 层理裂缝内摩擦角, （° ）; Ω f— — 求解域; Δ tm— — 第m步内的自适应时间步的步长, s; Δ wm— — 当前迭代步的缝宽变化, m; Δ un, Δ us— — 相邻节点之间法向和切向相对位移, m。

编辑 刘恋

参考文献
[1] XU W, THIERCELIN M, GANGULY U, et al. Wiremesh: A novel shale fracture simulator[R]. SPE 132218, 2010. [本文引用:1]
[2] MEYER B R, BAZAN L W. A discrete fracture network model for hydraulically-induced fractures: Theory, parametric and case studies[R]. SPE 140514, 2011. [本文引用:1]
[3] NAGEL N, IVAN G, MARISELA S N, et al. Simulating hydraulic fracturing in real fractured rock: Overcoming the limits of pseudo 3D models[R]. SPE 140480, 2011. [本文引用:1]
[4] WENG X, KRESSE O, COHEN C E, et al. Modeling of hydraulic fracture network propagation in a naturally fractured formation[R]. SPE 140253, 2011. [本文引用:1]
[5] BAO J Q, FATHI E, AMERI S. A coupled finite element method for the numerical simulation of hydraulic fracturing with a condensation technique[J]. Engineering Fracture Mechanics, 2014, 131(2): 269-281. [本文引用:1]
[6] SHIN D H. Factors controlling the simultaneous propagation of multiple competing fractures in a horizontal well[M]. USA: Prentice Hall, 2014. [本文引用:1]
[7] OLSON J E. Multi-fracture propagation modeling: Applications to hydraulic fracturing in shales and tight gas sand s[R]. ARMA 08327, 2008. [本文引用:1]
[8] 胥云, 陈铭, 吴奇, 等. 水平井体积改造应力干扰计算模型及其应用[J]. 石油勘探与开发, 2016, 43(5): 780-786.
XU Yun, CHEN Ming, WU Qi, et al. Stress interference calculation model and its application in volume stimulation of horizontal wells[J]. Petroleum Exploration and Development, 2016, 43(5): 780-786. [本文引用:1]
[9] 周彤, 陈铭, 张士诚, 等. 非均匀应力场影响下的裂缝扩展模拟及投球暂堵优化[J]. 天然气工业, 2020, 40(3): 82-91.
ZHOU Tong, CHEN Ming, ZHANG Shicheng, et al. Simulation of fracture propagation and optimization of ball-sealer in-stage diversion under the effect of heterogeneous stress field in a horizontal well[J]. Natural Gas Industry, 2020, 40(3): 82-91. [本文引用:1]
[10] 陈铭, 张士诚, 胥云, 等. 水平井分段压裂平面三维多裂缝扩展模型求解算法[J]. 石油勘探与开发, 2020, 47(1): 163-174.
CHEN Ming, ZHANG Shicheng, XU Yun, et al. A numerical method for simulating planar 3D multi-fracture propagation in multi-stage fracturing of horizontal wells[J]. Petroleum Exploration and Development, 2020, 47(1): 163-174. [本文引用:1]
[11] DAHI T A, OLSON J E. Numerical modeling of multi-strand ed hydraulic fracture propagation: Accounting for the interaction between induced and natural fractures[J]. SPE Journal, 2011, 16(3): 575-581. [本文引用:1]
[12] GORDELIY E, PEIRCE A. Implicit level set schemes for modeling hydraulic fractures using the XFEM[J]. Computer Methods in Applied Mechanics & Engineering, 2013, 266: 125-143. [本文引用:1]
[13] ZHAO X, YOUNG R. Numerical simulation of seismicity induced by hydraulic fracturing in naturally fractured reservoirs[R]. SPE 124690, 2009. [本文引用:2]
[14] ZANGENEH N, EBERHARDT E, BUSTIN R M. Investigation of the influence of natural fractures and in situ stress on hydraulic fracture propagation using a distinct-element approach[J]. Canadian Geotechnical Journal, 2015, 52(7): 926-946. [本文引用:2]
[15] NAGEL N B, SANCHEZ M A, ZHANG F, et al. Coupled numerical evaluations of the geomechanical interactions between a hydraulic fracture stimulation and a natural fracture system in shale formations[J]. Rock Mechanics and Rock Engineering, 2013, 46(3): 581-609. [本文引用:2]
[16] ZOU Y, ZHANG S, MA X, et al. Numerical investigation of hydraulic fracture network propagation in naturally fractured shale formations[J]. Journal of Structural Geology, 2016, 84: 1-13. [本文引用:2]
[17] MIEHE C, HOFACKER M, SCHÄNZEL L M, et al. Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 294(1): 486-522. [本文引用:1]
[18] 刘国威, 李庆斌, 左正. 相场断裂模型分步算法在ABAQUS中的实现[J]. 岩石力学与工程学报, 2016, 35(5): 1019-1030.
LIU Guowei, LI Qingbin, ZUO Zheng. Implementation of a staggered algorithm for a phase field model in ABAQUS[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(5): 1019-1030. [本文引用:1]
[19] 柳占立, 庄茁, 王涛, 等. 页岩水力压裂的关键力学问题[J]. 固体力学学报, 2016, 37(1): 39-49.
LIU Zhanli, ZHUANG Zhuo, WANG Tao, et al. The key mechanical problems on hydraulic fracture in shale[J]. Chinese Journal of Solid Mechanics, 2016, 37(1): 39-49. [本文引用:1]
[20] 张士诚, 郭天魁, 周彤, 等. 天然页岩压裂裂缝扩展机理试验[J]. 石油学报, 2014, 35(3): 496-503.
ZHANG Shicheng, GUO Tiankui, ZHOU Tong, et al. Fracture propagation mechanism experiment of hydraulic fracturing in natural shale[J]. ActaPetroleiSinica, 2014, 35(3): 496-503. [本文引用:1]
[21] 许丹, 胡瑞林, 高玮, 等. 页岩纹层结构对水力裂缝扩展规律的影响[J]. 石油勘探与开发, 2015, 42(4): 523-528.
XU Dan, HU Ruilin, GAO Wei, et al. Effects of laminated structure on hydraulic fracture propagation in shale[J]. Petroleum Exploration and Development, 2015, 42(4): 523-528. [本文引用:1]
[22] 衡帅, 杨春和, 曾义金, 等. 页岩水力压裂裂缝形态的试验研究[J]. 岩土工程学报, 2014, 36(7): 1243-1251.
HENG Shuai, YANG Chunhe, ZENG Yijin, et al. Experimental study on hydraulic fracture geometry of shale[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(7): 1243-1251. [本文引用:1]
[23] ZOU Y, ZHANG S, ZHOU T, et al. Experimental investigation into hydraulic fracture network propagation in gas shales using CT scanning technology[J]. Rock Mechanics and Rock Engineering, 2016, 49(1): 33-45. [本文引用:1]
[24] CHEN Z, JEFFREY R G, ZHANG X, et al. Finite-element simulation of a hydraulic fracture interacting with a natural fracture[R]. SPE 176970, 2017. [本文引用:1]
[25] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29(1): 47-65. [本文引用:1]
[26] ZOU Y, MA X, ZHANG S, et al. Numerical investigation into the influence of bedding plane on hydraulic fracture network propagation in shale formations[J]. Rock Mechanics and Rock Engineering, 2016, 49: 3597-3614. [本文引用:4]
[27] ZOU Y, MA X, ZHOU T, et al. Hydraulic fracture growth in a layered formation based on fracturing experiments and discrete element modeling[J]. Rock Mechanics and Rock Engineering, 2017, 50: 2381-2395. [本文引用:4]
[28] 马永生, 蔡勋育, 赵培荣. 中国页岩气勘探开发理论认识与实践[J]. 石油勘探与开发, 2018, 45(4): 561-574.
MA Yongsheng, CAI Xunyu, ZHAO Peirong. China’s shale gas exploration and development: Understand ing and practice[J]. Petroleum Exploration and Development, 2018, 45(4): 561-574. [本文引用:2]
[29] MUNJIZA A. The combined finite-discrete element method[M]. Chichester: John Wiley and Sons, 2004. [本文引用:1]
[30] JAEGER J C, COOK N G W, ZIMMERMAN R. Fundamental of Rock Mechanics[M]. 4th ed. Oxford: Blackwell Publishing Ltd. , 2009. [本文引用:2]
[31] HENG S, YANG C, ZHANGB , et al. Experimental research on anisotropic properties of shale[J]. Rock and Soil Mechanics, 2015, 36(3): 606-616. [本文引用:1]
[32] LEKHNITSKII S G. Theory of elasticity of an anisotropic body[M]. Russian: Mir Publishers, 1981. [本文引用:1]
[33] ZIMMERMAN R W, BODVARSSON GS. Hydraulic conductivity of rock fractures[J]. Transport in Porous Media, 1996, 23(1): 1-30. [本文引用:1]
[34] ADACHI J, SIEBRITS E, PEIRCE A, et al. Computer simulation of hydraulic fractures[J]. International Journal of Rock Mechanics & Mining Sciences, 2007, 44(5): 739-757. [本文引用:1]