页岩气储层四维地应力演化及加密井复杂裂缝扩展研究进展

朱海燕，宋宇家，唐煊赫

1 成都理工大学油气藏地质及开发工程国家重点实验室，成都 610059

2 西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500

0 引言

受美国页岩气革命的启发和影响，我国在页岩气勘探与开发上取得了阶段性成果，使我国油气资源结构，勘探开发理念和技术发生重大变化。自2009年，我国钻探了第一口页岩气调查井渝页1井开始，我国页岩气勘探开发工作快速推进，截至2018年底，我国页岩气探明地质储量约1.17×1012m3，页岩气产量108×108m3，约占全国天然气产量的6.74%[1-3]。与此同时，我国页岩气开发也逐步形成一套较为完善的方法体系，该体系主要采用“水平井+体积压裂”的开发工艺，配合“井工厂”的开发模式，一次性部署多口水平井，集中施工，集中投产[4-6]。我国页岩气储量大，产量增速快，但仍存在地质储藏环境复杂，开发程度低，气井递减快，开发难度大等问题[1-2,7]。因此，目前我国在继续勘探新储量的同时，还需保证现有开发区域产能建设，通过部署加密井或对老井进行重复改造来弥补产量递减，提高储层动用程度[8-10]。

地层岩石力学性质、应力状态、天然裂缝发育、非均质性及各向异性等特征对水力裂缝扩展影响至关重要[11-15]。在页岩气开采过程中，气井周围地层压力下降，老井周围储层存在不同程度的压降漏斗，储层压力下降扰动原地层应力状态，而储层应力随开采时间不断演化(即四维地应力)，导致加密井压裂裂缝及老井重复压裂裂缝扩展规律发生明显变化，从而影响气井压裂和气藏开发效果。如何准确计算页岩气藏长期开采过程中储层四维动态地应力场，并以此优化加密井和重复压裂井的压裂设计，是当前我国页岩气开发领域亟待解决的难题，同时也是目前国内外研究的热点和难点。

因此，本文详细论述了油气藏渗流—地质力学多物理场耦合方法，解释了页岩气开采过程中地层属性变化规律及其对页岩气井复杂裂缝扩展的影响，从而深化对页岩气储层四维动态应力场演化的认识，提高我国页岩气储层的压裂改造效果，加快页岩气藏开发进程。

1 油气藏渗流—地质力学多物理场耦合数值模拟方法研究

油气开发经验及实验研究表明，在油气藏开发过程中，油气开采、流体注入等措施会引起地层多孔介质中孔隙压力变化，导致储层岩石发生形变，而这种变化又反过来影响地层多孔介质中流体渗流规律。通常我们将流体与多孔介质之间的相互作用称为渗流—地质力学耦合作用，或流—固耦合作用，而建立准确的数学模型和合理的求解方式是研究这一问题的核心和关键。

1.1 油气藏渗流—地质力学耦合求解形式

流体和多孔介质的相互作用研究起始于Terzahgi提出的著名的有效应力理论，认为有效应力变化会引起多孔弹性介质变形，而有效应力与其总应力和孔隙压力有关[16]。Biot[17]最早在Terzaghi的一维流动—应力耦合理论基础上提出了三维多孔介质弹性力学模型，而后Geertsma[18]提出孔隙和岩石体积变化理论，并讨论了地应力变化对岩石弹性形变和孔隙体积的影响。岩石孔隙压力减小，有效应力增大，岩石和孔隙体积压缩使岩石渗流性降低。1952年，Fatt等[19-20]最早通过实验证实了应力变化对岩石渗透率的影响，反映了应力变化会对岩石中流体渗流产生影响。1964年，Benham和Hoyle[21]系统总结了热应力的理论模型，认为温度变化会影响介质应力变化，该理论为流—固—热耦合研究提供理论依据。1977年，Brownell等[22]从平衡方程角度系统总结了多孔介质中多相流体渗流与热力学间相互作用的数学模型。1986年，Lewis等[23]撰写出版了第一部针对多孔介质变形的专著，其中首次将多孔介质的渗流数学模型和应力应变数学模型及其求解进行了总结，并长期对总结结果进行更新。

从20世纪90年代开始，国内外许多学者在此基础上对多物理场油气藏渗流—地质力学耦合模型进行完善，用于研究储层压实、地面沉降[24-25]、注水驱替引起的地层形变和天然裂缝开启[26]等问题。随着SAGD、CO2埋存、地热能开采、水合物开采等技术的出现与大量应用，需要在渗流—地质力学耦合的基础上，考虑温度对流体性质、多孔介质应力形变的影响，以及注入流体与地层岩石、流体发生化学反应等，在此基础上建立流—固—热三场耦合[27-32]和流—固—热—化四场耦合模型[33-34]，但其耦合的核心依然在于渗流场和应力场之间的耦合。

目前，渗流—地质力学耦合的数值方法按照耦合求解形式可划分为全耦合、顺序耦合和单向耦合3种，此外，还存在一种拟耦合形式。上述4种耦合求解形式的特征及差别如表1所示：

