深层干热岩储层具有温度高、硬度大和基质渗透率低等特点，需采用水力压裂等方法改造致密干热岩储层，构建复杂的裂缝网络体系，提高循环工质换热体积，实现高效循环取热^[1–6]。干热岩水力压裂是一个多物理场耦合问题，需要考虑位移、流体压力、温度等多个物理量之间的相互影响^[7–17]。H. Hofmann等人^[18]通过研究发现，复杂裂缝网络有助于提高干热岩储层开采时的经济效益。Zhang Yanjun等人^[19]提出，水平井分段压裂技术可提高循环工质的换热效率，从而提高裂缝型干热岩储层的开采效益。谢紫霄等人^[20]通过干热岩压裂物理模拟试验，发现水力裂缝沟通天然裂缝后能显著提高裂缝渗流系数，实现地热能的高效开发。干热岩储层裂缝扩展过程中受多因素综合影响，且温度影响不容忽略。

相场法是一种基于系统能量最小化原理判断裂缝扩展的有限元方法，在处理多物理场耦合问题时具有网格处理简单、结果准确性高和计算量小等优点，近年来被广泛应用于水力压裂数值模拟领域^[21–26]。侯冰等人^[23]基于相场法建立了流−固−化耦合酸压裂缝扩展模型，分析了孔洞结构对酸压裂缝的影响。N. Noii等人^[25]提出适用于正交各向异性的相场裂缝描述方法，对相场裂缝建模进行了初步分析。

相场法在处理复杂地层多物理场耦合裂缝扩展模拟方面具有显著优势。基于相场法建立了流−固−热耦合裂缝扩展数值模型，采用自适应网格方法平衡计算成本与计算精度，通过解析解验证模拟结果的准确性；在此基础上，研究了干热岩储层水力裂缝的起裂与扩展规律，分析了地质与工程因素对干热岩储层水力裂缝扩展的影响，以期为深层干热岩储层水力压裂提供参考。

1.   深层干热岩储层水力压裂数值模型

基于相场法构建流−固−热耦合裂缝扩展数值模型时，其核心思想是引入无量纲相场参数ϕ描述空间点的物理状态。ϕ的取值在 [0,1]区间，ϕ= 0 表示材料断裂，ϕ= 1 表示材料未断裂， 0 < ϕ < 1 表示材料处于裂缝过渡区域。该参数可以描述裂缝形态，同时可作为指示参数耦合多物理场方程。

1.1   控制方程

1.1.1   相场−位移控制方程

在热弹性多孔介质的相场裂缝模型能量泛函基础上^[24–25]，将裂缝扩展问题转化为系统能量最小化问题并求解，可得相场−位移控制方程强形式：

式中：u为位移，m；k为正则化参数；ε(u)为应变；ε为正则化参数；α_B为Biot系数；α_T为热膨胀系数，℃^－1 ；p为流体压力，Pa；p₀为初始压力，Pa；$\theta$为温度，℃；$\theta_0$为初始温度，℃；C为界面热应力常数，Pa/℃；σ⁺(u)为裂缝驱动力，Pa；σ^－(u)为裂缝非驱动力，Pa；在二维问题中，K_d=μ+λ，μ，λ为Lame系数；G_c为临界能量释放率，Pa·m。

1.1.2   温度场控制方程

基于相场参数的定义，可构建指示函数：

基于局部热平衡假设，同一空间点处温度值具有唯一性。引入指示函数后，得到温度场控制方程的有限元强形式：

式中：C_f为压裂液热容，kJ·m³/℃ ；C_r为储层热容，kJ·m³/℃ ；K_f为压裂液热传导系数，W/(m·℃)；K_r为地层岩石热传导系数，W/(m·℃)； v为流体流速，m/s；Q_i为注入热源，J。

1.1.3   流体压力控制方程

对于流体场，其物质守恒方程为：

式中：ρ_f为液体密度，g/cm³；φ^*为孔隙度；q为单位时间内注入的压裂液质量，g。

计算流体流速时忽略重力影响，表达式为：

式中：K为渗透率，D；η为液体黏度，mPa·s。

式中：$\varphi_0^*$为初始时刻孔隙度；M为Biot系数；K_T为孔隙热膨胀系数，℃⁻¹。

