非对称应力作用下干热岩水力压裂裂缝扩展机理研究

曹 涵， 朱珍德， 周露明

（1.河海大学岩土力学与堤坝工程教育部重点实验室，南京 210098；2.河海大学江苏省岩土工程技术工程研究中心，南京 210098）

近年来能源安全问题已经逐渐引起人们的重视［1］. 地热能是一种可再生性的天然清洁能源，因为其具有低污染、低排放的特性，所以对于地热资源的开发利用引起了人们的高度关注. 现阶段，人类能够利用的地热能大多储存于高温高压、低孔低渗的干热岩中. 人们利用增强型地热系统（Enhanced Geothermal Systems，EGS）将这类干热岩通过人工手段改造为高渗透率的岩体，以此加热低温液体开采地热能. 由于干热岩高温高压、低孔低渗的特性，压裂液主要在水力压裂形成的裂缝网络（简称缝网）中流通，而建立增强型地热系统就是为了人为形成裂缝网络，以便更高效地提取干热岩中的地热能［2-5］.

增强型地热系统水压致裂过程中的热破裂是其所面临的重要问题，研究干热岩开发过程中的热流固耦合损伤（Thermal-Hydraulic-Mechanical-Damage，THM-D）问题有重要的科研价值. 李士斌等［6］分别分析了地层在拉伸破坏和剪切破坏下的裂缝扩展规律，就不同地应力形态对裂缝扩展的影响进行了初探. 袁志刚［7］通过研究实现了流固耦合损伤模型的有限元求解，同时对煤岩体水压致裂及裂缝扩展的影响因素进行了分析. Homand-Etienne和Houpert［8］对高温处理后的岩体试样进行了电镜扫描分析，同时分析了岩体微观结构的损伤状况对岩体物理力学性能的影响. 综合来看，现有的大量研究采用的都是热流固耦合损伤模型，且均为单井筒模型，同时关于热应力在水力压裂裂缝扩展中的作用的研究相对较少.

为研究非对称应力作用下干热岩双井筒模型的水力压裂裂缝扩展机理，本研究依据线弹性力学、能量守恒定律、质量守恒定律和弹-脆性损伤力学等理论建立干热岩水压致裂的热流固耦合损伤模型，然后使用有限元方法对非对称地应力和热应力影响下的双井筒模型的水力压裂裂缝扩展机理进行数值模拟研究.

1 热流固耦合损伤模型

1.1 应力场控制方程

根据弹性损伤本构关系，考虑有效应力、位移和热应力时线弹性体的静力平衡方程［9-10］如公式（1）所示.

式中：G为岩体剪切模量，Pa；v为岩体的排水泊松比；ui和Fi(i=x,y,z)分别为位移和体积力在i方向的分量；αp,i为水压力作用项；K′αTT,i为热应力项；K′为介质的排水体积模量，K′=2G(1+v)/3(1-2v)，Pa；αT为岩体的热膨胀系数，℃-1.

1.2 温度场控制方程

假设岩体与孔隙流体之间处于热平衡状态，则考虑热传导作用、能量守恒和体积变化的温度场控制方程［11］如公式（2）所示.

式中：λM=λs(1-φ)+λ1φ，λs和λ1分别为岩体和孔隙流体的热传导系数，W/（m·K）；(ρC)M为含流体的多孔介质的热容，(ρC)M=ρsCs(1-φ)+ρ1C1φ，kJ·m3/K；C1为孔隙流体的比热系数，kJ/（kg·K）；Cs为岩体的比热系数，kJ/（kg·K）；T0为零应力状态下的参考温度，K；ρ1为孔隙流体的密度，kg/m3；k为岩体的渗透率，m2；t为时间，s.

1.3 渗流场控制方程

基于流体的质量守恒方程以及状态方程，假设岩体与孔隙流体之间处于热平衡状态，则考虑变形和温度的渗流场控制方程如公式（3）所示.

公式（1）～（3）共同构成了热流固全耦合的控制方程组.

1.4 损伤分析

1.4.1 破坏准则

当岩体的应力状态达到最大拉应力准则(F1≥0)时，发生拉伸损伤；当岩体的应力状态达到摩尔-库伦准则(F2≥0)时，发生剪切损伤，F1和F2的表达式分别如公式（4）和公式（5）所示［12-13］.

式中：F1和F2分别为两个表示岩体应力状态的函数；fc和ft分别为岩体的单轴抗压强度和单轴抗拉强度，Pa.