表1 不同耦合求解形式对比Table 1 Comparison of coupling solutions

(1)全耦合模型

Lewis等[24]和Gutierrez[25]最早提出了全耦合模型理念，该模型通过某一离散方法，同步求解渗流公式和地质力学公式中所涉及到的变量。常用的离散方法主要包括有限差分法和有限元法。其中有限差分法离散算法相对简单，计算效率较高，其更适用于区域较为规则的线性问题[35-36]。而有限元法则相对复杂，可用于解决非线性问题。同时，该方法单元体划分较为灵活，对于复杂岩石形变计算较为适用[37-38]。虽然全耦合能够保证应力平衡方程和流动质量平衡方程的协调性，但是由于渗流方程和应力方程同步求解时需要对大型矩阵进行计算，不仅收敛性较差，而且计算时间也相当长[39]。因此，全耦合分析一般适用于较小尺度模型的分析。

(2)顺序耦合模型

顺序耦合也称迭代耦合，Settari等[40-41]为解决全耦合模型收敛性差等问题，提出顺序迭代耦合有限元模型。该模型分为渗流模型和地质力学模型，两者需要分开建立并求解。由于渗流模型必须要在时间域上求解，而地质力学模型则只需要关注不同边界条件变化对应力平衡的影响，进而计算应力形变即可，因此，在求解过程中，需要根据研究需求将渗流模型的时间域划分成多个时间段，每一段作为一个时间步，在一个时间步内，首先求解渗流模型，并将渗流模型的计算结果(孔隙压力、温度等)传递给地质力学模型进行平衡步分析，计算本时间步内固相介质的形变及应力变化，结合应力敏感性模型，更新孔渗参数并将其传递给渗流模型，进行下一时间步的计算。

由于渗流模型和地质力学模型各自独立计算，因此只要保证2个模型分别收敛就能得到计算结果，同时由于不需要进行全耦合中的大型矩阵计算，其计算效率大大提高。经过对比分析发现，在同一网格系统内，顺序迭代耦合可达到与全耦合模拟结果相近的效果，这证实了顺序耦合模型也可具备良好的计算精度[39-40]。另外，Tortike等[42]在其基础上，考虑了热力影响。Chin等[43]提出平行力学求解模型，对迭代耦合求解方法进行优化，进一步提升了耦合计算精度。

(3)单向耦合模型

单向耦合也称显式耦合，在渗流—地质力学耦合数值建模方法发展的初期，Fung[27]、Koutsabeloulis[44]、Minkoff等[45]就提出单向耦合模型。相对于全耦合和顺序耦合，单向耦合以其求解简便迅速的优势被广泛关注。该模型计算方法是在渗流模型中完成整个分析时间内的求解后，将计算结果作为地质力学模型在不同时间点的平衡条件，进行应力平衡求解，从而得到对应时间点介质的应力形变情况。由于单向耦合不考虑固相介质应力形变对渗流参数的影响，因此其计算精度相对较差。目前单向耦合多用于井筒失稳分析及储层地质力学环境的变化趋势预测等。

(4)拟耦合模型

早期在求解渗流—地质力学耦合问题时，为了提高计算效率，将孔隙压缩性简化为与孔隙压力有关的函数，在渗流求解过程中根据孔隙压力变化即可计算孔渗参数变化和岩石变形情况，而不需进行复杂的力学计算[46]。由于该模型准确性严重依赖于孔隙压缩性经验公式，因此，对于均质和各向同性地层，其耦合模拟计算结果较为准确[47]。

不同的耦合求解方法，其迭代计算方式不同，因此各自计算效率、计算精度和适用性会有所差异[39]。通常全耦合方法被认为计算精度最高，但计算效率最低，迭代耦合、单向耦合和拟耦合方法的计算精度依次降低，而计算效率依次升高。虽然在某些模型中，或某些优化方法作用下，迭代耦合，甚至于拟耦合，能同时保证较高的计算精度和计算效率，但针对实际问题，还需根据模型特征、计算效率和精度需求等综合考量应该使用哪种算法[29,48]。

1.2 油气藏渗流—地质力学耦合求解方法

近年来，商业求解器逐渐在渗流—地质力学耦合模拟中发挥出强大的优势。相比于自编程序需要同时兼顾渗流模型和地质力学模型所付出的巨大工作量，商业求解器能够提供适用于大部分研究需求的渗流模型及地质力学模型，这极大简化了建模流程，提高了分析效率。对于特殊储层或研究需求，现有商业求解器模型无法满足要求时，只需针对储层的特征对商业求解器进行二次开发，写入自己提供的数学模型，并利用商业求解器进行求解计算即可。