流体压力控制方程的强形式为：

1.2   自适应网格

模型中以相场参数为指示参数，引入自适应网格技术，在保证模拟精度的前提下，减少模型的计算量。自适应网格算法流程如图1所示。

图  1  自适应网格处理流程

Figure  1.  Adaptive mesh processing procedure

下载: 全尺寸图片 幻灯片

自适应网格可对同一时间步内的裂缝扩展状态进行多次求解，但与全局细分网格相比节点数量少、计算简单，计算精度仍保持良好^[26]。

1.3   多物理场耦合求解

模型共包括固体场模块、温度场模块和流体场模块等3个模块，模型中包括4个解变量，分别为相场参数（ϕ）、位移（u）、温度（T）和流体压力（p）。在固体场模块中，采用全耦合方法求解相场参数与位移；利用温度场模块、流体场模块分别求解温度和压力值，并采用固定应力迭代的方式与固体场耦合。具体流程如下：

1）在任意时间步tⁿ⁺¹，给定上一个时间步的解θⁿ，uⁿ，ϕⁿ和pⁿ。

2）基于固定应力法，构建迭代序列l=0，1，2…，定义初始解变量u_l=0=uⁿ，ϕ_l=0=ϕⁿ，p_l=0=pⁿ，θ_l=0=θⁿ。

3）给定θⁿ，uⁿ，ϕⁿ，pⁿ，u_l，ϕ_l和p_l，求解式（4），得到θ_l+1。

4）给定θⁿ，uⁿ，ϕⁿ，pⁿ，u_l，ϕ_l和θ_l，求解式（8）获得p_l+1。

5）给定θⁿ，uⁿ，ϕⁿ，pⁿ，p_l和θ_l，求解式（1）和式（2），可得{u_l+1，ϕ _l+1}。

6）若max{||p_l+1−p_l||，||θ_l+1−|θ_l||，||u_l+1−u_l||，||ϕ_l+1−ϕ_l||}≤TOL_fs，则迭代收敛（TOL_fs为固定应力法容差，取10 ^－8），令pⁿ⁺¹=p_l+1，θⁿ⁺¹=θ_l+1，uⁿ⁺¹=u_l+1，ϕⁿ⁺¹=ϕ_l+1；若迭代不收敛，则令u_l+1=u_l，ϕ_l+1=ϕ_l，p_l+1=p_l，θ_l+1=θ_l，重复步骤3）—5），直至满足收敛条件。

7）检查网格系统是否满足计算精度要求，若不满足，则更新自适应网格系统。

8）输出θⁿ⁺¹，uⁿ⁺¹， ϕⁿ⁺¹和pⁿ⁺¹作为时间步tⁿ⁺¹的解变量，并进行下一个时间步的计算。

1.4   初始条件与边界条件

模拟二维条件下的裂缝扩展过程，忽略重力影响与垂向应变。模型初始状态及边界条件如图2所示。

图  2  相场裂缝模型初始状态

Figure  2.  Initial condition of phase-field fracture model

下载: 全尺寸图片 幻灯片

模型长30 cm，宽30 cm。设置天然裂缝等结构时，将对应位置相场参数固定为0，其余岩石基质区域处相场参数设置为1，进行模拟分析，根据相场参数判断裂缝扩展情况。设定压裂液温度为20 ℃，岩石弹性模量为40 GPa，泊松比为0.20，岩石密度为2.50 g/cm³，Biot系数为0.80，渗透率为1 mD，地层岩石热容为2.0 MJ·m³/℃ ，压裂液热容4.2 MJ·m³/℃ ，地层岩石热传导系数为0.63 W/(m·℃)。

1.5   解析解验证

相场裂缝模型中，最小网格尺寸受网格细分次数控制。采用自适应网格方法，全局细分5次、局部细分3次时，最小网格尺寸为1.17 mm；采用全局网格细分方法，全局细分8次时，最小网格尺寸同样为1.17 mm。固定边界位移，分别采用全局细分网格与自适应网格对裂缝扩展进行模拟，模拟结果如图3所示。

图  3  固定边界位移条件下干热岩水力压裂数值模拟