选用弹-脆性损伤模型对岩体进行损伤判断，当岩体的应力状态满足F2＜0 且F1≥0 时，岩体发生拉伸损伤，损伤变量D可用公式（6）表示［14］.

式中：λ为抗拉强度的残余系数，λ=ftr/ft0；ft0为岩体的单轴抗拉强度，Pa；ftr为残余强度，Pa；εt0为岩体发生拉伸损伤时的拉应变；εtu为极限拉应变，εtu=ηεt0；在复杂应力条件下，ε=ε3，即为第三主应变.

当岩体的应力状态满足F2≥0 且F1＜0 时，岩体发生剪切损伤，损伤变量D可用公式（7）表示.

式中：εc0为岩体发生剪切损伤时的压应变.

1.4.2 损伤对模型参数的影响

根据弹性损伤理论和应变等效原理，岩体修正后的弹性模量［15-16］可用公式（8）表示.

式中：E0和E分别为岩体发生损伤前和发生损伤后的弹性模量，Pa.

岩体的应力状态影响着岩体的孔隙度，其关系式［17］可用公式（9）表示.

式中：φ0为岩体的初始孔隙度，干热岩可取值为0.01；φ为岩体当前的孔隙度；φr为高压缩状态下孔隙度的极限值；αφ为孔隙度的应力敏感系数，αφ=5.0×10-8Pa-1；σˉv为平均有效应力，其计算公式［17］如公式（10）所示.

式中：σ1、σ2、σ3分别为最大主应力、中主应力和最小主应力，Pa；p为孔隙水压力，Pa；α为Biot系数.

当岩石发生损伤后，其孔隙度和渗透性的演化规律较为复杂，可根据公式（11）描述其关系.

式中：αD为损伤对渗透率的影响系数，αD=5.0；k0为初始渗透率，m2；k为受到应力作用后的渗透率，m2.

虽然目前研究学者还没有在损伤变量对岩体热物性参数的影响方面形成统一的认识，且相关的理论和试验研究也比较少，但可以肯定的是，损伤必然会引起岩体热传导系数的升高. 本研究假定损伤是以指数关系的形式影响岩体热传导系数的［11］，其关系式如公式（12）所示.

式中：αλ为损伤对岩体热传导系数的影响系数；λs(T,D)为损伤发生后岩体的热传导系数；λs(T)为损伤发生前岩体的热传导系数.

2 数值模拟程序

由于热流固三场和损伤之间的耦合作用非常复杂，且考虑到热流固三场耦合控制方程组的高度非线性以及岩体的多个物性参数都处于动态的变化过程中，因此一般会使用有限元方法对热流固三场耦合控制方程组进行求解，其中较常用的有限元数值模拟软件为COMSOL，其易于进行耦合方程的建立，且具有较高的计算精度，能高效快速地求解非线性热流固三场耦合控制方程组.

本研究首先采用COMSOL软件建立并求解热流固三场耦合控制方程组，然后根据损伤对岩体物性参数的影响对各参数进行修正，之后再对热流固三场耦合控制方程组进行求解，如此重复，直到不再发生损伤为止.

3 非对称应力下水力压裂的数值模拟研究

3.1 数值模型的建立

为了研究非对称应力作用下干热岩双井筒模型的水力压裂裂缝扩展机理，建立了干热岩双井筒数值模型，如图1 所示. 模型尺寸为1 m×1 m 的正方形，双井筒的直径均为0.1 m，且为非中心对称. 模型的左边界和下边界固定位移，右边界和上边界施加水平地应力，双井筒筒壁施加水压力. 考虑到岩体的非均质性，可假定岩体的弹性模量服从Weibull 分布［18-20］. 因为本研究采用的双井筒模型为非中心对称的，所以地应力的分布情况将极大地影响裂缝扩展形态.故本研究分别在两种工况条件下对干热岩双井筒模型的水力压裂过程进行模拟，以得出非对称应力作用下干热岩双井筒模型的水力压裂裂缝扩展机理. 两种工况分别是：最大地应力方向垂直于双井筒连线方向（记为工况1）；最大地应力方向平行于双井筒连线方向（记为工况2）. 在进行数值模拟计算时，均固定初始地应力大小，为确定不同情况下水力压裂的起裂压力，以每步0.5 MPa 的速度增大孔壁的注水压力，直至模型发生破坏. 表1为数值模拟中采用的相关参数.