目前应用较广、认可度较高的耦合系统主要有 以CodeBright[49]、DYNAFLOW[50]、FEHM[51]和COMSOL[52]等为代表的全耦合模拟器，以TOUGH和FLAC3D耦 合[53]、ECLIPSE或CMG与ABAQUS耦合[54]、ATHOS和VISAGE耦合[55]等为代表的顺序耦合模拟器组合，还有以ECLIPSE与VISAGE耦合[56]、NUFT与LDEC耦合[57]为代表的单向耦合模拟器组合，另外CMG软件平台中融合了岩石压缩模块，可用于进行拟耦合计算。表2为不同渗流—地质力学耦合模拟器组合的数学模型结构、离散求解方法及适用性等对比情况，可以看出，各求解器组合之间存在较大差异，可根据计算需求选择合适的耦合求解方法。

然而，目前的渗流—地质力学耦合模型仍存在诸多不足：(1)针对裂缝性储层，主流渗流模型中需要增加一套网格系统对裂缝进行等效处理，同时由于大型求解器中的渗流模型采用正六面体网格，这种等效方式会使得天然裂缝的尺寸及形态描述不够准确。而以离散元或边界元为代表的耦合方式虽然能够准确模拟裂缝形态，但是却无法考虑复杂基质渗流；(2)面对复杂的耦合问题，顺序迭代耦合可以提供灵活高效的求解思路，但其一般需要依赖两套以上求解器分别处理渗流模型和地质力学模型，这势必存在两套模型之间的参数传递与数据耦合，难以避免由于网格算法而产生的误差；(3)受各自数学模型和计算方法限制，现有耦合模型多基于对实际情况的简化和假设，这势必产生诸多误差和矛盾，因此，需要针对不同的实际问题，优化耦合模型方法。

2 页岩气藏渗流与地质力学耦合的四维地应力演化研究进展

页岩气井生产会产生地层孔隙压力变化，引起四维地应力演化。页岩气储层天然裂缝发育，地层流体运动和岩石变形机理复杂，开采过程中气藏渗流和岩石变形的耦合模拟还需考虑基质、微裂缝及宏观裂缝的多介质、多尺度和多维度问题[58-60]。

2.1 页岩储层四维地应力演化模拟方法研究进展

页岩储层四维地应力演化是一个复杂的多孔介质渗流—地质力学耦合过程，其源于Terzahgi提出的著名的有效应力理论，而Biot[61]及Nur等[62]在此基础上提出了有效应力系数，被称为Biot有效应力理论(式1)，并建立三维固结理论[63]，奠定了岩石力学基础，一直沿用至今。

在有效应力理论中，biot有效应力系数的确定极为重要。Geertsma[18]和Skempton[64]根据介质线弹性形变理论，提出一种广泛认同的计算方法，如式2所示。并且他们针对该公式提出了相应的测量方法，即：在岩心三轴压缩试验中，当孔隙压力为0时，即可测得岩石体积模量；而当孔隙压力与围压相同时，即可测得岩性中固体颗粒体积模型。Todd和Simmons[65]对有效应力系数的计算方法进行了扩展，如式3所示，这很大程度上提高了有效应力系数测试的灵活性，其中Q可以用纵波波速等参数取代。随后有学者研究发现用纵、横波波速及计算得到的体积模量、泊松比等同样可以计算得到有效应力系数[66]，这被称作动态有效应力系数[67]。同时，部分学者认为渗透率也可作为有效应力系数的计算参数[68-69]。而通过测量体积应变方式计算得到的有效应力系数，被称为静态有效应力系数[70-72]。对于有效应力系数取值，一直颇有争议，而目前普遍认为，其变化范围为0～1之间，且动态和静态有效应力系数基本保持一致[73-74]。

其中：Kbulk为岩心体积模量，GPa；Kgrain为岩心固体颗粒体积模量，GPa；Q为任意检测的物理量；P为孔隙压力，MPa；σeff为Terzahgi有效应力，即围压与孔隙压力之差，MPa。

有效应力系数变化受多种因素影响，包括岩石组分，孔隙结构、孔隙度等，其具有明显的各向异性[75-79]。而页岩的有效应力系数通常小于砂岩及碳酸盐岩，其取值范围约0.5～0.8，但在不同地区和储层中，其取值具有较大差异[80-82]。同时，有效应力系数还受环境因素影响，通常孔隙压力增大、围压减小、温度减低会使其有效应力系数增大[70,74,81-82]。因此，在页岩气藏生产开发过程中，地层状态变化对有效应力系数的影响不容忽视。但该规律尚未得到一个普适性公式，因此目前在页岩储层渗流—地质力学耦合模拟过程中，通常将其视为常数。另外，岩体在高围压或高孔隙压力下表现为弹塑性特征，部分学者正在进行该方面的研究和探索，但目前尚未取得统一性认识[83-84]。

页岩地层天然裂缝发育，传统连续性模型已不能满足计算需求，需要在模拟方法上加以改进，而目前改进模型主要分为连续介质模型和离散裂缝模型[58-60]。