Figure  3.  Numerical simulation of hydraulic fracturing in HDR reservoir under fixed boundary displacement

下载: 全尺寸图片 幻灯片

无限大线弹性物体内有一条缝内压力均匀分布的裂缝，其裂缝宽度的计算公式为^[27]：

式中：w为裂缝宽度，mm；p_f 为裂缝内流体压力，Pa；c为裂缝半缝长，mm；y为裂缝内某点到裂缝中心点的距离，mm；$\nu$为岩石泊松比；E为岩石弹性模量，GPa。

将流体压力代入式（9），计算裂缝宽度解析解，并与不同网格方案下裂缝宽度对比，结果如图4所示。

图  4  沿y轴方向缝宽分布情况

Figure  4.  Fracture width distribution along y-axis direction

下载: 全尺寸图片 幻灯片

模拟结果表明，裂缝宽度的数值解与解析解吻合程度较好，自适应网格与全局细分网格的模拟结果相接近。自适应网格方法节点少，算力消耗低，计算相场裂缝模型时有显著优势。

2.   地质与工程因素对裂缝扩展的影响

2.1   温度

温度场的变化会造成岩石体积变化，在系统内部互相约束，从而产生热应力^[28–29]。为研究温度对水力裂缝起裂与扩展的影响，设置初始地层温度分别为20，100和200 ℃，压裂液注入温度均为20 ℃，排量均为20 mL/min，边界应力分别为10和7 MPa，裂缝长度变化趋势如图5所示。

图  5  不同地层温度下干热岩储层水力压裂缝长变化趋势

Figure  5.  Variation trend of hydraulic fracture length in HDR reservoir with different formation temperature

下载: 全尺寸图片 幻灯片

模拟过程中，由于压裂液温度低于地层温度，热应力呈现为拉应力，且地层温度越高，热应力越大；由于水力裂缝以拉伸缝为主，热应力在裂缝扩展过程中起辅助作用。压裂液排量相同时，地层温度越高，裂缝扩展速度越快，起裂压力越低。真实干热岩储层中，岩石不同矿物热膨胀时存在各向异性，热应力方向与大小更为复杂，更容易产生局部应力集中；若地层中微裂缝发育，热应力易激活并诱导微裂缝起裂，提高水力裂缝复杂程度，且热应力越大，微裂缝激活程度越高。

2.2   压裂液排量

为探究不同压裂液排量对水力裂缝扩展的影响，设置压裂液排量分别为10，20和30 mL/min，地层初始温度均为200 ℃，压裂液温度均为20 ℃，边界应力分别为10和7 MPa，注液30 s后，模拟结果如图6所示。

图  6  不同压裂液排量下干热岩储层水力压裂模拟结果

Figure  6.  Simulation results of hydraulic fracturing in HDR reservoir under different pumping rates of fracturing fluids

下载: 全尺寸图片 幻灯片

从图6可以看出，压裂液排量为10 mL/min时，缝长为11.63 cm；压裂液排量为20 mL/min时，缝长为20.26 cm；压裂液排量为30 mL/min时，缝长为27.48 cm。研究表明，随着压裂液排量增加，液体憋压速度更快，裂缝周边岩石温度降低更明显。干热岩储层水力裂缝扩展受流体压力和热应力的共同驱动，裂缝扩展速度更快，裂缝宽度也随之增大，更有利于后续循环介质的注入与热交换。

2.3   天然裂缝夹角

为探究天然裂缝对水力裂缝扩展的影响，设置天然裂缝与水平最大主应力方向的夹角分别为90°，45°和20°，地层初始温度均为200 ℃，压裂液温度均为20 ℃，压裂液排量为20 mL/min，边界应力分别为10和7 MPa，模拟结果如图7所示。

图  7  不同天然裂缝夹角下干热岩储层水力压裂模拟结果

Figure  7.  Simulation results of hydraulic fracturing in HDR reservoir with different angles of natural fractures

下载: 全尺寸图片 幻灯片

水力压裂裂缝起裂后沿水平最大主应力方向扩展，遇到天然裂缝后发生转向并沟通天然裂缝。天然裂缝与水平最大主应力方向夹角为90°时，水力裂缝转向不明显，沟通天然裂缝后在天然裂缝两侧同时起裂并扩展，扩展角度与最大水平地应力夹角约为45°，2条裂缝呈分离趋势。岩石温度在水力压裂裂缝与天然裂缝沟通处降幅较大，主要原因为压裂液在该区域流动情况更复杂，传热效率高。

