含三维深埋裂纹脆性岩石水力压裂数值模拟研究

2019-08-24陈礼婧

水力发电 2019年5期

陈礼婧，张敏

(1.广西理工职业技术学院土木工程(建筑材料)学院，广西崇左532200；2.江西科技师范大学建筑工程学院，江西南昌330013)

0 引言

裂纹与节理作为岩体的固有属性，一直被水利、岩土、石油勘探工程界关注。对深埋地下工程来说，岩体常常处于复杂的应力状态[1]，含裂纹岩体常常容易发生破坏。针对受压岩体裂纹的起裂与扩展，国内外许多学者进行了相应的探讨，Griffith[2]最早基于玻璃提出了近代断裂经典理论；L.Obert、Hoek.E[3- 4]对双轴条件下的裂纹扩展规律进行了试验研究；于庆磊等[5]对含中心穿透型裂纹的巴西圆盘进行了室内试验，研究了不同角度预制裂纹在轴压下的裂纹扩展演化规律。但是，以上研究仅仅是针对穿透型裂纹(二维情况)来说，事实上，裂纹往往以三维形态存在于岩体中，且三维裂纹扩展的形态与规律与二维裂纹相比有较大的不同，将三维裂纹简化为二维裂纹进行理论分析及数值模拟会丢失掉很多三维信息[6-7]。因此，对三维裂纹在复杂荷载下的裂纹起裂与扩展研究一直是学术界的研究热点。

水力压裂是近年来发展的新技术[8-9]，广泛利用于水利水电工程、石油和天然气工业、地热资源开发、核废料储存及地应力测量中[10-12]。基于以上应用，许多学者对水力致裂进行了一定的探讨，M. A. Kayupov等[13]对含孔洞的立方体试件进行了注水破坏试验；徐世烺[14]对比了楔入式紧凑拉伸试件在不同注水压下的裂纹扩展规律。但是，以上研究均未摆脱二维的范畴，对三维深埋裂纹在围压与水压共同作用下的裂纹扩展规律还缺乏一个客观的认识。

基于此，本文根据文献[1]提出的水力压裂裂纹扩展影响因子理论，运用三维裂纹分析软件Franc3D，同时考虑最大、最小主应力以及不同水压对三维深埋裂纹扩展的影响，对不同围压、不同内水压下的裂纹扩展过程及裂纹前缘应力强度因子的变化规律进行了数值模拟，为相应含裂纹岩体的裂纹扩展特性的认识提供参考。

1 计算理论

1.1 M积分

对于Ⅰ型、Ⅱ型、Ⅲ型的应力强度因子KⅠ，KⅡ，KⅢ，用M积分可以表示为[15]

(1)

(2)

(3)

(4)

式中，E为弹性模量；v为泊松比；KⅠ，KⅡ，KⅢ，分别为Ⅰ型、Ⅱ型、Ⅲ型应力强度因子；(1)、(2)为2种独立的线弹性状态；a、b、c分别为M积分的Ⅰ型、Ⅱ型、Ⅲ型应力强度因子分量。利用式(2)～(4)，便可计算出三维裂纹尖端的3个应力强度因子。

1.2 最大周向拉应变准则

Franc3D软件基于最大周向拉应变理论计算裂纹的起裂方向。对裂纹尖端来说，应力分量在极坐标下可以表示为

(5)

式中，σr为径向应力；σθ为环向应力，τrθ为切应力；θ为裂纹尖端角度。

对于线弹性材料，应力应变则满足胡克定律，即

(6)

式中，εr为径向应力；εθ为环向应力；γrθ为切应力。

最大周向拉应变准则认为，裂纹尖端沿着环向拉应变最大的地方扩展。同时，环向拉应变达到最大时裂纹开始扩展。

1.3 裂纹起裂方向

最大周向拉应变准则下，裂纹起裂方向满足以下关系

(7)

(8)

根据式(8)便可计算出最大周向拉应变准则下的裂纹扩展方向。

2 模型建立

2.1 计算模型

试件尺寸为标准立方体试件尺寸，长×宽×高为50 mm×50 mm×100 mm，其中裂纹尺寸为2a=20 mm、2b=15 mm(a、b分别为裂纹的长、短半轴长度)。材料的弹性模量E为16 GPa，泊松比μ为0.21，密度ρ为2.3 g/cm-3。计算边界条件为底部约束3个方向的位移，顶部施加应力边界为7 MPa，侧边施加应力边界为3 MPa，裂纹内表面施加水压边界。计算模型及网格见图1。

图1 计算模型及网格

2.2 水力压裂裂纹扩展影响因子的定义

为综合考虑最大与最小主应力以及裂纹面内水压的影响，根据康红普院士的研究成果，定义为

(9)

式中，D为水力压裂裂纹扩展影响因子；P为裂纹面的内水压力，σ1、σ3分别为最大、最小主应力。本文依次取D值为0、1、5，10和50进行分析计算。

2.3 模型验证

为说明本文计算的合理性，对比李术才等[16]运用树脂及中空裂隙注水进行双轴水力压裂试验，文献[16]的内水压力为10 MPa、σ1=10 MPa、σ2=2 MPa(即按照本文的定义，D=0)。裂纹扩展形态与本文数值模拟结果对比见图2。从图2可知，两者结果高度相似，即裂纹上下尖端出现了翼形包裹状裂纹，深色染剂部分即为翼形包裹裂纹的典型形状。证明了Franc3D软件在计算三维裂纹扩展方面有着独特的优势。

