我国页岩油气资源丰富，开发潜力巨大。近年来，以水平井+水力压裂开发工艺为基础、以地质工程一体化为核心的立体建模开发技术的应用，大大提升了页岩油气的钻探开发效率。但压裂时地层孔压变化与天然断层裂缝、压裂缝网共同作用形成扰动应力场，井眼周围应力实时改变^[1–2]，甚至出现应力反转，影响压裂缝网扩展方向和形态，改变油气从岩石基质向井筒的流动效率；同时，还会引起地层形变，激活断层^[3–6]，使断层失稳错动引发套管变形^[7]。N.P. Roussel等人^[8–9]结合位移不连续法和有限元法，发现应力场动态变化会导致孔隙体积与岩体介质压缩，对流体流动和裂缝扩展产生影响。Zhang Guangqing等人^[10]发现地应力差和起裂角会促使新裂缝的方向突破应力屏障后移动至最大水平应力方向。Su Yu等人^[11]发现压裂时局部地应力动态增大导致套管根部有效应力增加，更容易发生损坏。Lu Zongyu等人^[12]认为压裂地层有效应力减小与断层滑移是造成套管变形的主要原因。因此，准确分析储层不同阶段地应力场的实时状态^[13]，是理解和预测页岩储层压裂过程中裂缝扩展行为、防止地层损伤、优化压裂设计的关键。

针对页岩储层压裂与生产过程中地应力动态变化问题，国内外学者先后提出了多种渗流−应力耦合模拟方法，建立了数值模型。P. Hatchell^[14]、J. Herwanger^[15]、A. Onaisi等人^[16]将储层模拟器与地震驱动的地质力学模型耦合，构建了考虑生产或注入过程时间效应的四维动态地应力模型。D. Yang等人^[17]和Guo Xuyang等人^[18]采用有限元方法，研究了不同水力裂缝参数及岩石力学参数对页岩储层应力变化的影响。张滨海等人^[19]基于有限差分算法，提出了四维地应力渗流−应力耦合建模分析方法。张洪军等人^[20]利用ECLIPSE−VISAGE耦合提出了特低渗透油藏的动态裂缝模拟方法，开展了四维应力场耦合模拟分析。朱海燕等人^[21–22]采用ECLIPSE−ABAQUS交叉迭代耦合方法，建立了渗流−应力耦合四维地应力演化及复杂裂缝扩展的多物理场模型。李奎东^[23]以涪陵页岩气立体井组为例，建立了渗流−应力交叉耦合的动态地质力学模型，模拟了老井压裂复杂裂缝扩展及压后长期开采四维地应力演化过程。然而，上述四维地应力数值模拟方法存在明显不足：1）多采用交叉迭代或单向耦合方法，只提高了计算效率，过程中多假设理想量，计算精度不足；2）渗流与地质力学模型采用不同的处理器独立处理并相互赋值，在耦合模拟过程中由于网格和地质特征，物理场参数映射易产生误差，且需要更新一个物理场参数作为另一个物理场的计算条件，模拟过程繁琐且精度受限。

基于上述现状，笔者综合考虑页岩层理的非均质性、地质力学参数及压裂缝网扩展的复杂性，提出了基于有限差分算法的压裂动态地应力渗流−应力全耦合建模与演化模拟方法，以整合地质建模软件的地层描述能力和有限元软件的多场耦合计算优势，实现三维地质力学模型到四维动态地应力有限元模型的有效转换，建立区域渗流与地质力学共同控制全耦合模型，并在多孔介质中同时使用不同的控制方程描述物理场，实现多物理场效应在同一次迭代过程中的全耦合和动态响应。同时，结合现场施工数据系统分析了压裂过程中地应力场时空演化规律，以期为页岩储层压裂方案设计与现场施工提供参考。

1.   页岩储层压裂渗流−地质力学耦合理论

1.1   地应力演化全耦合模拟方法

页岩储层压裂地应力动态演化渗流−应力全耦合分析方法，综合考虑了地层非均质性、各向异性、物性、地质力学参数及实际压裂制度，可实现从初始地应力场到压裂引发应力扰动的全过程模拟，分析真实地层条件下压裂地应力场的动态变化规律。该方法首先利用井层数据、地震数据、测井资料等，以岩石力学与物性参数实验数据为校核点，构建包含物性及地质力学参数的压裂区三维精细化地质模型。针对压裂缝网扩展，利用压裂微地震监测点云数据处理，准确校正压裂改造区缝网扩展范围与空间位置。其次，基于三维静态地质力学模型提取压裂区初始地应力场、孔隙压力分布以及三维属性参数场作为耦合模拟的材料属性与初始条件。