图1 干热岩双井筒数值模型示意图Fig.1 Schematic diagram of numerical model of hot dry rock double wellbore

表1 数值模拟中采用的相关参数Tab.1 Relevant parameters used in numerical simulation

3.2 数值模拟的结果分析

温度差是引起热应力的最主要因素，为了模拟热应力对干热岩双井筒模型水力压裂裂缝扩展的影响，在工况1的条件下，分别模拟了基岩温度为100、150、200 ℃时干热岩双井筒模型的水力压裂过程，其最终的损伤区域分布如图2所示.D=0表示岩体未发生损伤，0＜D≤1表示岩体发生了拉伸损伤，-1≤D＜0表示岩体发生了剪切损伤.

研究表明，当基岩温度为100、150、200 ℃时，岩体的破裂压力分别为35、32、30 MPa，说明随着岩体和注入压裂液之间温差的上升，岩体的破裂压力逐步降低. 由图2可知，工况1下，基岩温度为100 ℃时，裂缝首先沿着最大地应力方向发展，随后低温压裂液由裂缝进入岩体内部，低温压裂液与高温的裂缝壁面接触，在两者的温差作用下产生了热应力，热应力的存在加剧了后续的裂缝扩展，最终形成裂缝网络；基岩温度为150 ℃和200 ℃时，岩体在裂缝壁面产生了更多的分支裂缝，这使得整个裂缝网络更为复杂. 从图2还可以看出，当岩体和注入压裂液之间的温差升高到一定程度时，两个井筒压裂产生的裂缝网络会在损伤发展后期发生交汇，并在缝网交汇处迅速产生新的裂缝，这说明温度更高的岩体会产生范围更大、更广的裂缝网络. 综上可知，当最大地应力方向垂直于双井筒连线方向时，随着基岩温度的增加，岩体和注入压裂液之间的温差会逐渐增加，从而使热应力也随之增加，而热应力越大，岩体越容易发生损伤.

图2 工况1基岩温度对干热岩水力压裂损伤区域分布的影响Fig.2 Influence of different bedrock temperatures on the distribution of hydraulic fracturing damage area of hot dry rock under working condition 1

为了模拟注水压力对干热岩双井筒模型水力压裂裂缝扩展的影响，在工况2的条件下，分别模拟了注水压力为3.5、20、43 MPa（对应的迭代步数分别为7、40、86步）时干热岩双井筒模型的水压致裂过程，其最终的损伤区域分布如图3所示. 由图3可知，工况2下，随着注水压力的增大，岩体首先沿着最大地应力方向发生破坏，随着裂缝的不断发展，双井筒内侧区域中产生的裂缝发生交汇，裂缝扩展后期，双井筒外侧的裂缝网络逐渐趋于稳定，但双井筒内侧缝网交汇区的裂缝继续发展，最终导致岩体在双井筒内侧形成的损伤区域明显大于其在双井筒外侧形成的损伤区域，由此可见，双井筒内侧裂缝网络的交汇对后续裂缝的发展有明显的促进作用.

图3 工况2注水压力对干热岩水力压裂损伤区域分布的影响Fig.3 Influence of different water injection pressures on the distribution of hydraulic fracturing damage area of hot dry rock under working condition 2

4 结论

本研究建立了描述岩体细观结构的热流固耦合损伤模型，并对非对称应力和热应力作用下干热岩双井筒模型的水力压裂裂缝扩展过程进行了数值模拟研究. 得到如下结论：

1）当最大地应力的方向垂直于双井筒连线方向且基岩温度较低时，裂缝首先沿最大地应力方向发展，且此时两个井筒周围产生的裂缝网络互不干扰；随着基岩温度的升高，岩体破裂所需的注入压力降低，岩体的损伤区域面积则随之增大，当基岩温度升高到一定程度时，两个井筒周围的缝网发生交汇，在交汇处裂缝会继续扩展，且裂缝在交汇处的扩展速度快于其在无交汇损伤区域的扩展速度.

2）当最大地应力的方向平行于双井筒连线方向且注水压力较低时，裂缝同样首先沿着最大地应力方向发展，随着注水压力的增大，裂缝会在双井筒内侧区域发生交汇，内侧交汇区的裂缝扩展范围明显大于双井筒外侧.

3）地应力场对双井筒模型的围岩破裂形态影响显著，且最大地应力的方向直接决定了裂缝发展的方向和大致范围.