0 引言
1 数学模型
1.1 模拟过程
燃爆-水力复合压裂数值模拟过程包含燃爆应力冲击造缝、燃爆气体扩缝与水力压裂3个阶段。在模拟时，整个模拟过程中使用同一套初始网格与边界条件。首先进行燃爆应力冲击阶段的模拟，对该阶段形成的裂缝形态进行统计，将其作为下一阶段的初始裂缝，以此实现阶段之间的衔接。气体扩缝阶段的初始压力取应力冲击时的平衡压力，经过两个阶段的模拟形成燃爆裂缝。水力裂缝与上一阶段的衔接过程采用同样方式，导入燃爆裂缝作为初始裂缝，最终实现燃爆-水力复合压裂数值模拟。
1.2 CDEM基本原理
1.3 固体本构模型
采用线弹性-断裂本构模型描述岩石基质的力学行为，通过增量法来表征线弹性本构[18]：
对于岩石基质中的块体单元，均满足应力平衡方程，应力平衡方程的单元矩阵形式为[19]：
由几何关系得到节点应变：
由应力应变关系获得节点应力：
每个单元节点的加速度根据该节点在前一个时间步长的不平衡力计算。然后，采用显式欧拉差分法迭代求解下一时间步长的位移场，迭代格式为：
1.4 界面模型及裂缝扩展准则
相邻有限元单元之间通过接触连接，相邻单元接触点的相对位移与弹簧力的关系满足胡克定律[20]：
当两个单元之间的法向与切向任一弹簧出现断裂，说明裂缝起裂。用于判断界面断裂的条件如下：
最大拉应力准则，
Mohr-Coulomb准则，
当条件满足（7）式时，弹簧发生拉伸破坏，此时弹簧力修正为Fn=0；当条件满足（8）式时，弹簧发生剪切破坏，此时弹簧力修正为 。当弹簧发生断裂之后，所产生裂缝的开度等于弹簧两端节点的位移差：
1.5 渗流计算模型
通过块体单元之间的裂隙单元来描述流动，认为裂隙单元是具有一定宽度的平行板，其中流体流动满足立方定律。气体和液体的渗流场基本控制方程和求解方法相同，气体渗流时假设为理想气体且不考虑气体的黏度变化，液体渗流时假设其为微可压缩流体。
裂隙渗流的基本微分方程式为[21]：
裂隙单元节点的不平衡流量由流经裂隙单元节点的流量和施加在裂隙节点的源汇项计算得到[22]：
1.6 爆源模型
应力波冲击阶段：由于燃爆的时间非常短暂，尺度为毫秒级别，因此忽略燃爆过程中的热交换，将其看作绝热膨胀过程。压力变化使用朗道与斯坦纽科维奇公式计算[23]：
气体扩缝阶段：在燃爆应力波冲击作用之后，在爆源空间内仍存在大量高压气体，将燃爆高压气体的压力作为压力边界条件施加到井筒周围裂隙单元节点处可计算气体向地层中的流动。忽略爆炸空间内气体温度的变化，由理想气体状态方程可知：
Δt后爆源空间内气体压力为：
1.7 模型耦合关系
燃爆压裂模拟过程涉及固体本构模型、界面模型、渗流模型与爆源模型，水力压裂模拟过程涉及固体本构模型、界面模型、渗流模型，各模型相互耦合。
2 模型求解与验证
2.1 模型求解
CDEM采用动态松弛技术进行显式迭代求解，该方法通过在动态计算中引入阻尼项，使得初始不平衡的振动系统逐渐衰减到平衡位置。该方法的基本计算流程如图3所示。
2.2 模型验证
2.2.1 与物模实验对比验证
2.2.2 气体流动模型验证
2.2.3 与KGD（Khristianovic-Geertsma-de Klerk）模型理论解对比验证
通过与物模实验结果、气体流动模型模拟结果及KGD理论解的对比，充分验证了本文所建模型的可靠性。
3 燃爆-水力复合压裂数值模拟
3.1 直井压裂地质模型
根据上述数学模型，建立尺寸20 m×20 m的小尺寸直井压裂二维模型，开展燃爆压裂造缝后的地应力分布模拟，分析水力裂缝起裂扩展机理。同时建立尺寸为100 m×40 m的大尺寸直井压裂二维模型，进行燃爆-水力复合压裂模拟，分析燃爆-水力裂缝的展布规律与影响因素。两个模型均在中间设置一直径为0.16 m的圆孔作为井筒，同时也作为爆源。模型网格为三角形网格，在井眼附近划分小尺寸网格，由井眼向外网格尺寸逐渐增大。模型基本参数参考2 000~3 000 m实际页岩气储层数据，爆源参数参考甲烷气体数据（见表1）。
表1 燃爆-水力复合压裂模型参数表 |
| 类别 | 参数 | 数值 |
|---|---|---|
| 岩石基质 | 密度 | 2 600 kg/m3 |
| 抗拉强度 | 5 MPa | |
| 弹性模量 | 40 GPa | |
| 泊松比 | 0.2 | |
| 内摩擦角 | 30° | |
| 内聚力 | 10 MPa | |
| 原始地层压力 | 30 MPa | |
| 渗透率 | 0.1×10−3 μm2 | |
| 地应力 | 最大水平主应力 | 48 MPa |
| 最小水平主应力 | 40 MPa | |
| 爆源参数 | 爆速 | 400 m/s |
| 爆热 | 5×107 J/kg | |
| 爆源密度 | 50 kg/m3 | |
| 泵注参数 | 密度 | 1 000 kg/m3 |
| 黏度 | 1 mPa·s | |
| 排量 | 3 m3/min |
3.2 燃爆压裂后地应力分布
图7为采用小尺寸模型模拟燃爆压裂后，水力压裂初期压裂液完全进入燃爆裂缝时地层最小水平主应力与水平应力差的分布。由图7a可看出，在裂缝延伸方向与σH方向接近的裂缝的尖端附近区域，最小水平主应力降低幅度最大，故水力压裂时该区域的燃爆裂缝会率先起裂、扩展。图7b显示，裂缝延伸方向与σh方向接近的裂缝附近，主要为低应力差分布区域；裂缝延伸方向与σH方向接近的裂缝附近，主要为高应力差分布区域。说明在水力压裂初期，延伸方向接近σH方向的燃爆裂缝会沿着σH方向扩展，而与σH方向存在一定偏离角度的燃爆裂缝，会先沿着初始的方向进行扩展，随后再逐渐向σH方向偏转，这说明在燃爆压裂造缝后再进行水力压裂，可以同步开启多方位的燃爆裂缝并延伸，实现更充分的压裂改造。
3.3 水力-燃爆复合压裂缝展布规律
4 影响因素分析
为定量描述燃爆-水力复合压裂后储层的破裂程度，采用σH、σh方向的最大破裂长度与改造面积指标评价储层改造效果。同时，使用σH、σh方向的最大破裂长度与破裂度评价燃爆压裂裂缝形态。最大破裂长度定义为某方向上裂缝长度的最大值；改造面积定义为σH方向最大破裂长度、σh方向最大破裂长度分别作为长轴、短轴形成的椭圆面积；破裂度定义为破裂单元面积占模型所有单元总面积的比值。
采用表1中的基础数据，对地应力差、燃爆压裂峰值压力、燃爆压裂增压速率、水力压裂排量与压裂液黏度5个参数进行单因素分析，探究不同因素对燃爆-水力复合压裂裂缝扩展的影响。为方便进行单因素分析，设定峰值压力200 MPa、增压速率720 GPa/s、排量3 m3/min、黏度1 mPa·s、地应力差为8 MPa为模拟对比基础值，单因素影响模拟时，仅改变其中一项参数值。
4.1 地应力差
改变地应力差为2，8，14，20 MPa，模拟不同地应力差下燃爆-水力复合压裂裂缝形态（见图12）。由图可知，随着地应力差增大，σh方向发生破裂的阻力增大，燃爆-水力复合裂缝更倾向于沿σH方向扩展，燃爆-水力复合裂缝在σH方向破裂长度增大，在σh方向破裂长度明显减小，燃爆-水力复合裂缝整体形态变窄。
由图13可知，随着地应力差的增加，燃爆-水力复合压裂改造面积明显减小，当地应力差由2 MPa增至20 MPa时，σh方向最大破裂长度由23.8 m降至3.9 m，σH方向最大破裂长度由46.5 m增至70.5 m，改造面积由869 m2降至215 m2。燃爆压裂对地应力差的敏感性较小，随地应力差的增加，最大破裂长度的变化趋势与燃爆-水力复合压裂基本一致，但变化幅度要小得多。分析原因，燃爆压裂的持续时间短，瞬间能量高，能在极短的作用时间内压制地应力的控制，故地应力差的影响较小。实施燃爆-水力复合压裂，地应力差较小时，燃爆压裂可在σh方向形成较多的裂缝，破裂度较高，水力压裂时，水力裂缝沿燃爆裂缝末端扩展延伸，且不同方向裂缝扩展时不易相互沟通，最终形成多条沿不同方向延伸的复合压裂裂缝；地应力差较大时，燃爆裂缝在σh方向上的破裂长度与地层破裂度均减小，水力裂缝在σh方向上的破裂长度受限，在地应力作用下，水力裂缝在近井地带迅速偏转，最终聚并成一条近似沿最大主应力方向扩展的主缝，大幅降低燃爆-水力复合压裂的改造面积。
4.2 燃爆压裂峰值压力
改变燃爆压裂峰值压力为90，200，310，420 MPa，模拟峰值压力对压裂效果的影响（见图14、图15）。由图可知，随峰值压力的增大，燃爆压裂在σh方向上的裂缝长度明显增加，σH方向上的裂缝条数也明显增加。说明峰值压力越高，燃爆压裂在近井区各个方向产生的径向裂缝越长、微裂缝越多，更有利于后续水力压裂形成复杂缝网。低峰值压力时，由于燃爆压裂在近井区形成的裂缝数量相对较少、裂缝长度相对较小，燃爆-水力复合压裂仅形成沿最大主应力方向延伸的一条主缝，效果较差；高峰值压力时，燃爆压裂在近井区形成的裂缝数量多、裂缝长，水力裂缝沿燃爆裂缝扩展延伸，虽然同样受地应力差的影响形成偏转，但在σh方向上仍有较大的改造范围，储层改造效果较好。
由图15可知，燃爆压裂过程中，随着峰值压力的增大，在σH、σh方向上造成的最大破裂长度、破裂度均有不同程度增加，这有利于后续水力压裂的改造。燃爆-水力复合压裂过程中，随着峰值压力的增大，σh方向最大破裂长度随之增大，峰值压力从90 MPa增至420 MPa，σh方向最大破裂长度从4.4 m增至16.7 m，增幅279%；峰值压力对σH方向最大破裂长度的影响甚微。随着峰值压力的增大，改造面积增加，峰值压力从90 MPa增至420 MPa，改造面积从235 m2增至881 m2，增幅275%。
4.3 燃爆压裂增压速率
改变燃爆压裂增压速率为120，420，720，1 020 GPa/s，模拟增压速率对压裂效果的影响（见图16、图17）。可以看出，燃爆压裂过程中，增压速率增加，燃爆裂缝在σH方向和σh方向上的最大破裂长度均小幅增加；从改造效果来看，随增压速率的增加，燃爆压裂的破裂度增幅较大，说明增压速率越大，燃爆压裂形成缝网越复杂，水力裂缝更容易沿不同方向延伸，获得更大的改造范围。燃爆-水力复合压裂过程中，随着增压速率的增加，复合裂缝在σH方向的最大破裂长度有所增加，但当增压速率超过420 GPa/s后，σH方向的最大破裂长度反而有所降低，而σh方向的最大破裂长度呈小幅增加的趋势；从形成裂缝的复杂程度来看，增压速率越高，裂缝改造面积越大，缝网越复杂。从图17b的统计数据看，增压速率从120 GPa/s升至1 020 GPa/s，裂缝在σH方向上的最大破裂长度减少2.1 m，σh方向的最大破裂长度增加2.2 m；改造面积增加99 m2，增幅21%，增加幅度相对有限，影响程度远低于峰值压力。
4.4 水力压裂排量
改变水力压裂排量为1，3，6，9 m3/min，模拟排量对压裂效果的影响（见图18、图19）。由图可知，增大水力压裂排量，σh方向的最大破裂长度增大，σH方向的最大破裂长度减小。排量增大，裂缝转向半径增大，同时高缝内净压力可促使不易起裂的燃爆裂缝继续扩展，裂缝沿不同方位扩展更均衡。排量为1 m3/min时，生成4条长裂缝，而排量提高为9 m3/min时，生成长裂缝12条，这说明大排量下更多的燃爆裂缝在水力压裂液的作用下继续扩展，产生更多数量的长裂缝。与6 m3/min排量相比，虽然9 m3/min排量下裂缝沿各方向扩展更均衡，σh方向的最大破裂长度有所增大，平均缝宽也更大，但σH方向的最大破裂长度减小，改造面积有所下降，说明以改造面积为优化目标时，排量存在最优值，该条件下最优排量为6 m3/min。
4.5 水力压裂液黏度
改变水力压裂液黏度为1，5，15，25 mPa·s，模拟压裂液黏度对压裂效果的影响（见图20、图21）。由图可知，黏度增大，σh方向的最大破裂长度增加，σH方向的最大破裂长度减小；黏度增大，缝内净压力提升，裂缝转向半径增大，有利于裂缝沿原方位均衡扩展延伸，储层改造面积增大，但黏度与改造面积并非呈完全正相关关系，当黏度超过15 mPa·s时，改造面积不再继续增加。分析原因认为，在较高缝内净压力的作用下，裂缝宽度显著增大，如黏度为15 mPa·s时，平均缝宽为1.45 mm，当黏度增大到25 mPa·s时，平均缝宽增加至1.95 mm，在相同的模拟时间下，压裂液的液量是恒定的，因此裂缝宽度的增加影响了裂缝的延伸距离。
5 结论
燃爆-水力复合压裂可有效结合燃爆压裂近井造缝复杂和水力压裂远井深穿透的优势，形成多条深穿透的长裂缝，大幅提升压裂改造效果，优于单一燃爆压裂或水力压裂。
地应力差增大，燃爆-水力复合压裂改造面积减少，地应力差较小的储层，燃爆-水力复合压裂改造效果更好；高峰值压力、高增压速率均有利于增加燃爆裂缝的最大破裂长度与破裂度，进而增加燃爆-水力复合压裂的改造面积，提高改造效果，其中峰值压力的影响更显著。
提高水力压裂排量、压裂液黏度会提升缝内净压力，利于激活沿各个方向分布的燃爆裂缝，增大裂缝转向半径，产生更多的长裂缝，有效增大储层改造面积，但储层改造面积与排量、黏度并非完全正相关，存在临界值，当超过临界值后，改造面积减小。
符号注释：
A——弹簧面积，m2；c——黏聚力，Pa；C——单元阻尼矩阵，kg/s；Cf——流体压缩系数，Pa−1；D——燃爆速度，m/s；D——弹性矩阵，Pa；E——弹性模量，Pa；F——单元节点外力列阵，N；Fn，Ft——弹簧上的法向力和切向力，N；G——剪切模量，Pa；h——地层厚度，m；i，j——分量编号；k——裂隙单元编号；K——体积模量，Pa；K——单元刚度矩阵，kg/s2；Kn，Kt——法向、切向刚度，N/m；Ks——地层渗透率，m2；l——模型长度，m；L——界面等效厚度，m；M——摩尔质量，kg/mol；M——单元质量矩阵，kg；n——迭代步数；n12——接触面的单位法向量，无因次；Ne——与节点相连的裂隙单元数；p——流体压力，Pa；p0——平均爆轰压力，Pa；pa——大气压力，Pa；pb——爆生气体的实时压力，Pa；pc——临界压力，Pa；pf0——初始爆源空间内气体压力，Pa；pf1——经Δt时间后爆源空间内气体压力，Pa；pg——气体压力，Pa；pi——原始地层压力，Pa；pin——入口压力，Pa；pout——出口压力，Pa；q——裂隙节点流量，m3/s；qin——单元流入流量，m3/s；Q——源、汇，s−1；Qf——流量，m3/s；Qw——单位质量物质的爆热，J/kg；Qw0——井筒周围裂隙单元流量之和，m3/s；r——径向距离，m；R——气体常数，J/（mol·K）；t——时间，s；Δt——间隔时间，s；T——抗拉强度，Pa；Tf——地层温度，K；Tfa——地层平均温度，K； Tg——气体温度，K；u——位移矩阵，m；u1，u2——弹簧两端节点位移，m； ——节点速度列阵，m/s； ——节点加速度列阵，m/s2；Δun——法向位移，m；Δut——切向位移，m；V——单元体积，m3；V0——气体初始体积，m3；Vb——爆生气体体积，m3；Vc——临界体积，m3；Vf——爆源空间体积，m3；w——裂缝宽度，m；x——从起点沿渗流方向到计算点的距离，m；X——横坐标方向的长度，m；Y——纵坐标方向的长度，m； ——气体压缩因子，无因次；γ——高压绝热指数，无因次；γ1——低压绝热指数，无因次；δij——克罗内克记号，无因次；ε——节点应变，无因次；Δεij——应变张量增量，无因次；η——气体导压系数，m2/s；Δθ——体应变增量，无因次；μ——流体黏度，Pa·s；ρg——气体密度，kg/m3；ρw——爆源空间内的装药密度，kg/m3；σ——节点应力，Pa；σh——最小水平主应力，MPa；σH——最大水平主应力，MPa；σij——应力张量，Pa；Δσij——应力张量增量，Pa；σij（t0），σij（t1）——t0与t1时刻的应力张量，Pa；σn——法向应力，Pa；σt——切向应力，Pa；φ——内摩擦角，（°）。