连续介质模型是将孔隙介质分为基质孔隙介质和裂缝孔隙介质，并针对两种介质分别构建平衡方程和连续性方程等，设置两种介质之间的物质及能量传导关系。1996年Lewis和Ghafouri[85]最早提出基于有限元方法的双孔模型，对天然裂缝进行了等效处理。Jalali和Dusseault[86]对双孔模型中基质和裂缝孔隙介质应力模型进行完善，将两者有效结合。Wei和Zhang[87]在此基础上，提出适用于煤层的三孔双渗有限体积模型，并考虑吸附引起的形变影响。但由于天然裂缝展布发育差异较大，该类模型无法准确模拟天然裂缝非均匀变化。

而离散裂缝模型是将岩体处理为非连续介质，用特殊单元表征宏观裂缝，允许流体同时在基质和裂缝中流动，分别计算基质和裂缝非线性变形情况[88]。该模型数值模拟方法主要包括离散元法[11,89]、扩展有限元法[90-91]和有限元法[92-93]。其中离散元方法是在宏观离散裂缝基础上，对基质网格进行非结构化网格剖分，其计算量大，计算精度难以保证，同时离散元模型无法准确模拟基质中流体复杂流动情况[94]。而有限元方法则是在原有连续网格基础上，引入cohesive黏聚力区域模型，从而实现离散裂缝渗流和形变效果[95]，该方法计算准确度高，但依然需要进行非结构网格划分，计算效率较低。扩展有限元则不需要把裂缝当作为几何实体，该方法中裂缝面网格与结构内部的几何形状或物理界面无关，从而减小网格划分难度，但该方法主要用于二维空间，三维模型依然在探索[90,96]。

上述模型多采用全耦合方法，较适用于均质模型，但对于页岩非均质性储层，采用迭代耦合方法计算时效性较高。因此，本课题组针对页岩储层的复杂特征，提出了一种有限差分渗流模拟器和有限元地质力学模拟器相结合的交叉迭代耦合模拟方法(图1)，研究生产过程中地应力变化及其对孔渗参数的影响[97-99]。该方法中可加入离散裂缝模型，模拟天然裂缝和老井水力裂缝展布，真实反映储层实际非均质性和各向异性特征，裂缝属性可通过Oda法进行等效处理，转化为连续介质模型属性，参与迭代耦合计算[100]。

图1 渗流—地质力学迭代耦合模拟方法[99]Fig.1 Modeling workflow of flow and geomechanical iterative coupling[99]

迭代耦合模型中不同求解器之间网格数据结构存在差异(图2a)，为解决该问题，Samier等[101]提出相同网格条件下不同求解模型间数据转换方法。但由于地层结构复杂，不同求解模型网格划分方式多样(图2(b)、(c))，为了保证各模拟器的网格划分灵活性和计算准确性，课题组自主开发的网络映射子程序，通过球形自适应搜索算法(图3)，最大程度满足不同模拟器之间网格数据转化需求[97]。

图2 有限差分模型和有限元模型网格对比[97]Fig.2 The grid differences between FD model and FE model[97]

图3 有限差分模型和有限元模型网格间数据转换[97]Fig.3 Data communicating between the FD grid in the flow model and the FE grid in the geomechanical model[97]

页岩气藏流体运移机理复杂，储层中同时存在基质孔隙、微裂缝及宏观裂缝等多种孔隙介质系统，现有模型对此做了大量的简化。连续介质模型注重于宏观尺度上基质孔隙和微裂缝渗流—地质力学耦合模型，而忽略的宏观裂缝影响。而离散裂缝模型则重点研究宏观裂缝及基质，或微观尺度上微裂缝及基质之间的渗流及形变规律，忽略了不同尺度下微裂缝、宏观裂缝的非均匀变形规律及其对储层力学状态的影响。

2.2 页岩储层四维地应力演化规律研究

页岩储层渗透率低、气井产量递减较快，国内外已开始重新调整开发方案，优化井网，挖掘储层产能。其中最为有效的方式为老井重复压裂和部署加密井，而地层应力状态严重影响压裂改造裂缝扩展，因此，页岩储层四维地应力演化规律对页岩气藏开发至关重要。

Guo等[102-103]、Sangnimnuan等[104]、Gupta等[105]采用有限元方法、Roussel等[106]采用有限元与边界元方法、Safari等[107]采用边界元和有限差分方法模拟研究了页岩气藏简单两翼水力裂缝水平井生产过程中地应力变化情况。研究结果表明：开发速率越快，地层孔隙压力和最小水平主应力下降越快。由于页岩储层渗透率较低，应力变化主要集中在水力裂缝附近。初始水平两向应力越小，两井之间(水力裂缝尖端位置)越易发生水平两向地应力反转。当井间距离减小时，孔隙压力及地应力下降速率加快。同时，水力裂缝形态也会对地应力分布产生影响(图4)，而杨氏模量差异对地应力影响较小。另外，水平井空间位置、水力裂缝空间延展等也会影响地应力大小和方向分布变化。相邻两井布井层位相同时，更容易引起井间水平应力转向。而对位于不同层位两口邻井，水力裂缝高度和长度变化引起的井间干扰会使空间应力变化更加复杂。

图4 不同水力裂缝形态下孔隙压力及地应力分布变化[103]Fig.4 Pore pressure and stress distribution with different fracture geometries after production[103]

