页岩油立体开发诱发的储层地应力动态响应特征

郭旭洋 金衍 林伯韬

1.中国石油大学(北京)油气资源与探测国家重点实验室；2.中国石油大学(北京)石油工程学院

在我国原油对外依存度逐年攀升的背景下，页岩油大规模建产已经成为保证能源安全的一项重要手段，鄂尔多斯盆地、准噶尔盆地等都已成为我国页岩油的主要产区［1］。由于页岩基质渗透率低，自然产能无法满足工业化生产需求，常通过水平井和水力压裂等手段进行商业化开采［2-4］。Eagle Ford、Bakken 等北美地区的页岩油开采已有一些成熟经验［5-6］，而我国页岩油具有分布面积和规模相对较小、但累计厚度大的特征。因此，北美页岩油开采更倾向于在同一埋深内进行横向上的大规模开发，而我国页岩油开采更倾向于在水平面上较小的范围内充分进行立体开发。页岩油立体开发的关键是尽可能增加压裂缝网与低渗储层的接触面积，增加流体的流动能力。立体开发需要在纵向上进行多层水平井布井［7］，使得垂向上水平井井间距较小，投产后会造成储层压降，进而通过压力波诱发储层地应力的复杂演化。储层地应力对立体开发中投产井的邻近井位的钻完井、压裂都有直接影响。因此，压降诱发的储层地应力响应特征就成为页岩油高效立体开发需要研究的一个关键问题。

金衍等［8］、邓金根等［9］认为综合水压致裂法、Kaiser 效应法、差应变法、多级子测井方法可以较好确定页岩层的原始地应力。在此基础上，常可通过建立数学模型的方法对开采过程中的地应力动态演化进行分析。由Biot 等［10］创立并发展的孔隙弹性力学理论是描述储层压降诱发的地应力演化这一物理过程的基础，该理论阐述了一维至三维空间下饱和多孔介质中总应力、有效应力和孔隙压力的关系。据此，科研人员针对页岩特征不断完善相关方法和理论，研究了低渗基质中的流体流动特征、吸附作用、扩散现象、应力敏感性、孔隙结构等相关因素的作用，考虑了基于动量平衡的页岩骨架变形效应［11-13］。Gupta 等［14］通过水平面内地应力演化计算结果确定了最优井工厂加密井井距，并指导了重复压裂。郭旭洋等［15］认为，生产压降诱导的储层地应力演化可能导致井工厂的水平井压裂出现复杂裂缝形态，进而造成裂缝串通，干扰邻近水平井生产，影响Eagle Ford 页岩油井工厂的开发效率。Roussel 等［16］则提出了一种有限元流固耦合模型，表征储层二维平面内生产压差作用下的孔压、地应力变化规律，并通过水力压裂模拟证明了储层开采诱发的地应力演化会显著改变同一井工厂内邻近水平井的裂缝密切割形态。

现有成熟的页岩层地应力表征方法仅适用于未经开发的原始状态，进入开采状态后的地应力演化过程一般通过数学模型进行计算。目前已有的数模研究主要针对水平面内横向尺度情况，较少针对立体情况下纵向延展较大的页岩油储层。然而我国页岩油储层累计厚度大且采用立体开发的情况较多，因此针对性地结合现场数据建立三维计算模型，开展立体开发过程中的储层地应力响应特征研究。

1 储层地应力响应数学模型

生产压差造成的压降诱发页岩油储层中地应力演化涉及到2 个力学问题：渗流力学和固体力学。其中，渗流力学描述水平井生产过程中地下流体在页岩中的流动过程；固体力学描述流体流动造成的岩石骨架变形过程。传统意义上的研究局限于单独研究压力场(渗流力学问题)或单独研究应力场(固体力学问题)，因而无法表征生产诱发地应力演化这一压力场、应力场协同作用的物理过程。

Biot 孔隙弹性理论［11］指出，多孔介质的孔隙压力和应力状态可以通过Biot 系数和有效应力概念联系起来，建立了应力作用下岩石变形与流体压力的关系为

式中，σij′为有效应力，Pa；σij为总应力，Pa；α为Biot 系数；δij为狄拉克德尔塔函数；p为孔隙压力，Pa；Kd为岩石排水模量，Pa；Ks为岩石固相体积模量，Pa。式（1）和式（2）说明了孔隙压力和应力的求解是对模型展开研究的关键参数。

为了表征储层开采造成的压力场演化，基于质量守恒建立的流动方程为［17］

式中，mi为i相(油或水)在体积元内的运移质量，kg/m3；fi为i相在面积元上的流动质量，kg/(m2· s)；n为单位法向向量；t为时间，s；qi为水平井设计产量，kg/(m3· s)；Ω为储层内方程的域；Γ为方程求解区域的边界。式(3)表征了质量守恒过程中质量累计项、质量流动项、汇/源项的平衡。流动项可进一步表示为

