1. 引言
开采低渗岩石内油气时,常采用水力压裂技术进行储层增产改造。然而,水力压裂的应用亦伴随着对水资源的大量消耗、污染地下饮用水层、增大地震活动概率等环境风险[1]-[3]。为克服水力压裂技术存在的诸多问题,各国学者正在积极探索替代方法。无水压裂技术因其无水、无残留等特点,成为代替水基压裂液的首选方案[4]。
围绕影响气体压裂增产效果的因素,国内外学者开展了大量的试验研究。Alpern和Gan等[5] [6]对聚甲基丙烯酸甲酯材料进行了以H2O、CO2和N2为压裂流体介质的压裂试验,结果发现采用氮气压裂的破裂压力最小并且压裂产生的裂隙最复杂,试样内部形成了高比表面积的螺旋形裂隙网络。Gomaa等[7]利用低渗透率页岩研究了压裂液流体粘度和注入速率对破裂压力和裂纹扩展的影响,结果表明破裂压力和裂纹的复杂性与流体粘度呈很强的指数关系,流体粘度越大,破裂压力越大,而裂纹的复杂性随流体粘度的增大而减小;Li等[8] [9]通过理论分析与实验探究了H2O、CO2和N2压裂对页岩破裂压力和裂缝拓展形态的影响,发现与其他压裂液相比,CO2压裂会产生较为复杂的压裂模式;Cai等[10]通过室内实验测试了液氮对岩石的压裂效果。发现在相同的喷射压力下,液氮射流在磨料颗粒加速和增压方面的表现优于水射流。李畅等[11]在不同围压下对无烟煤开展水、ScCO2压裂试验,结果表明ScCO2压裂形成的裂纹数量更多,影响范围更大。
破裂压力是在外部流体作用下岩石内部产生裂纹时所对应的流体临界压力。目前对岩石破裂压力的理论主要集中在岩石破裂准则上,其中最大拉应力准则和应力强度因子法则应用最为广泛。Hubbert和Willis [12]基于最大拉应力准则,提出了预测不可渗透岩石裂纹启裂压力的H-W模型;黄荣樽[13]考虑了上覆岩层产生的垂直应力对破裂压力的影响,对H-W模型进行修正;Zhang [14]等认为H-W模型计算的是岩石启裂压力而不是破裂压力,增加了额外输入能量进行修正。Haimson和Fairhurst [15]基于均质、各向同性以及小变形的假设,考虑了压裂液滤失,给出了可渗透岩石裂纹启裂压力的H-F模型。李文魁[16]通过在井筒周围分布对称单裂缝的假设上,结合I型应力强度因子提出了岩石的破裂压力计算式;Zhang等[17] [18]综合考虑了地应力、注射压力和孔隙压力,基于应力强度因子法则,给出了破裂压力的解析解。
本文基于叠加原理建立低渗岩石气体压裂力学模型,选取氮气、空气、二氧化碳三种气体,利用MatDEM离散元软件构建了带孔的二维岩石数值模型,模拟了围压与恒定气体压力耦合作用下,岩石内部裂纹扩展、孔隙压力分布及能量演化过程,研究不同气体对低渗岩石裂纹扩展及能量演化特征的影响。
2. 裂纹扩展与能量演化理论
2.1. 裂纹扩展
裂纹扩展过程:初始时,岩石内部的裂隙很小,或处于紧闭的状态,在钻孔开始施加压力时,流体被压入岩石裂隙内部,裂隙周边的受力状态开始发生变化。
当裂隙的周围受力达到临界条件时,裂隙开始的位置,也就是裂隙尖端开始产生较大的裂隙,这个过程就是裂纹启裂的过程,接下来会导致裂纹不稳定,进一步发生裂纹的扩展。
裂纹的尖端主要受围岩压力、流体压力、孔隙压力三种力的作用,三种力叠加作用,导致裂纹的扩展。引入应力强度因子,只要裂纹尖端的应力强度达到一定的数值就会发生扩展,反之裂纹的扩展就不会发生。
根据叠加原理,可以把这个复杂的受力简化为图1所示的情况:
假设
为应力强度因子之和,
为围岩压力引起的应力强度因子,
为流体压力引起的应力强度因子,
为孔隙压力引起的应力强度因子。在塑性力学中,按照受力情况,可以分为三种不同类型的裂纹,在此次岩石裂纹扩展的实验中,裂纹的类型为I类型裂纹,即受到的是垂直于裂纹面的拉应力的作用。
围岩压力作用下:
(1)
公式中
为裂纹尖端到是岩石中心的距离,
为形状因子。形状因子与
、
有关,其中
为裂纹的长度。
流体压力作用下:
(2)
公式中,
为流体压强,
为形状因子,与裂纹的长度和岩石的尺寸有关系。
孔隙压力作用下:
(3)
公式中
为钻孔半径。综合上述公式可以得到:
(4)
如果
超过了岩石的临界应力强度因子时,裂纹开始在围岩压力、流体压力、孔隙压力的作用下扩展。根据上述公式,就可以计算出岩石开始破裂的时间和破裂时的受力情况。根据叠加原理,应力强度因子公式如下:
(5)
可以用图1形象的表示叠加原理。根据叠加原理,可以把这个复杂的受力简化为如图1所示的情况:
Figure 1. Stress simplification diagram of gas fracturing in low-permeability rocks based on superposition principle
图1. 基于叠加原理的低渗透岩石气体压裂应力简化图
2.2. 能量演化理论
岩石材料产生裂纹,发生破坏的本质是能量驱动下的岩石部分结构的改变。从岩石开始受力到产生裂隙并稳定,这个过程中能量演化贯穿始终,发挥着重要的作用。通过气体压力、围岩压力的作用,向岩石内部输入能量,最初阶段,弹性能在岩石的内部积聚,塑性能、断裂热等通过产生裂隙的方式被消耗。当强度超过岩石的强度极限时,内部积聚的弹性能被释放出来转化为岩石的动能。从能量演化的角度出发,可以更有利于研究岩石的破坏本质。
数值模拟实验要研究岩石在受到外力的情况下,能量演化的规律,即在施加气体压强、围岩压力的条件下,记录不同时间段岩石内的弹性能以及耗散能的演化规律。岩石能量的演化规律可以分为四个阶段:能量输入、能量积聚、能量耗散、能量释放。可以将能量简化为弹性能和耗散能两部分。即根据能量守恒,可以得到如下公式[19]:
(6)
为外界对岩石的能量输入,
为岩石内部的弹性能,
为岩石由于产生裂隙所耗散的能量。
MatDEM离散元软件中,系统的机械能是由动能、弹性势能、重力势能组成的。动能是颗粒运动产生的能量,计算公式如下:
(7)
弹性势能是颗粒之间切向弹簧和法向弹簧的应变能相加,计算公式如下:
(8)
重力势能是重力的存在产生的能量,计算公式如下:
(9)
在离散元系统中,动能是一个瞬态值,最终会通过阻尼振动转换为阻尼热和摩擦热[20]。系统的耗散能可认为是由阻尼热、断裂热、摩擦热相加所得。其中,阻尼热是作用在颗粒上的阻尼力消耗系统中动能所产生的;断裂热是颗粒间的弹簧发生断裂时失去的弹性势能之和;摩擦热是因颗粒间产生相互滑动而产生的能量。
3. 数值模拟流程与方案
3.1. MatDEM原理
MatDEM软件是由南京大学自主研发的一款高性能离散元软件。采用原创的矩阵离散元计算方法和三维接触算法,可以实现数百万颗粒的快速离散元模拟[21]。
离散元法是一种不连续数值模拟方法,从离散元颗粒之间的关系开始分析,分析接触模型,并建立力学模型,根据牛顿第二定律对岩石中的颗粒进行模拟仿真。
离散元法是用来讨论模拟由颗粒组成的物质的相互作用以及运动形式,比如力学、矿业领域的岩石。从宏观角度解析,岩石是连续的;但是从微观角度来看岩石是由颗粒、缝隙组成的,是非连续的,离散的。离散元可以建立颗粒之间的关系,可以更为有效的模拟岩石的非连续性,在岩土、矿业等领域发挥着很大的作用。离散元为研究岩土的力学性质,提出了一套从微观角度入手的解决方案。
离散元是通过堆积并胶结很多具有指定力学性质的颗粒来构建真实的岩石模型,可以假设颗粒之间存在一根横向弹簧,并假定弹簧的弹性系数为
,弹簧的断裂位移为
。设颗粒之间的法向力为
,法相变形为
,法向刚度为
,且
。
假设颗粒之间存在一根竖向弹簧,并假定最大剪切力为
,颗粒之间的初始抗剪力为
,颗粒之间的摩擦系数为
。设颗粒之间的剪切力为
,剪切变形为
,切向刚度为
,且
。
1) 法向力(线弹性模型):
弹簧恒定工作存在压缩
、无变形
、拉伸
三种情况:
(10)
弹簧破坏
:
弹簧破坏。
2) 剪切力(线弹性模型):
弹簧恒定工作存在压缩
、无变形
、拉伸
三种情况:
(11)
弹簧破坏
:
弹簧在切向上的破坏准则为莫尔–库仑准则:
,
为初始抗剪力,即在不受法向压力时材料能承受的最大剪切力,最大抗剪力与初始抗剪力有关。当切向力大于最大剪切力的时候,切向的弹簧断裂,此时颗粒间只有滑动摩擦力
。
经过上述的分析可以得到每个颗粒的受力情况,在第一个时间步内,通过牛顿第二定律
,单元合力与颗粒的质量之比得到每个颗粒在第一个时间步时的加速度,进而得到加速度增量,与初始速度相加得到下一个时间步开始时的初始速度,并通过上述分析,得到颗粒新的受力情况。随着时间步的推移,实现迭代的过程。
3.2. 模拟流程
模拟过程分为三个步骤:堆积建模、赋予颗粒材料属性、迭代算法进行数值模拟。
1) 建立堆积模型
为讨论中间带孔的二维正方形岩石,在围岩压力
、
,气体压力
的作用下,岩石的裂纹扩展及能量演化规律,建立如下模型进行数值模拟(图2)。
Figure 2. Gas fracturing model establishment in low-permeability rocks
图2. 低渗岩石气体压裂模型建立
第一步建立一个正方体的模型箱,之后的步骤均在模型箱内进行。并通过sampleW、sampleL、sampleH定义模型箱的长宽高、ballR等指令定义颗粒的平均半径等其他单元属性。
第二步进行堆积建模,为了使得模拟结果与真实情况更接近,在施加载荷前进行堆积建模,模拟真实情况下的重力沉积过程。岩土颗粒在重力沉积后可能会经历压实的作用,使用compactSample指令,使得上压力板产生载荷并压缩模型(图3)。
第三步挖孔,以便于气体压强的施加。通过find指令,找出以正方形中心为圆心、给定数值为半径的圆所覆盖的所有单元,并通过delElement指令删除圆内的单元(图4)。
Figure 3. MatDEM numerical model of low-permeability rock
图3. 低渗岩石MatDEM数值模型
Figure 4. Gas fracturing borehole configuration
图4. 气体压裂钻孔
2) 岩石材料设置
第四步对产生的模型进行材料赋值,如表1所示。离散元紧密堆积模型的宏观力学性质与微观力学性质之间存在关系。可以通过输入宏观力学性质:杨氏模量、泊松比、抗压强度、抗拉强度、内摩擦系数,通过MatDEM材料训练窗口得到微观力学性质:法向刚度、切向刚度、断裂位移、初始抗剪力、摩擦系数。得到围观性质后通过上述公式,进行迭代计算。宏微观转换公式[13]如下:
(12)
(13)
(14)
(15)
(16)
(17)
Table 1. Material parameters of low-permeability rock formation
表1. 低渗岩石材料参数
岩石材料参数 |
数值 |
杨氏模量 |
5 GPa |
泊松比 |
0.2 |
抗拉强度 |
2 MPa |
抗压强度 |
20 MPa |
内摩擦系数 |
0.6 |
密度 |
2600 kg/m3 |
单元平均直径 |
0.6 mm |
3) 迭代计算
第四步,施加载荷并进行数值模拟。在施加气体压强时,由于在每一个时间步后会进行平衡处理,使用单元编号加载气体压强无法恒定工作。利用find指令找到距离圆中心
距离内的颗粒,利用setBallPressure指令对颗粒周围孔隙施加气体压强(图5)。
Figure 5. Gas pressure application on peri-borehole fractures in rock particles
图5. 对钻孔周围颗粒间裂隙施加气压
3.3. 数值模拟方案
选择边长为5 cm的岩石,在正中间挖出一个半径为2.5 mm的钻孔,用于施加气体压力。设置围岩压力为0.1 MPa,气体压力为恒定的5 MPa,进行三轴压缩数值模拟。通过宏微观转换公式,输入宏观性质:弹性模量E、抗压强度C、抗拉强度T、泊松比
、密度
、内摩擦系数
转换为微观力学性质:法向刚度
、切向刚度
、断裂位移
、最大抗剪力
、摩擦系数
,材料宏观参数如表1所示。
在数值模拟过程中,记录裂隙的形态及数量、裂纹扩展特征、孔隙压力以及能量演化等数据。
在区别氮气、二氧化碳、空气时,运用范德瓦尔斯方程对参数进行调整。范德瓦尔斯方程,简称范氏方程,是荷兰物理学家范德瓦尔斯于1873年提出的一种实际气体状态方程。范氏方程是对理想气体状态方程的一种改进,特点在于将被理想气体模型所忽略的气体分子自身大小和分子之间的相互作用力考虑进来,以便更好地描述气体的宏观物理性质。
由于MatDEM软件内嵌的关于跟该不同气体的方式:是通过调节参数,再带入公式。所以将范德瓦尔斯方程拟合为一元一次方程
的形式,其中
为密度,P为气体的压强。通过拟合的结果,计算出直线的斜率K和截距C。不同气体K和C的数值见表2。
Table 2. Physical parameters of fracturing gases
表2. 压裂气体参数
不同气体 |
K |
C |
氮气 |
0.0000129397 |
0.6987042159 |
二氧化碳 |
0.0000217789 |
−0.8676968004 |
空气 |
0.000013429 |
−0.5663 |
计算出不同气体的数值后,再通过MatDEM软件p.fluid_k=K;p.fluid_c=C的命令对不同的气体进行设置。
4. 数值模拟结果分析
4.1. 不同气体压裂的裂纹扩展特征
本章计算分析了在不同气体压裂的过程中,岩石中产生的裂纹数量。考虑到裂纹产生的不同机理,计算分析了张拉裂纹、剪切裂纹和总裂纹的数量。图6(a)~(c)分别是氮气、二氧化碳、空气压裂的过程中裂纹数量随时间变化的曲线。
通过对三组曲线进行分析,可以发现在气体压裂过程中产生的张拉裂纹的数量占总裂纹数量的99%以上,几乎没有剪切破坏。不同介质对裂纹萌生、扩展速率及最终稳定长度的影响存在明显差异。在相同时间内氮气压裂产生的裂纹数量最多,最少的是空气压裂产生的裂纹数量。氮气压裂在0~0.2 ms内裂纹数量迅速增加至270条,然后持续增长至330条,裂纹扩展阶段能量释放过程较为平缓。二氧化碳压裂在0~0.2 ms时曲线变化与氮气压裂相同,迅速增加至270条,但之后裂纹数量几乎不再增加最终达到280条,裂纹产生初期气体膨胀能释放更为剧烈。空气压裂的裂纹数量在6 ms内达到205条,其曲线形态呈现两阶段特征:0~0.3 ms快速增加至155条,然后曲线近似水平,在1.2 ms时突然增加至205条,之后裂缝数量几乎不再增加。
(a) 氮气压裂的裂纹变化曲线
(b) 二氧化碳压裂的裂纹变化曲线 (c) 空气压裂的裂纹变化曲线
Figure 6. Crack change curves Under different fracturing fluids
图6. 不同压裂流体压裂过程中的裂纹变化曲线
图7(a)~(c)为不同压裂流体在压裂过程中的裂隙分布情况,可以发现氮气压裂相较二氧化碳压裂产生的裂隙数量更多,在试样钻孔区域更为集中,使得裂隙相连接形成较大规模的裂纹,并从钻孔中心扩展至右上方处。空气压裂的裂纹数量最少,仅沿钻孔至右上方发育较小规模的裂缝。
(a) 氮气压裂的裂隙分布
(b) 二氧化碳压裂的裂隙分布 (c) 空气压裂的裂隙分布
Figure 7. Fracture distribution characteristics under different fracturing fluids
图7. 不同压裂流体压裂过程中的裂隙分布
4.2. 不同气体压裂的孔隙压力特征
图8(a)~(c)为氮气、二氧化碳及空气压裂过程中孔隙压力场的空间分布特征。
氮气压裂形成的孔隙压力场呈现显著的空间局域化特征,钻孔中心附近存在高压核区(Pmax = 3.5 MPa)。高压区从钻孔中心呈椭圆形扩展,分别向右下方扩散和向上方然后压力降为2.7 MPa扩散至右上方,在高压区边界压力迅速降至0.5 MPa以下。这种集中式压力分布源于氨气的高粘度和低压缩性抑制了流体向远端扩散,促使能量优先积聚于微裂隙尖端,驱动密集张拉裂纹的持续扩展。
二氧化碳压裂在钻孔周围(半径为0.007 m)至右上方的区域内,孔隙压力分布较为均匀,大小在3.2 MPa左右,在高压区边界压力骤降至0.5 MPa以下。这种宽域低压缓变模式与二氧化碳的低粘度和高扩散系数密切相关,导致流体快速渗入微小孔隙但难以维持局部高压,削弱了裂纹扩展的持续性。
空气压裂的压力场呈现长条型分布特征,从钻孔周围(半径为0.005 m)均匀扩散至右上方,最大压力值4.7 MPa,这种均匀化压力分布使得钻孔周围缺乏大范围高压区,对应裂纹萌生数量最少,验证了空气压裂的低效性。
(a) 氮气压裂的孔隙压力
(b) 二氧化碳压裂的孔隙压力 (c) 空气压裂的孔隙压力
Figure 8. Pore pressure evolution during fracturing with different gases
图8. 不同压裂流体压裂过程中的孔隙压力
4.3. 不同气体压裂的能量演化特征
图9(a)~(c)分别展示了氮气、二氧化碳及空气压裂过程中能量随时间演化的曲线。
氮气压裂相同时间内产生的总能量最高,能量曲线呈现分阶段特征,分为初期快速生成阶段和稳定增长阶段。在0~0.2 × 10−3 ms时能量产生速率最快,达到0.4 × 10−3 J后,产生能量的速率变低,以近于恒定的速率在7.9 × 10−3 ms内从0.4 × 10−3 J增加至2.7 × 10−3 J,表明散热路径受限,易引发局部热应力集中。
二氧化碳压裂的能量曲线主要分为快速生成阶段和稳定阶段,其总能量在0.2 × 10−3 ms内迅速增长至4.8 × 10−4 J,然后逐渐增加至5.2 × 10−4 J,最终保持稳定。
空气压裂的总产能量最低,其总产能量在0~0.5 × 10−3 s内迅速攀升至2.1 J,然后维持稳定,到1.1 × 10−3 ms时迅速升高至3.2 × 10−4 J,而后产生的能量保持稳定。
(a) 氮气压裂的能量曲线
(b) 二氧化碳压裂的能量曲线 (c) 空气压裂的能量曲线
Figure 9. Heat curves during fracturing with different fracturing fluids
图9. 不同压裂流体压裂过程中的能量曲线
根据上述分析,氮气压裂过程中会产生更多的裂隙并消耗更多的能量。从流体特性角度来看,氮气的黏度高于二氧化碳。式18和式19 [22]为流体损失速率v1与损失系数k1、流体黏度μ的关系式。损失系数会随着流体黏度的增大而减小,即黏度较高的流体具有较低的流体损失率。
由于氮气的高黏度特性,在裂纹尖端区域,其流体损失量大幅降低。面对微小裂缝,氮气凭借较低的流体损失率,能够更有效地促使裂缝扩展。在加载恒定的气体压力下,氮气压裂效果总体上略优于二氧化碳压裂。这一结论不仅为优化压裂工艺提供了理论依据,同时也在非常规油气资源开发领域具有重要的实践指导价值。
(18)
(19)
5. 结论
1) 氮气压裂在低渗岩石中产生裂纹数量最多,裂纹扩展呈现两阶段特征,初期快速扩展后持续平缓增长,形成高密度螺旋形裂隙网络;二氧化碳压裂次之,裂纹萌生初期能量释放剧烈但后续扩展受限;空气压裂效果最差,裂纹扩展呈现间歇性突增特征。张拉裂纹占比超99%,表明气体压裂以拉伸破坏为主导机制,氮气的高粘度特性通过降低流体损失率有效促进微小裂隙扩展,是其优于其他气体的关键因素。
2) 氮气压裂形成局域化高压核区,压力梯度显著,利于能量积聚驱动裂纹持续扩展;二氧化碳压裂孔隙压力分布均匀但缺乏有效压力梯度,削弱了裂纹扩展动力;空气压裂呈现长条型低压分布,能量分散导致压裂效率低下。氮气的高粘度和低压缩性抑制流体扩散,维持裂隙尖端局部高压状态,是其形成复杂裂隙网络的重要物理基础。
3) 氮气压裂能量转化效率最高,体现为弹性势能持续转化为断裂能的稳定耗散过程;二氧化碳压裂能量生成呈两阶段特征,初期能量快速释放后进入平衡态;空气压裂能量演化呈现间歇性跃升,反映其不连续的裂纹扩展模式。流体粘度通过控制能量积聚速率(氮气 > 二氧化碳 > 空气)和耗散路径(氮气局域化耗散/二氧化碳全域耗散)共同决定压裂效果,为优化无水压裂工艺提供理论依据。
基金项目
国家自然科学基金项目资助(编号:52174091,52204113);国家重点研发计划项目资助(编号:2022YFC3004600)。
NOTES
*第一作者。
#通讯作者。