Ren等[108]和Moradi等[109]结合离散裂缝模型，考虑了天然裂缝的影响。研究表明：天然裂缝可形成高渗流通道，使其周围的孔隙压力快速降低。同时，天然裂缝也会影响改变岩体形变规律，如图5所示，横向裂缝会阻隔岩体形变趋势，而使岩体沿裂缝面发生错动。

图5 注采过程中地层位移变化分布[108]Fig.5 Formation displacement distribution after production and injection[108]

本课题组结合有限差分和有限元方法，综合考虑页岩气储层的非均质性和地质力学参数的各向异性特征，以及老井复杂水力裂缝系统(图6)，研究了四川盆地涪陵页岩气藏生产过程中四维地应力的演化规律[97,100]。研究表明：受气井生产影响，随着开采程度的提高，地层孔隙压力减小，有效应力增大，三向主应力减小(图7(a))，且主要集中在水力裂缝改造区域，而垂向两向应力差及水平两向应力差则增大(图7(b))。该区域最大主应力方向发生轻微偏转，且主要集中在井筒趾部和跟部的水力裂缝边缘，井间区域应力方向无明显偏转，该规律与Hwang等[110]及Rod等[111]研究结果相同。另外，对于非均质页岩地层，页岩水平主应力、压实程度、孔渗条件等与孔隙压力变化呈明显的非线性关系，图8为不同地层条件下，孔隙压力与地层压实程度的变化关系。

图6 水平井水力裂缝模拟结果[100]Fig.6 Simulated complex fractures of horizontal wells

图7 生产后地层应力分布与变化[100]Fig.7 Stress distribution and change after production[100]

图8 地层压实程度与孔隙压力关系[97]Fig.8 Relationship between the reservoir compaction and pore pressure[97]

页岩储层地质特征及压裂工艺复杂，在四维地应力演化模拟过程中需考虑地层天然裂缝、非均质性和各向异性特征，以及老井水力裂缝影响，因此，目前模型还需逐渐完善。生产过程中页岩储层应力与孔隙压力变化规律已基本形成统一认识，但由于模拟方法及初始地质工程条件方面存在差异，在水平两向主应力分布变化规律方面仍存在分歧，特别是井间区域水平两向主应力是否存在转向可能等问题。因此，在此方面还需深入研究探讨。

3 页岩气加密井复杂裂缝扩展研究

3.1 油气藏渗流—地质力学耦合的复杂裂缝扩展模型研究

地层力学状态变化直接影响水力裂缝扩展形态，因此需要在常规静态条件下的裂缝扩展模型基础上结合油气藏渗流—地质力学耦合模型。1980年，Settari[112]和Hagoort等[113]首次考虑油藏渗流引起的储层孔隙压力和地应力变化，建立了油藏渗流—地质力学与压裂裂缝扩展耦合的数值模型。随后，Ji等[114]建立了油藏渗流与地质力学全耦合的裂缝扩展有限元模型，实现均质油藏注采诱导地应力场动态变化条件下的单条平面裂缝扩展模拟。

近年来，由于非常规油气开发的需要，许多学者开展了加密井裂缝扩展的数值模拟研究。2012年，Gupta等[105]首次将老井生产效应与加密井水力压裂裂缝扩展相结合，利用有限元模型分析了老井生产引起的地应力变化对加密井压裂裂缝扩展的影响。随后，Roussel等[106]、Kumar和 Ghassemi[9]、Rezaei等[115]结合有限元与位移不连续方法(DDM)、Safari等[10]结合有限差分和DDM方法，分析了均质地层条件下，老井生产过程中的地应力变化，以及在其变化下加密井的水力裂缝转向和非均匀扩展，但模型中加密井水力裂缝均为简单两翼裂缝。Zhang等[11]在离散元模型中引入离散裂缝网络(DFN)，研究了老井生产对加密井复杂水力裂缝扩展的影响，并系统讨论了裂缝性储层中水力裂缝起裂和扩展，以及天然裂缝与水力裂缝、多条水力裂缝之间影响。为了描述真实的储层地质特征，Huang等[116]、Marongiu-Porcu等[117]通过集成“地质力学—微地震监测—气藏渗流—离散天然裂缝网络”多物理场模型，提出了一套加密井复杂裂缝扩展模拟的地质工程一体化解决方案，但其渗流—地质力学耦合模型较为简单，水力裂缝扩展为拟三维状态，无法真实表达页岩地层水力压裂复杂裂缝三维扩展。

本课题组[98,100]针对涪陵页岩气的真实储层特征及生产参数，提出了集成“离散天然裂缝网络—老井压裂复杂裂缝扩展—老井气藏渗流与应力耦合—加密井裂缝扩展”的综合建模方法(图9)，研究了页岩气藏加密井复杂裂缝扩展机理。该方法根据地震、测井、岩心观察及实验测试等资料，分别建立了精细三维地质网格属性模型、地质力学模型及离散裂缝模型，并在此基础上建立基于离散裂缝的水力压裂复杂裂缝扩展模型，采用Oda方法将离散裂缝属性等效转化为连续网格属性。通过渗流—地质力学耦合模型计算老井生产过程中地层属性变化情况，并在此基础上模拟加密井裂缝扩展情况。