式中，ρi为i相密度，kg/m3；k为渗透率，μm2；μi为i相黏度，Pa· s ；p为孔隙压力，Pa。

为了表征页岩岩石变形造成的应力场演化，采用准静态假设，根据动量平衡建立应力平衡方程为

式中，σ为应力张量；ρb为固体密度，kg/m3。

式(5)描述了三维空间内各方向是主应力和切应力的分量与外界作用的应力平衡状态。根据无穷小变换假设，位移和应变的关系表示为

式中，ε为应变；∇s为对称梯度；u为位移，m。

得到孔隙介质流动和岩石骨架应力平衡方程后，需要通过Biot 孔隙弹性理论将2 个物理场进行流固耦合，具体通过孔隙度变化实现

式中，φ为孔隙度；εv为体积应变。

综合式(3)至式(7)则可以求解压降诱发的三维地应力演化。对时间采用有限差分方法进行隐式欧拉法求解，以确保数值解的准确性和收敛性，对空间离散采用伽辽金方法(Galerkin Method)进行有限元分析，并采用Newton-Raphson 法处理系统的非线性特征。页岩油储层水力裂缝在三维模型中采用局部网格加密(Local Grid Refinement)和储层改造体积(Stimulated Reservoir Volume)的方法进行表征。

2 三维建模

F 油田页岩油储层厚度较大，适合立体开发手段开采。研究区域的流体和储层物性见表1。前期已投产分段分簇水力压裂水平井的某段裂缝微地震监测数据如图1 所示，绿色线段代表水平井段，紫色代表共计13 个微地震事件及其位置，灰色轮廓代表该段压裂的储层改造监测区域(缝网)，作为对单段压裂规模的测算，对应的裂缝几何形态见表2。重点研究储层三维空间内的孔压和地应力演化。

表1 模型流体、储层参数Table 1 Model fluid and reservoir parameters

图1 某段压裂裂缝微地震数据Fig.1 Microseismic data of one certain hydraulic fracture

表2 单段裂缝几何形态Table 2 Geometry of single fracture

根据上述储层、裂缝数据，建立如图2 所示三维数值计算模型，研究水平井一段裂缝投产后的储层压降和地应力演化，重点研究压力场、应力场在纵向随时间发生立体演化的过程。三维建模中，通过2 条离散裂缝及SRV(储层改造区域15 m×180 m×15 m)，对图1 中单段水力压裂裂缝进行表征。在纵向上，研究的储层厚度为缝网高度的3 倍左右，以表征地应力响应在纵向尺度上的变化特征。SRV 位于模型正中心，裂缝高度在模型中z方向位置是15 m至30 m。模拟缝网生产时，使用固定井底流压5 MPa，生产时长为2 年。通过模拟，可以获得三维空间内孔压、地应力的演化特征，进而进行分析。

图2 考虑裂缝的三维计算模型Fig.2 3D calculation model considering fractures

3 计算结果与分析

缝网生产2 年后三维模型不同深度层的孔隙压力分布情况如图3 所示，5 个深度层裂缝均未穿透。第1 层(a)深度为0～5 m，孔隙压力分布为12.3～13.0 MPa；第2 层(b)深度为5～10 m，孔隙压力分布为11.8～13.0 MPa；第3 层(c)深度为10～15 m，孔隙压力分布为9.8～13.0 MPa；第4 层(d)深度为15～20 m，孔隙压力分布为5.1～13.0 MPa；第5 层(e)深度为20～25 m，孔隙压力分布为5.0～13.0 MPa。(a)～(e)层对应的深度为三维模型中的z方向对应的高度，而不是埋深。

每个深度层平面内初始最大主应力方向与裂缝方向相同，生产诱发储层压降后引起的最大主应力方向发生变化，与裂缝方向产生的夹角用 θp表示为

式中，θp为主应力方向变化角度，rad；τxy为剪应力，Pa；σx为x方向主应力，Pa；σy为y方向主应力，Pa。σx、σy和 τxy均因为缝网生产而发生变化，各个深度层内主应力角度也会因而发生变化。

图3 不同深度层孔隙压力分布Fig.3 Pore pressure distribution in different layers

图3 中，远离SRV 的区域，孔压变化几乎可以忽略，这是由基质极低的渗透率阻碍了生产井压降向储层远处传播而造成的。这些区域的最大主应力转向也几乎可以忽略，均为指向初始最大主应力方向，这是由于压降的压力波几乎未传播到这些区域。各层的SRV 均可观测到孔压和最大主应力的角度变化。(c)～(e)层由于被水力裂缝穿过，造成的压降和主应力转向较为明显，而且这种趋势越接近水力裂缝中间层越明显，并在(e)层达到最明显。(a)、(b)层由于未被水力裂缝穿过，虽然可以观测到一定的压降和最大主应力转向，但是与(c)～(e)相比较为轻微。需要注意的是，根据理想均质假设，主应力转向数值计算结果均以裂缝为中心呈左右对称，但由于绘图采集的样本点不以裂缝左右对称，使得图中主应力转向没有严格按裂缝左右对称。