为保证现场适用性和计算精度，采用了全耦合计算方法，通过自编程提取转化精细化地质力学模型中的地质构造和网格属性数据，重构有限元地质模型。结合渗流−应力耦合控制方程设置物理场与边界条件，利用等效材料处理对断层与压裂缝网（微地震数据）穿过的网格赋予渗流材料属性，在渗流耦合场模拟计算时采用孔渗参数应力敏感性模型，实时更新区域渗流物性。同时，以初始孔隙压力与地应力场作为有限元地质力学模型初始力学边界条件^[24]，以渗流模型计算所得压裂渗流场与孔隙压力变化数据作为地应力耦合模型边界条件，最终建立真实地层尺度下包含地质力学属性与渗流物性的页岩储层四维动态地应力模型，结合现场施工数据研究压裂行为对井区地层压力与地应力的影响规律。

1.2   渗流−应力耦合数学模型

由于页岩地层的特殊性，本文采用等温达西渗流−固体力学全耦合数学模型，假设孔隙介质中流体为弱可压缩性流体，多孔介质为小变形与各向同性弹性材料，全耦合模拟时将多物理场控制方程组合，构建固体力学应力场与渗流流体场耦合联合控制方程，由流体控制、岩体控制以及孔渗实时变化方程共同描述。

1.2.1   岩体控制方程

考虑到页岩的岩性特征，采用小变形与各向同性弹性假设下的多孔介质变形，结合Biot等效弹性模型，建立多孔介质有效应力与本构方程：

式中：${\boldsymbol{\sigma}}$为总应力张量，Pa；${{\boldsymbol{\sigma}} }_{\mathrm{e}}$为有效应力张量，Pa；${\alpha }_{{\mathrm{B}}}$为biot系数；$p$为孔隙流体压力，Pa；${\boldsymbol{I}}$为笛卡尔度规张量；${K}_{\mathrm{d}}$为排水后多孔介质整体体积模量，Pa；${K}_{\mathrm{s}}$为骨架材料的无孔隙固体体积模量，Pa。

将弹性体本构方程的应力张量换成有效应力张量，把孔隙压力的贡献添加到总应力张量中，可得储层岩体应力场本构方程。联合小变形多孔弹性模型的控制方程组，岩体控制方程采用替换后的多孔介质应力场本构方程、平衡微分方程、几何方程：

式中：${\boldsymbol{C}}$为弹性张量；${\boldsymbol{\varepsilon}}$为小应变张量；$\nabla$为矢量微分算子，1/m；$\nabla \cdot {\boldsymbol{\sigma}} \mathrm{为}\nabla$和${\boldsymbol{\sigma}}$张量的点积运算，Pa/m；${{\boldsymbol{F}}}_{\mathrm{V}}$为体积力矢量，N；${\boldsymbol{U}}$为位移矢量，m；$\nabla {\boldsymbol{U}}$为$\nabla \otimes {\boldsymbol{U}}$的简写，即矢量的张量积运算；$(\nabla {\boldsymbol{U}}{)}^{\mathrm{T}}$为$\nabla {\boldsymbol{U}}$张量的转置。

1.2.2   流体渗流控制方程

求解渗流场需构建流动方程描述流体在孔隙和裂缝中的流动行为，假设压裂时储层中流体为微可压，满足达西渗流定律，且不考虑温度变化，流体控制方程即渗流场达西动量方程：

式中：${\boldsymbol{u}}$为渗流速度矢量（孔隙中流体在表征单元体积内流动的平均值），m/s；${\boldsymbol{K}}$为渗透率张量，在各向同性渗流中${\boldsymbol{K}}=K\boldsymbol{I}$，m²；K为渗透率，m²；$\mu$为动力黏度（由牛顿内摩擦定律定义），Pa∙s；$\nabla {\boldsymbol{p}}$为压力梯度，Pa/m；$\rho$为流体密度，kg/m³；${\boldsymbol{g}}$为重力加速度，m/s²。

不考虑重力时:

根据流体在孔隙中流动时质量守恒，可以得到连续性方程：

式中：$t$为时间，s；$\phi$为多孔介质孔隙度；${Q}_{m}$为质量源或汇，kg/s。

整合动达西定理与连续性方程可得:

1.2.3   渗流应力敏感性控制方程