图9 加密井复杂裂缝模拟流程[100]Fig.9 Complex fracture propagation simulation workflow of infill well[100]

现有加密井裂缝扩展模型主要有2种类型，一是考虑油藏渗流—地质力学全耦合的均质各向同性储层模型，该模型无法考虑复杂裂缝扩展问题；二是“地质模型—老井压裂—气藏渗流与地质力学耦合模拟—加密井复杂裂缝扩展模拟”多模型交叉迭代的非均质性和各向异性储层模型，复杂裂缝的扩展均是基于解析或简化的裂缝本构模型，这些均不能真实反映页岩储层流体运移、岩石形变及复杂裂缝扩展机理和特征，需要在连续介质模型和离散裂缝模型基础上，深入探究复杂裂缝扩展数理模型，考虑模型尺度影响，优化计算时效。

3.2 页岩气加密井压裂裂缝扩展规律研究

在页岩气藏开发过程中，受老井生产诱导应力场干扰影响，老井水力裂缝改造区边缘位置可能会产生地应力方向偏转，甚至水平两向应力反转，图10中绿色区域为模拟所得应力转向区[10,105]。由于加密井井筒位置处地应力偏转程度相对较小，同时射孔方向对水力裂缝扩展存在诱导作用，加密井水力裂缝在初期扩展时依然沿原最大水平主应力方向，当裂缝进入应力转向区内时，扩展方向逐渐发生偏转，阻碍加密井裂缝扩展进入老井压裂改造区，如图11所示[10,106,115,118]。而地层应力变化范围和变化程度与老井生产时间及产量有关，因此，从老井生产到加密井压裂之间存在一个时间窗口，在窗口期内实施加密井压裂能够取得较好的压裂效果[10]。

图10 生产后地层最大水平主应力大小及方向分布[10]Fig.10 The magnitude and direction distribution of the maximum horizontal principal stress after production[10]

图11 水平地应力差对加密井裂缝形态的影响[10]Fig.11 Effect of horizontal stress difference on fracture propagation path[10]

页岩气井的生产会使其压裂改造范围内孔隙压力和三向主应力降低，加密井压裂过程中受“Frac-hit”效应影响，水力裂缝产生非对称扩展，使其更倾向于向老井孔隙压力下降区域扩展[11,119-120](图12)。这一现象已经在现场试井、示踪剂测试及压裂微地震监测结果中得到了验证[120-122]。加密井压裂“Frac-hit”效应受多方面因素影响，如：井距、老井生产程度、地层物性条件、地层应力状态、天然裂缝发育、老井压裂时间、水力裂缝复杂性、压裂液类型、施工参数等，因此，老井与加密井压裂裂缝扩展规律应该综合多方面因素进行研究[123-124]。

图12 不同加密井压裂时间水力裂缝扩展模拟结果[119]Fig.12 Simulation results of hydraulic fracture propagation of infill well in different fracturing timing[119]

页岩储层层理和天然裂缝发育，页岩基质脆性较强，水力压裂过程中，水力裂缝受其诱导改变扩展方向，产生复杂裂缝网络。周健和陈勉等[125-126]通过真三轴水力压裂物理模拟实验验证天然裂缝对水力裂缝扩展方向的影响。杨春和团队[127-128]通过物理模拟实验发现水力压裂过程中页岩层理会使水力裂缝等人工裂缝发生止裂、分叉、穿过或转向等现象(图13)，最终演变成5种模式：顺从、先顺从后转向、贯穿闭合、贯穿开启及多裂缝，如图14所示。邹雨时等[129]研究了不同天然裂缝分布对水力裂缝扩展的影响，郭天魁等[130]和Tan等[131]研究了不同应力状态、注入速率、液体黏度下，水力裂缝在多天然弱面或层理页岩中的扩展规律。实验结果表明：低应力差、高注入速率、较小的液体黏度有利于沟通天然弱面或层理，形成复杂裂缝网络。因此，在页岩气压裂模拟中，单一裂缝并不能准确表征水力裂缝展布。另外，对于老井开发对新井水力裂缝扩展影响的实验研究较少，Bruno[132]通过实验研究表明了岩石局部注水会改变其应力状态，使水力裂缝发生偏转(图15)，而老井生产作用对加密井水力裂缝扩展影响目前还未得到实验验证。

图13 裂缝网络横切面示意图[128]Fig.13 Schematic section of fracture network[128]

图14 水力裂缝与天然裂缝相互作用过程[128]Fig.14 Interaction process between hydraulic fracture and natural fracture[128]

图15 不同注水条件下裂缝扩展形态对比[132]Fig.15 Comparison of fracture propagation geometry under different injection conditions[132]

3.3 页岩气加密井四维地应力演化的微地震屏障效应