图3 中各个深度层最大主应力 σH的分布情况如图4 所示。受立体开发影响，σH的变化主要集中在SRV 区域以及SRV 沿y方向延伸数10 m 的区域。σH演化最明显的区域是SRV 中的(e)层，计算得到(e)层的 σH最大值为23 MPa。类似地，如图5 所示是各个深度层最小主应力 σh的分布情况。σh的变化也是主要集中在SRV 区域内，在(e)层 σh最低低至20 MPa。在SRV 区域外，在x方向和y方向均有σh的演化。(a)～(b)层由于未被水力裂缝穿过，压降不明显，但(c)～(e)层压降造成的应力扰动较明显地传递到(a)～(b)层，因此造成少部分区域最小主应力在生产2 年后上升的现象。

图4 生产2 年后各个深度层最大主应力分布Fig.4 Distribution of the maximum principal stress in each formation at different depth after 2 years’ production

图5 生产2 年后各个深度层最小主应力分布Fig.5 Distribution of the minimum principal stress in each formation at different depth after 2 years’ production

如图6 所示为SRV 内和SRV 外 σx和 σy随时间的变化，总体上呈现随时间递减的规律，这主要是由于孔压下降造成的：由于总应力由有效应力和孔压组成，孔压下降会造成总应力显著下降。但SRV 外σx在(a)～(c)层初始阶段呈现上升趋势，这是由于初始阶段这些层位的压降诱发的有效应力上升无法被有限的压降抵消。此外，越接近SRV 中间层 σx和σy随时间的演化就越明显，这也可以说明SRV 对立体开发造成的储层应力动态演化具有增强作用。

图6 σx 和σy 时间变化曲线Fig.6 Time dependent curves of σx and σy

图3 至图6 的结果表明，立体开发造成的储层地应力大小和方向的演化与纵向深度有显著的关系。SRV 穿过的层位中地应力变化的响应更为明显，SRV 未穿过的层位中虽然地应力大小和方向的变化有限，但是依然造成了一定的地应力场改变，仍可能影响到立体开发井网中纵向上邻近井位的钻完井和压裂效果。

4 层间应力干扰分析

立体开发过程中水平井分段分簇裂缝造成的不同层位中的应力干扰是页岩油现场开发关注的一个重点，可以通过研究纵向尺度上的水平最小主应力σh和地层破裂压力pf来表征。破裂压力表达式为

式中，pf为破裂压力，Pa；σh为最小主应力，Pa；σH为最大主应力，Pa；St为岩石抗拉强度，Pa。在层间应力干扰分析中，重点研究经过2 年立体开发后相关变量与初始值对比的变化量Δσh和 Δpf。

如图7 所示为水平井所处x-z垂直平面内的Δσh分布。立体空间内，该垂直平面与水平井井筒在同一平面内。结果显示分段分簇裂缝生产导致的层间最小主应力干扰可以延伸至SRV 的厚度以外，直至储层的顶端和底端，并在这些区域造成最高可达1 MPa 的应力值增加。在SRV 内，由于压降较为明显，导致最小主应力值可降低最多3 MPa。立体维度下，这种层间干扰可以扩散至SRV 以外15 m 处。这说明如果在图示水平井位置的深浅±20 m (总厚度40 m)以内布水平井，均可能受到图示水平井造成的应力干扰。

图7 x-z垂面内不同层位间最小主应力扰动情况Fig.7 Disturbance of minimum principal stress between different horizons in the x-z vertical plane

如图8 展示，水平井所处x-z垂直平面内的Δpf分布。图示水平井对立体层位造成的地层破裂压力干扰主要集中在与SRV 同宽度的范围内，在厚度上则波及至储层顶端和底端，造成了地层破裂压力增加。这意味着在图示水平井深浅±20 m 且与SRV 同宽度的区域内进行立体开发时，造成地层拉伸破坏所需的井内流体压力升高，且最高可升高4 MPa，使地层更难发生破裂。此外，在x方向上，靠近SRV 左右边缘较近(＜30 m)的区域内，则由于干扰作用导致地层破裂压力下降，最多可下降1 MPa，说明此区域地层拉伸破坏所需流压降低，更易发生地层破裂。

图8 x-z垂面内不同层位间地层破裂压力扰动情况Fig.8 Disturbance of formation fracturing pressure between different horizons in the x-z vertical plane

5 结论

(1)投产后，在纵向上会影响3 倍SRV 厚度内储层的地应力大小和方向。该影响在SRV 中间层最大、在SRV 远端最小。SRV 对立体开发造成的储层地应力演化具有增强作用。

(2)总应力在演化过程中总体上递减，但是在压降诱发的有效应力增加大于压降的区域，可能出现暂时的总应力增加。

(3)水平井生产造成的应力干扰可波及至深浅±20 m (总厚度40 m)区域内，对最小主应力和地层破裂压力造成干扰，增加或降低地层破裂难度。