天然裂缝与水平最大主应力方向夹角为45°时，水力裂缝转向效果较明显，沟通天然裂缝后在天然裂缝两侧起裂并扩展，左右两侧扩展速度差异明显，左侧扩展速度更快，这是因为该侧与水平最大主应力方向夹角较小，扩展速度更快。另外，岩石温度在水力裂缝与天然裂缝右侧夹角处降低幅度较大，主要是该侧天然裂缝与水力压裂裂缝距离较近，传热效率更高。

天然裂缝与水平最大主应力方向夹角为20°时，天然裂缝诱导水力压裂裂缝产生转向现象。水力压裂裂缝沟通天然裂缝后在天然裂缝两侧起裂并扩展，左右两侧扩展速度差异明显，左侧扩展速度更快，裂缝沿水平最大主应力方向继续扩展。岩石温度在水力压裂裂缝与天然裂缝右侧夹角处降低幅度较大。

基于上述模拟结果可知，水力压裂裂缝遇到天然裂缝后易沟通天然裂缝，天然裂缝与水平最大主应力夹角越小，转向现象越明显。水力压裂裂缝与天然裂缝沟通位置流动情况复杂，传热效率高。水力压裂裂缝再次起裂方向受到地应力与天然裂缝方向的共同作用，在水平最大主应力方向夹角较小的方向扩展迅速，在水平最大主应力方向夹角较大的方向扩展较慢，且该侧岩石降温幅度更大。

2.4   多天然裂缝干扰下裂缝扩展

为探究多天然裂缝干扰下干热岩储层水力裂缝扩展力学行为，设置不同地应力与排量进行模拟。单侧天然裂缝条数为2条，天然裂缝与水平最大主应力方向夹角为45°，地层初始温度为200 ℃，压裂液温度为20 ℃，压裂液排量分别为20 和50 mL/min，边界应力差分别为3和8 MPa，模拟结果如图8所示。

图  8  不同天然裂缝密度下干热岩储层水力压裂模拟结果

Figure  8.  Simulation results of hydraulic fracturing in HDR reservoir with different densities of natural fractures

下载: 全尺寸图片 幻灯片

地层中存在多条天然裂缝时，水力压裂裂缝优先与距离最近的天然裂缝发生相互作用。地应力差为3 MPa、压裂液排量为20 mL/min时，水力压裂裂缝在靠近天然裂缝时轻微转向并开启初条天然裂缝，随后在天然裂缝尖端二次起裂；地应力差为8 MPa、压裂液排量为20 mL/min时，水力压裂裂缝靠近天然裂缝无明显转向现象，开启初条天然裂缝后在裂缝尖端二次起裂；地应力差为8 MPa、压裂液排量增加至50 mL/min时，水力压裂裂缝贯穿初条天然裂缝并继续扩展，由于水力压裂裂缝贯穿天然裂缝后会消耗水力能量，因此仅能贯穿单条天然裂缝。

3.   结　论

1）基于相场法建立了流−固−热多场耦合裂缝扩展数值模型，根据能量最小化原理计算裂缝形态，无需额外添加裂缝扩展准则，模型搭建简单，可用于干热岩储层水力裂缝扩展模拟研究。

2）干热岩储层水力压裂过程中，热应力体现为拉应力，有助于增加裂缝开度，在裂缝扩展过程中起辅助作用。压裂液与地层温度差越高，排量越大，热应力影响效果越显著；若地层中微裂缝发育，热应力易激活并诱导微裂缝起裂，提高水力压裂裂缝的复杂程度。

3）地层中存在多条天然裂缝时，水力压裂裂缝优先与距离最近的天然裂缝发生相互作用。水力压裂裂缝遇到天然裂缝后易沟通天然裂缝，天然裂缝与水平最大主应力方向夹角越小，地应力差越小，转向现象越明显。水力压裂裂缝再次起裂时沿天然裂缝方向起裂，且向水平最大主应力方向偏转。