本课题组[98,100]针对涪陵页岩气的真实储层特征及生产参数，采用“离散天然裂缝网络—老井压裂复杂裂缝扩展—老井气藏渗流与应力耦合—加密井裂缝扩展”综合模拟方法，初步探索了某开发平台在特定井间距条件下的加密井复杂裂缝扩展规律。研究发现：老井生产会引起其水力裂缝网络范围内孔隙压力和三向主应力减小，水平两向应力差增大，但应力场方向并未发生明显转向，同时，加密井位置处地层力学状态变化不大。受地应力变化影响，加密井水力裂缝依然沿原水平最大主应力方向扩展，但裂缝形态与老井相比发生较大变化。加密井水力裂缝主要集中在井筒周围，特别是层理缝扩张形成的分支扩张裂缝；越靠近老井，加密井裂缝数量越少，且主要为构造缝扩张形成的分支裂缝，层理缝难以开启；而老井水力压裂时，地层保持初始力学状态，其裂缝扩展较为均匀(图16)。这使得相同施工液量下，加密井水力裂缝扩展范围小于老井。

图16 水力压裂裂缝复杂裂缝模拟结果[100]Fig.16 Hydraulic fracturing simulation results[100]

基于此，本课题组首次将该现象定义为加密井“微地震事件屏障效应”[98,100]。如图17所示，微地震事件在加密井井周附近分布情况与老井压裂施工类似，但靠近前期生产井的微地震事件震级和数量将发生明显的减弱，甚至未检测到微地震事件，加密井在前期生产井改造区域内未形成复杂裂缝。在该区域上存在类似一道屏障将微地震事件阻挡在外，因此将这一现象称之为微地震事件“屏障效应”。与微地震事件“屏障效应”相对应，靠近前期生产井的区域称之为微地震事件“屏障区”。

图17 微地震事件屏障效应示意图[100]Fig.17 Schematic of the microseismic events barrier eff ect[100]

不仅如此，北美Eagle Ford和Bakken页岩气区块的现场压裂施工微地震反演结果中观测到了类似的现象[133]，图18为一开采了14年的某页岩井组加密井压裂微地震事件反演结果，该井组所在区块储层厚度 20～60 m，渗透率 1×10-4～5×10-4mD，孔隙度6%～8%，含水饱和度0.3，初始水平两向主应力分别为56 MPa和55 MPa，初始孔隙压力48 MPa，岩石弹性模量15 GPa，且天然裂缝发育，在物性及地质力学特性方面与涪陵页岩气藏具有一定的相似性。在旁边仅有一口老井压裂并生产的情况下，实施加密井压裂。从压裂后的微地震事件点分布状态可以看出，微地震事件点从加密井井筒向老井逐渐靠近，当到达老井生产压降区(压裂改造区)时，微地震事件点数量突然急剧减少，与涪陵发现的现象相似。同时，由于老井生产区域内地应力整体下降，水力裂缝更容易在该区域内扩展，从而导致了加密井裂缝非对称现象，该现象恰是本研究所提出的“加密井微地震屏障效应”的一种特殊表现形式。

图18 加密井水力压裂微地震事件分布[133]Fig.18 Hydraulic fracturing microseismic events of infill wells[133]

4 页岩气藏地应力演化及调整井复杂裂缝扩展发展方向展望

(1)页岩气藏水平井重复压裂和加密井压裂复杂裂缝扩展模拟

随着国内主要页岩气区块进入中期开发阶段，原开发区域内老井的产能迅速下降，除勘探新的区块或向深部层位下探外，老井重复压裂及部署加密井是提高采收率的重要措施。然而，由于老井经历了数年的开采，储层地质力学特征发生了一定程度的变化。在此环境下，无法准确认识其裂缝扩展机理，就有可能导致加密井在压裂过程中与老井压窜，或造成重复压裂裂缝无效扩展，影响重复压裂及加密井改造效果，甚至产生老井套管变形等安全事故。

储层压力、地应力等地质力学参数随着页岩气开采不断演化，致使储层条件非均匀变化，这是加密井压裂裂缝扩展模拟与老井初次压裂的最大不同。同时，老井重复压裂时，初次压裂裂缝影响也不容忽视的。因此，准确预测老井重复压裂及加密井复杂裂缝扩展，需要弄清页岩气藏长期开采条件下裂缝性储层地质力学参数的动态演化。

目前尚未对页岩重复压裂及加密井压裂裂缝扩展机理形成较为统一的方向性认识。现有油气藏渗流—地质力学耦合的裂缝扩展模型尚不能准确反映页岩地层天然裂缝发育、非均质性及各向异性特征，不同尺度下基质、微裂缝及宏观裂缝的流体运移及地质力学属性变化机理模型尚不完善，对该条件下复杂裂缝扩展问题未形成一套较为完善的数值模型和模拟方法，因此，针对这些问题还需进一步探索和研究。

(2)页岩气储层多层立体开发中复杂裂缝空间干扰机理研究