多孔介质渗流场和应力场控制方程都含有共同孔隙压力变量，可建立渗流−应力全耦合。但压裂时压力边界条件的改变与流体渗流会导致地层孔压的重新分布，对于多孔弹性介质，孔隙压力及应力条件的变化会使岩体骨架发生变形，导致孔隙度变化。根据Kozeny−Carman公式，渗透率也相应地发生变化，渗流场再次改变，与孔压形成新一轮耦合。因此，需建立孔隙度、渗透率应力敏感性模型，对渗流参数实时更新，其耦合过程如图1所示。

图  1  渗流−应力耦合过程

Figure  1.  Seepage-stress coupling process

下载: 全尺寸图片 幻灯片

孔隙度、渗透率应力敏感性模型可表示为^[25]：

式中：${\varepsilon }_{\mathrm{v}}$为体积应变；${K}_{0}$为初始渗透率，m²。

1.2.4   达西渗流−应力耦合

联立流体控制、岩体控制以及孔渗应力敏感性控制方程（井周渗流场、岩体biot应力场以及应力敏感动态方程）共同构成页岩地层流−固耦合方程组，分析了未知数与方程个数，可知方程是封闭的。然后，给定相应的定边界条件和初始条件，建立了封闭的页岩储层流−固耦合方程，开展了页岩地层应力场动态展布规律研究。

2.   基于渗流−地质力学耦合的四维地应力模型建立

2.1   三维地质力学模型的建立

整理储层、井层等资料数据，重点针对某页岩油储层水平井段层位划定研究区。该层天然裂缝断层发育、有机质丰度高，易于形成压裂缝网。区域储层致密，地层基质孔隙度主要在2.0% ～ 6.0%之间，渗透率为0.006 6～0. 594 0 mD，最大水平地应力、最小水平地应力与垂向应力梯度分别为2.30，2.03和2.56 MPa/100m，垂向应力>最大水平地应力>最小水平地应力，属于正断层控制范围。

根据压裂时数值模拟的精度与运行条件要求，选定构造网格水平方向与垂向平面尺寸为25 m×25 m×5 m，以室内试验数据为校核点开展三维构造建模与单井属性数据区域化，选择指数模型的变差函数分析井区属性主方向和变程，采用序贯高斯方法对粗化数据进行井间插值，完成目的层三维精细化地质模型的建立，如图2所示。

图  2  压裂区地质力学模型与属性体三维展布

Figure  2.  Geomechanical modelling of fracturing regions with 3D spreading of property bodies

下载: 全尺寸图片 幻灯片

2.2   有限元地质力学模型的建立

2.2.1   压裂区地质构造转换及网格构建

基于压裂区精细地质力学模型进行有限元模型的构建。通过计算机图像学算法处理压裂微地震点云数据，在宏观模型上采用小定向边界盒方法，将各压裂段抽象为一个个压裂区，构建压裂缝网空间几何域（见图3）。解析地质构造为三角网格（TS格式），经自编程序转化为STL文件格式，转化流程如图3所示，将几何零件导入构建目的层构造几何模型，还原真实地层。为减少边界效应与薄板效应，外围增设地面及1～3倍长度围岩，组合形成一个更大的几何模型（见图4（a））。

图  3  微地震点云数据监测与压裂区构造的几何构建

Figure  3.  Microseismic point cloud data monitoring and geometric construction of the fracturing zone

下载: 全尺寸图片 幻灯片

图  4  压裂区有限元地质模型

Figure  4.  Finite element geological model of the fracturing zone

下载: 全尺寸图片 幻灯片

压裂过程中应力状态不断传递变化，外部围岩与内部储层网格划分相同精度易造成耦合网格数量过大，难以收敛。为平衡模拟计算效率、精度与稳定性，网格划分采用基质面—面扫掠自适应四面体网格与断层处三角网格处理，整体外部地层使用粗网格构建，压裂区与井筒网格扫掠局部加密细分（见图4（b）），实现区域针对性高效计算。

2.2.2   压裂区有限元模型属性参数场构建

基于压裂区精细化属性地质模型，进行地质网格到点的属性转换，形成带属性的点云数据，提取属性三维插值函数（如渗透率K(x,y,z)）对目的层有限元网格进行材料属性赋值，如图5所示（其中围岩孔渗参数为致密低值）。

图  5  有限元网格属性参数赋值方法

Figure  5.  Assignment method of finite element mesh property parameters

下载: 全尺寸图片 幻灯片

对断层裂缝和压裂缝网区域并未构建复杂的域与几何边界，而是采用等效材料属性形式体现，根据现场数据与经验公式等效设置为高孔隙度、渗透率和低弹性模量，最终通过插值函数将精细化属性模型转换到有限元网格赋值，重构压裂区有限元材料属性模型，如图6所示。