图2 结果对比

3 计算结果分析

3.1 不同工况下的裂纹扩展过程

水力压裂裂纹扩展影响因子D分别为0、1、5、10和50下的裂纹扩展过程见图3。从图3可以看出，不同D下的裂纹扩展规律有较大的差异：

(1)D较小时(对应D=0时)，裂纹面的水压力值较小，裂纹呈现典型的翼形包裹状裂纹扩展规律，即预制裂纹尖端的裂纹扩展与原裂纹平面呈现较大的转角，最终翼裂纹沿着最大主应力方向发展(本文中沿着应力边界为7 MPa，即试件的上下方向)，最终贯穿试件，试件发生破坏。

(2)D适中时(对应D=1、5、10)，预制裂纹先发生翼裂纹扩展，值得注意的是，此时翼裂纹与原裂纹面的夹角较小，且随着水力压裂裂纹扩展因子的变大，翼裂纹与原裂纹面的夹角越来越小，即逐渐趋向于沿着原裂纹面进行自相似扩展，翼裂纹扩展到一定程度时，裂纹面逐渐沿着平面扩展，且随着D的变大，扩展平面与原裂纹平面的夹角也越来越小。

(3)D较大时(对应D=50)，裂纹扩展与其他情况不同，预制裂纹尖端没有经历翼裂纹的扩展直接呈现自相似扩展，最终试件的断裂面也是沿着原裂纹面方向，值得注意的是，对裂纹自相似扩展，裂纹在未穿透试件表面一直呈现圆弧状扩展，直到裂纹穿透表面，裂纹前缘才会变直。

图3 不同D下的裂纹扩展过程

3.2 预制裂纹尖端特征点裂纹扩展规律

为定量监测不同工况下的裂纹特征点在每个计算步的扩展长度，根据对称性，选择如图4所示的裂纹特征点(点A和点B，分别对应于长轴端点与短轴端点)进行裂纹扩展长度的定量监测，图中数字为每一扩展步的关键点标识，C、D为A、B轴对称点。裂纹扩展的总长度与计算时步的关系见图5。

图4 裂纹特征点

图5 裂纹扩展长度随计算时步的变化

从图5可知，A点扩展长度随计算时步呈线性增大，在不同工况下，裂纹在前期(第1个计算时步)扩展长度无明显的差异，在裂纹的扩展后期(第2个计算时步后)，随着D的变大，扩展长度越大，但不同工况下A点的扩展长度之间差异较小。B点扩展长度随计算时步几乎呈线性增大，在第1个计算步之前，不同工况下的裂纹扩展长度无明显差异，但在裂纹扩展后期差异明显出现，即D越大，裂纹扩展总长度增长得越快。需要注意的是，B点总体上的裂纹扩展总长度要大于A点，这与数值模拟中观察到的“裂纹前后首先贯穿试件形成贯穿性裂纹”这一现象相吻合。

3.3 裂纹前缘应力强度因子分析

为研究不同工况下裂纹前缘的应力强度因子变化规律，A点为起点，定义初始裂纹前缘ABCD方向总长度为单位1(见图4)，并对3个应力强度因子KⅠ、KⅡ、KⅢ进行归一化处理，即KⅠ、KⅡ、KⅢ同时除以KR，KR可表示为如下形式[17]

(10)

Q=1+1.464(a/b)1.65

(11)

式中，σ0为模型边界上的拉伸应力。由式(10)、(11)即可得出归一化的应力强度因子随裂纹尖端路径的变化规律。

为体现Ⅰ型应力强度因子不同工况的差异，取y轴为对数坐标轴，而Ⅱ型、Ⅲ型应力强度因子的y轴则为正常坐标轴，绘制归一化应力强度因子变化规律，见图6。由图6可知：

(1)对于Ⅰ型应力强度因子来说，由于内水压力的作用，使得Ⅰ型应力强度因子都大于0，即裂纹面处于张拉的状态。Ⅰ型应力强度因子在裂纹的短半轴(图4中B点和D点)取得最大，而在裂纹的长半轴(图4中A点和C点)取得最小，这也解释了裂纹短半轴的扩展速率要大于长半轴，与数值模拟裂纹首先贯穿试件的前后表面相吻合。同时，D越大，Ⅰ型应力强度因子也越大，说明裂纹面水压力的增大会加剧裂纹的张开效应，由数值模拟中D较大时裂纹呈现自相似扩展规律可知，Ⅰ型应力强度因子是水力压裂情况下的裂纹扩展方向的主控因素。

(2)对于Ⅱ型应力强度因子来说，Ⅱ型应力强度因子的绝对值在裂纹的长轴顶点处取得最大，说明裂纹的长半轴顶点首先发生Ⅱ型裂纹破坏，结合数值模拟中翼形包裹状裂纹出现在裂纹上下尖端，说明Ⅱ型裂纹是翼形包裹状裂纹发生的主控因素。同时注意到，在预制裂纹尖端的短半轴，Ⅱ型应力强度因子达到最小为0，而Ⅰ型应力强度因子达到最大，说明裂纹的短半轴处发生纯Ⅰ型裂纹扩展。D越大，总体上Ⅱ型应力强度因子的绝对值越大，但不同工况下Ⅱ型应力强度因子的绝对值差异不大。