目前我国主要页岩气田一期已初步完成大量的加密井部署(同一储层平面加密布井)，今后几年，龙马溪组上部页岩气藏储层将成为开发的重点，因此，我国页岩气将面临下部储层和上部储层同时开发(简称立体开发)的情况[6,134]。在该开发模式下，下部储层开发改变了下部储层的孔隙压力、地应力，对上部储层的压裂将产生额外的诱导干扰。该问题涉及上下部储层精细地质建模、下部储层初次压裂、加密井压裂、气藏模拟等复杂的地质工程一体化交叉融合。如何实现上、下部储层充分改造，但不发生垂向压窜，是主要的技术瓶颈。

然而，国内外对该问题的研究较少，目前的模型及研究成果主要针对于单一均质油气藏模型，开展同层位内老井生产过程中地层属性变化及加密井裂缝扩展机理研究，尚未考虑老井生产对上下多地层属性动态影响，及其对不同层位加密井裂缝扩展规律的影响。因此亟待开展该方面研究，为我国页岩气的立体开发提供理论与技术指导。

(3)页岩气藏水平井重复压裂和加密井压裂时机优化

页岩气重复压裂及加密井压裂方案制定需以其产量和经济效益为评价标准，丰富的储量、良好的流体运移条件、最优的裂缝扩展是保证其改造效果的必然条件。由于老井改造及生产开发作用，地层条件变化难以直接观测，因此，准确模拟和预测地层属性变化成为压裂改造的核心和前提条件。在此基础上，选择最佳压裂时间可为重复压裂及加密井压裂效果提供有力保障。但该优选过程复杂，涉及多因素、多尺度、多维度的综合影响，而目前研究仅局限于其中某一环节，忽视了压裂生产全过程的综合考量；仅局限于部分静态参数，忽视了时间尺度的深入探讨。

因此，页岩气藏水平井重复压裂和加密井压裂需要同时考虑渗流、地质力学及裂缝扩展，结合地质和工程各因素，建立一套“地质精细解释模型+地质力学模型+水力压裂裂缝扩展模型+产能预测模型+经济评价模型”的综合模拟评价方法，从时间和空间维度上研究页岩气藏开发过程中各环节地层变化情况，为页岩气开发方案优选提供理论支撑。

(4)基于地质工程一体化的页岩气藏子母井压裂干扰及套管损伤机理

近年来随着页岩气藏体积压裂改造技术的广泛应用，压裂过程中的套管变形和损伤问题日益突出，其直接导致桥塞无法坐封到位或套管破损，影响后续层段施工，缩短页岩气井生命周期，严重制约页岩气藏开发。

水平井多级压裂过程中套管变形及损伤机理复杂，除考虑压裂过程中液体流动、岩石变形、裂缝及断层滑移等地层力学状态变化，还需研究该条件下固井水泥和套管等受力形变问题。另外，目前页岩气开发普遍采用“井工厂”开发模式，多井同步施工，段间、井间干扰不容忽视，地层属性变化存在累加效应，因此，该过程涉及多物理场、多介质、多维度耦合模拟分析。而目前模型大多针对于单一裂缝或断层滑移条件下套管及水泥环受力变形分析，未考虑复杂地质条件及多井、多压裂段复杂裂缝扩展的相互影响。对该问题，特别是套变风险点位置准确预测尚未形成统一认识[135-136]。

地质工程一体化研究集地层构造及特征分析、地质力学分析、井位部署、钻井、固井、压裂设计与施工分析于一体，以三维数值模型为载体，精细分析地层岩石力学特征，计算各工程施工过程中井筒及地层状态演化，为准确预测套变风险点和有效解决套管变形损伤问题提供理论依据。

5 结论

(1)目前油气藏渗流—地质力学耦合模型种类多样，商业求解器的开发与应用为流体运移和地质力学演化等问题研究提供了便利，但各方法的求解模型、计算时效及适用性不同，我们可根据研究需求进行选择及优化。

(2)页岩气藏地质特征复杂，储层非均质性强、天然裂缝发育且需水力压裂后投产，目前储层渗流—地质力学耦合模型均在不同程度做了假设和简化，不能准确反映地层基质、微裂缝及宏观裂缝多尺度条件下流体运移及岩体变形机理。

(3)页岩渗流—地质力学耦合的裂缝扩展模型需要考虑储层非均质性、天然裂缝等因素之外，还需要考虑气藏长期开发过程中储层地应力、孔隙压力等的动态非均匀演化，目前模型只针对其中某些方面进行详细分析，尚未形成较为完善的模拟方法，因此，目前需要综合考虑上述因素，开展相关研究，进一步明确页岩加密井及重复压裂井复杂裂缝扩展机理。

(4)页岩开发过程中的气藏渗流—地质力学耦合及裂缝扩展研究是多物理场、多维度、多尺度的耦合问题，需要探索地质工程一体化的解决方法，在此基础上，建议下一步开展复杂裂缝扩展研究、立体开发空间干扰研究、压裂时机及参数优化、子母井压裂干扰、套管变形及损伤机理研究，对我国页岩气藏的持续高效开发具有重要意义。