图  6  压裂区有限元属性正演

Figure  6.  Fracture zone finite element properties orthogonal

下载: 全尺寸图片 幻灯片

2.3   渗流−应力耦合物理场与边界条件构建

基于有限元围岩与区域材料参数、裂缝等效材料参数、初始地应力及孔隙压力，结合压裂水平井分段施工条件，设置耦合模型的力学物性参数、应力场外边界条件与渗流场内压边界条件，并在耦合场中设置孔隙度、渗透率的更新公式，在每个时间步实时更新孔隙度与渗透率，构建储层四维动态渗流−应力物理场，最终开展压裂过程中四维动态地应力渗流−应力全耦合研究。

2.3.1   原地应力场与耦合应力场边界条件

压裂前，地质构造在孔隙压力、初始有效地应力及构造应力作用下处于地层位移为0的平衡态，耦合模拟前需开展稳态地应力平衡计算还原地应力状态。地应力大小和方向由构造应力的方向、梯度和深度决定（见图7（a）），最大水平地应力方向约为NE90°。模拟时目的层位移场与应力场会发生变化，每次模拟需提取模拟结果坐标系下的应力张量场，作为外部应力添加到模型应力张量中再次平衡构造应力边界引起的模型的位移，按此步骤多次模拟直至位移为0，还原压裂前地应力与位移状态，如图7所示。最后将孔隙压力与初始地应力引入目的层作为初始地应力场边界条件，并设置构造应力边界作为耦合模拟的应力场外边界条件。

图  7  原地应力状态模拟与结果

Figure  7.  In-situ stress state simulation and results

下载: 全尺寸图片 幻灯片

三个地应力方向和大小显示，原地应力模拟结果范围与压裂前试验地应力大小误差在8%以内，应力机制为正断层类型，在断层处应力释放显示低值，误差较小，与目的层应力场分布相吻合。

2.3.2   渗流场物理场和边界条件

基于孔隙度ϕ(x,y,z)、渗透率插值函数K(x,y,z)，设置渗流场初始材料参数。边界条件采用内外边界双重约束：渗流场初始与压力外边界由三维孔隙压力插值函数p_p(x,y,z)确定，内边界条件根据压裂施工曲线设定为压裂工况下水平井筒压力（压裂井口油压p_co +压裂液液柱压力ρ_f）。其中，压裂液密度0.854 4 kg/L，动力黏度11.85 mPa∙s。压裂中井口油压动态变化（见图8），压裂缝网起裂至扩展的范围为68.5～74.1 MPa，结合井周边界面每点的高度，计算内边界面上每一点的压力形成内边界应力场，如图9所示。

图  8  压裂施工曲线

Figure  8.  Parameter diagram of fracturing construction curve

下载: 全尺寸图片 幻灯片

图  9  渗流场内外边界条件示意

Figure  9.  Schematic representation of the internal and external boundary conditions of the seepage field

下载: 全尺寸图片 幻灯片

3.   压裂四维地应力动态演化规律

3.1   压裂区地层孔隙压力变化特征

运用四维应力场模拟方法对目的层做切片剖面，分析压裂前后孔隙压力动态演化过程。压裂时受压裂液泵入与压裂缝网延伸的影响，在井筒内外压力梯度作用下，压裂液与压力沿高孔隙度、高渗透率的压裂缝网与断层裂缝逐渐向周围地层扩散传递，形成不稳定渗流过程，孔隙压力动态增加，与未波及地层对比形成“应力集中漏斗”，并随着压裂时间增长，应力集中漏斗的半径与压力波及范围不断扩大，如图10所示。从图10可以看出，水平压裂井周地层压力增幅最大。

图  10  压裂前后地层孔隙压力四维展布

Figure  10.  4D spreading of formation pore pressure before and after fracturing

下载: 全尺寸图片 幻灯片

结合压裂施工曲线，在井周压裂段地层取一小块区域对不同时间段的数据点进行分析，见图11。由于压裂各阶段物性参数、水力裂缝位置及扩展情况不同，孔压变化有所不同。压裂初期，取点位置地层压力为47 MPa左右，随着压裂液注入井筒内经过一定时间的憋压过程，井筒内压力快速增加但未达到水力缝网起裂点，孔隙压力缓慢增加，压力未向周围地层扩散，但井筒内壁压力与周围地层形成较大应力差。到达起裂点时岩石发生破裂形成压裂缝网，并持续扩展延伸联通天然裂缝断层，压裂液大量渗流进周围地层，地层能量得到补充，压力快速向周围传递，孔隙压力迅速增至78.3 MPa。随着压裂时间增长，裂缝被压裂液逐渐填充，压力传递至基岩，逐渐形成动态平衡，上升趋势趋于平稳，压裂末期孔隙压力稳定至88.5 MPa左右。

图  11  压裂过程中井周地层孔隙压力变化曲线

Figure  11.  Pore pressure change curve of the formation around the well during the fracturing process

下载: 全尺寸图片 幻灯片

3.2   压裂区地层地应力场变化特征

压裂区三地主应力与孔压变化趋势负相关。根据有效应力原理，孔隙压力上升导致有效应力场变化，地应力场相应减弱，同时孔隙压力变化越剧烈，孔隙压力梯度越大，应力变化越明显。上覆地层压力（垂向地应力）主要受压实作用影响，其他各项应力变化对其影响较小，因此仅展示水平两向主应力切片剖面变化趋势，如图12、图13所示。在压裂初期，起裂形成缝网后，受孔隙压力影响，水平主应力从井轴开始向两侧呈现快速递减趋势。随压裂时间延长，井筒和周围地层形成“压降漏斗”并逐渐向周围地层扩散波及，范围不断扩大，直至压裂中后期下降趋势趋于平稳，其中各压裂井近井筒区域的应力下降幅度最明显，如图14所示。截至压裂末期，最大水平主应力由45.5 MPa降至27.1 MPa；最小水平主应力由36.3 MPa降至13.6 MPa，垂向地应力与最大水平主应力的下降幅度均小于最小水平主应力。这是由于，该区域为正断层，孔隙压力变化会造成最小水平地应力比例上的更大下降，且根据裂缝扩展力学机制与岩石破裂遵循最小阻力路径的原则，当压裂液注入地层时，压裂裂缝更倾向于在垂直最小水平主应力的平面内张开（即沿${\sigma }_{\mathrm{H}}$方向延伸），形成主裂缝系统，这一过程中裂缝面张开会导致局部弹性应变能释放，裂缝壁面附近${\sigma }_{\mathrm{h}}$方向弹性应力场产生“应力卸载效应”，形成低应力屏蔽区。随着水平地应力进一步减小，岩石破裂在水平最小地应力方向阻力最小，压裂裂缝不断沿最小水平应力的方向扩展，压裂液沿优势通道快速渗流扩展孔隙压力扰动范围，加剧了${\sigma }_{\mathrm{h}}$的下降幅度。

图  12  压裂区最小水平地应力四维展布

Figure  12.  4D spreading of horizontal minimum principal stresses in fractured zones

下载: 全尺寸图片 幻灯片

图  13  压裂区最大水平主应力四维展布

Figure  13.  4D spreading of horizontal maximum principal stresses in fractured zones

下载: 全尺寸图片 幻灯片

图  14  压裂过程中井周地应力变化曲线

Figure  14.  Geostress variation curve in the area around the well during the fracturing process

下载: 全尺寸图片 幻灯片

地应力数值下降的同时，地应力方向也同样出现了变化，如图15所示。从图15可以看出，当“压降漏斗”向井筒周围扩散，遇到断层、裂缝发育及压裂缝网区域时，地应力的方向也会受影响发生偏转。这是由于，有效应力减小，岩石膨胀变形，地应力方向受挤压发生偏转，且随着孔压剧烈变化，方向偏转越明显。

图  15  压裂区地应力方向四维展布

Figure  15.  4D spread of the geostress direction in the fracture zone

下载: 全尺寸图片 幻灯片

4.   结论与认识

1）综合考虑页岩层理非均质性、地质力学特征、压裂缝网扩展复杂性等因素，提出了页岩储层压裂动态地应力渗流−应力全耦合演化模拟方法，用以指导现场压裂施工、二次压裂优化设计等。

2）地应力的变化集中分布在各生产井周区域，同时受到压裂缝网扩展与范围的影响。压裂初期，井周压裂缝网的高孔渗使井筒压力快速向周围地层传递扩散，近井区孔隙压力急剧升高，逐渐与周围地层形成“应力集中漏斗”；有效应力与三向地应力出现不同程度的减小趋势，形成相对称的“压降漏斗”，二者同步沿着压裂缝网不断扩大波及范围。至压裂中后期，裂缝被压裂液逐渐填充，压力传递到基岩，孔压与地应力场变化逐渐趋于平稳，其中最小水平主应力下降幅度最大。

3）在地应力大小改变的同时，地应力方向在裂缝断层附近会有偏转现象。