, 张砚琨, 赵岐, 梁力 东北大学 资源与土木工程学院，辽宁 沈阳 110819
收稿日期：2020-12-07
基金项目：中央高校基本科研业务费专项资金资助项目(N2001015)；国家留学基金委基金资助项目(201906085021)；国家自然科学基金资助项目(51474048)。
作者简介：李明(1980-)，男，辽宁沈阳人，东北大学副教授;
梁力(1955-)，男，辽宁丹东人，东北大学教授，博士生导师。
摘要：基于流固耦合理论，在基于渗透率的水力压裂(PHF) 模型基础上建立了扩展渗透率的水力压裂(extended permeability-based hydraulic fracture，EPHF) 模型.该模型假设岩石材料在开裂过程中渗透率的改变为平均有效应力的函数，采用多孔介质弹性本构关系和Drucker-Prager塑性模型描述岩石材料开裂前后的力学行为，孔隙流体采用达西定律描述.引入LET模型和混合黏滞系数考虑等效开裂区域内由压裂液与原液体混合引起的渗透率修正.分析了抗拉强度对关键位置的水压力和应力路径发展过程的影响，以及等效开裂区域和裂缝高度随注水时间的变化过程.结果表明：裂缝路径上的点(不含注水点) 的水压力先降低后增加至传播压力；当有效应力差和初始水压力同时保持不变时或当总应力不变时，起裂压力随有效应力比的增加而降低.
关键词：水力压裂多孔岩石材料起裂压力EPHF模型ABAQUS软件
Characteristics of Hydraulic Fracture in Homogeneous Porous Rock Material Based on EPHF Model
LI Ming
, ZHANG Yan-kun, ZHAO Qi, LIANG Li School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: LI Ming, E-mail: liming@mail.neu.edu.cn.
Abstract: Based on a coupled fluid-solid theory, an extended permeability-based hydraulic fracture(EPHF) model was developed according to the permeability-based hydraulic fracture(PHF) model. The proposed EPHF model assumed that the permeability of porous rock was a function of mean effective stress when a rock was fractured. The porous elastic constitutive law and Drucker-Prager plastic model were adopted to model the pre-failure and post-fracture process of rock material, respectively. The Darcy's law was used to describe the liquid phase. To consider the influence of injection fluid, the LET model was integrated to modify the permeability of fracture zone with the mixed viscosity of injection and original fluid being used. The development of pore pressure, stress path, equivalent fracture zone and fracture height for different tensile strength were studied by using the proposed model. The simulation results indicated that the pore pressure in fracture path(except injection point) first decreased and then increased to propagation pressure. If the effective stress difference and pore pressure were kept the same or if the total stress was constant, the breakdown pressure decreased when the effective stress ratio increased.
Key words: hydraulic fractureporous rock materialbreakdown pressureextended permeability-based hydraulic fracture(EPHF) modelABAQUS software
我国油气资源丰富，储存油气的岩石具有低渗透和难开采的特点，而水力压裂技术是开采此类油气的重要手段.水力压裂的分析方法主要包括理论求解、物理模型试验和数值模拟等.采用理论方法得到的平面及三维水力压裂模型有PKN，PKN-C，KGD，KGD-C，P3D和P3D-C模型等[1-2].理论模型虽然可以得到水力压裂相关问题的精确解，但是其求解过程均建立在严格的假设基础上，对于复杂因素影响下的求解过程仍较为困难.水力压裂的物理模型试验可以综合考虑多种影响因素，例如，Lu等[3]采用真三轴物理试验方法对比研究了水力压裂过程和气(氮) 压裂过程，并得到两种压裂条件下的起裂压力和裂缝传播等的特性.与理论求解和物理试验方法相比，丰富的水力压裂数值计算方法可以更有效地考虑多种复杂的影响因素.例如，Li等[4]基于流固耦合理论采用弥散裂缝模型探讨了大尺度具有随机分层特征的非均质岩石模型的水力裂缝动态发展过程；Gutierrez等[5]采用扩展有限元法研究了三层岩石材料在多点注水过程中储油层内压裂顺序引起的应力阴影效应；Shimizu等[6]基于离散元法研究了水力压裂过程受岩石材料的脆性(弹性模量和泊松比的函数) 以及初始应力的影响；Oliaei等[7]采用无网格伽辽金法研究了饱和多孔岩石材料的水力裂缝的起裂和传播过程；Zhuang等[8]基于比奥多孔弹性理论采用相场法实现水力压裂过程中裂缝与岩石材料分层界面相交过程的裂缝发展特征的模拟.
在数值计算方法中，弥散裂缝模型的求解过程无需引入真正的裂缝，易于实现且求解速度快.该法由Bazant等[9]提出用于求解混凝土裂缝发展过程.随后，Hu等[10]采用该法研究了岩体的水力压裂过程.Li等[11]基于弥散裂缝模型的思想，假设渗透率为平均有效应力的函数，在流固耦合理论基础上采用线弹性本构关系和摩尔-库仑塑性模型建立了基于渗透率的水力压裂(permeability-based hydraulic fracture，PHF) 模型，结合水平集法实现了非均质岩石材料建模与分段水力压裂过程的模拟[12].
本文在PHF模型基础上进一步开展研究，建立扩展基于渗透率的水力压裂(extended permeability-based hydraulic fracture，EPHF) 模型，实现饱和多孔岩石材料中水力裂缝发展过程的模拟，并探讨关键点位置的水压力和应力路径、等效开裂区域和裂缝高度随注水时间变化特征，并进一步开展起裂压力受初始应力状态、初始水压力和材料性质等因素的影响研究.
1 EPHF数值模型1.1 基本方程EPHF模型基于流固耦合理论，采用多孔介质弹性本构模型描述饱和岩石材料中固体骨架开裂前的力学行为，孔隙流体采用达西定律描述.其中平衡方程、渗流方程、变形方程和有效应力原理分别如式(1)~式(4) 所示.
 | (1) |
 | (2) |
 | (3) |
 | (4) |
岩石材料开裂前的力学行为由多孔弹性本构模型描述，见式(5)和式(6)：
 | (5) |
 | (6) |
水力压裂过程中岩石材料开裂后的等效开裂区域采用Drucker-Prager塑性模型描述，该模型可以计算由张拉和剪切引起的破坏，其p-q空间中的屈服方程F为
 | (7) |
1.2 水力压裂过程中的渗透率修正EPHF模型假设当岩石出现开裂时，渗透率随着平均有效应力的增加而迅速增大，对应的渗透率修正如式(8)~式(10)所示：
 | (8) |
 | (9) |
 | (10) |
沿用PHF模型中通过引入油水相对渗透率krw和混合黏滞系数μmix考虑压裂液与地层流体混合引起的渗透率改变[12]，即对式(8)计算得到的渗透系数采用式(11)进一步修正：
 | (11) |
饱和多孔岩石材料的两相相对渗透率模型可通过相对渗透率曲线拟合得到[15]，用于描述饱和度和相对渗透率关系的模型有Corey模型[16]、Chierici模型[17]和LET模型[18]等.与PHF模型中采用的单一参数描述相对渗透率曲线的Corey模型不同的是，在本文的研究中引入了包含三个描述曲线特征参数的LET模型计算相对渗透率，即
 | (12) |
图 1(Fig. 1)
 | 图 1 饱和度与相对渗透率关系(LET模型)Fig.1 Relationship between water saturation and relative permeability (LET model) |
1.3 EPHF模型在ABAQUS中的实现软件ABAQUS的二次开发功能提供了材料参数为场变量函数的基础功能.本文利用该功能增加了饱和岩石材料的渗透系数与平均有效应力的关系，从而采用该软件提供的用户子程序实现多孔岩石材料中水力压裂过程计算的EPHF模型.在该模型的实现过程中使用的关键用户子程序有UEXTERNALDB，SDVINI和USDFLD，分别用于岩石材料参数的读取与模型参数初始化、状态变量初始化和渗透率修正.使用软件ABAQUS提供的材料参数可设置为场变量的函数，EPHF模型共使用了8个场变量FIELD(1∶8) 和存储对应初值的状态变量STATEV(1∶8)，对应的参数依次为：对数体积模量κ，泊松比ν，抗拉强度Rm，渗透系数K，密度ρ，饱和度Sw，孔隙比e和混合黏滞系数μmix.各子程序对应的算法如表 1所示.
表 1(Table 1)
 表 1 实现EPHF模型的子程序算法Table 1 The algorithm of EPHF model in user subroutine
| 表 1 实现EPHF模型的子程序算法 Table 1 The algorithm of EPHF model in user subroutine |
2 饱和多孔岩石材料水力裂缝特征2.1 多孔岩石材料有限元模型建立以长度和宽度均为100 mm的二维岩石模型为例，探讨采用EPHF模型计算单点注水过程中的水力裂缝特征并与KGD模型的理论解对比.模型受到垂直和水平方向的初始应力分别为σv=70 MPa和σh=60 MPa，以及初始水压力pw=50 MPa.在模型中心设置宽度为1 mm的注水平面(注水点)，并施加持续时长为600 s的注水荷载，即注水速率qinj=0.001 mm/s.力学模型简图如图 2a所示，采用CPE8P单元建立的平面应变有限元模型如图 2b所示，其中单元宽度为1 mm，与注水宽度一致.计算中所采用的材料参数如表 2所示.
图 2(Fig. 2)
 | 图 2 力学模型与有限元网格Fig.2 Mechanical model and FEM meshes (a)—力学模型简图；(b)—有限元网格. |
表 2(Table 2)
表 2 材料参数Table 2 Material parameters
| 表 2 材料参数 Table 2 Material parameters |
2.2 水力裂缝特征以抗拉强度为4 MPa的情况为例，采用EPHF模型计算的水力裂缝等效开裂区域(即：受拉区域σ′m > 0) 随注水时间的发展过程如图 3所示.图 3中位置A，B和C分别为注水点、水力裂缝传播路径内点和等效开裂区域外点.
图 3(Fig. 3)
 | 图 3 Rm=4 MPa时裂缝等效开裂区域发展过程Fig.3 The development of equivalent fracture zone with Rm=4 MPa |
由于EPHF模型是弥散裂缝模型的一种，故而水力压裂的计算结果中无真实的裂缝.因此为获得裂缝的等效几何特征，假设裂缝的等效高度hf为等效开裂区域在y方向的范围，注水点位置的裂缝等效开度wf为等效开裂区域x方向外边界处的水平方向位移差(u+-u-)，其中lw大于等效开裂区域在x方向的宽度，如图 3所示.
基于KGD模型的水力裂缝高度hf、裂缝开度wf以及注水压力pw随注水时间变化的理论解可分别采用式(13)~式(15) 进行计算[19]：
 | (13) |
 | (14) |
 | (15) |
当抗拉强度Rm=4 MPa时，采用EPHF模型计算的水力裂缝等效高度、等效开度和注水点水压力随注水时间发展过程与对应的KGD模型的理论解分别如图 4~图 6所示.计算结果表明由EPHF模型计算的水力裂缝参数(裂缝高度hf、裂缝开度wf和注水压力pw) 与KGD模型的理论解保持良好的一致性.
图 4(Fig. 4)
 | 图 4 水力裂缝等效高度发展过程Fig.4 The development of equivalent fracture height |
图 5(Fig. 5)
 | 图 5 水力裂缝等效宽度发展过程Fig.5 The development of equivalent fracture width |
图 6(Fig. 6)
 | 图 6 注水点水压力随注水时间发展过程Fig.6 The development of injection pore pressure versus injection time |
2.3 关键位置的应力路径特征当抗拉强度分别为4，6和8 MPa时，图 3所示的三类位置A，B和C的应力路径发展如图 7所示(图例中的数值4，6和8均代表抗拉强度，S代表对应的应力路径起点).
图 7(Fig. 7)
 | 图 7 位置A，B和C的应力路径发展特性Fig.7 Development characteristics of stress paths at position A, B and C |
当Rm为4 MPa时，位置A，B和C的初应力状态(应力路径起点) 均在屈服面上，因此注水点A的平均有效应力p和等效剪应力q随着注水时间的增加沿着屈服面降低至开裂.位置B的p和q在开裂区域到达前沿着屈服面增加，到达后沿屈服面降低至开裂.位置C不处于裂缝区域内，受裂缝的影响较小，故p和q变化较小.当Rm为6和8 MPa时，随着注水时间的增加，注水点位置A的剪应力q略增加，平均有效应力p持续降低至屈服面，随后沿着屈服面继续降低至开裂.在裂缝到达之前，位置B的平均等效应力p维持不变而剪应力q增加直至屈服面，随后p和q沿着屈服面增加；当裂缝到达后p和q同时沿屈服面降低直至开裂.位置C仍未处于开裂区域内，故而与Rm为4 MPa情况相同.当抗拉强度分别为4，6和8 MPa时，位置A达到受拉状态的时间分别为15.95，9.42和9.42 s；位置B达到受拉状态的时间分别为302.4，319.4和382.3 s；位置C在整个注水过程中未开裂.
2.4 水压力发展过程与起裂压力不同抗拉强度下三类位置A，B和C对应的水压力发展过程如图 8所示.注水点A的水压力随着注水时间的增加而快速增加至起裂压力后降低至传播压力，对应为裂缝的起裂和传播过程.位置B为裂缝开裂区域经过点，当裂缝未到达该点时，水压力维持在初始状态；当裂缝达到该点时，水压力先降低后增加至传播压力，即该位置开裂完成.在注水过程中由于开裂区域未经过位置C，故而其水压力一直保持为初始水压力水平，变化较小.
图 8(Fig. 8)
 | 图 8 位置A，B和C的水压力发展特性Fig.8 Development of pore pressure at position A, B and C |
依据图 8的水压力变化曲线可以得到起裂压力和传播压力.结果表明，两者均随着抗拉强度的增加而增加.当前计算条件下，抗拉强度4，6和8 MPa对应的起裂压力分别为97.1，103.6和110.6 MPa.当假定岩石对压裂液分别为不可渗透和可渗透时，基于多孔弹性理论的起裂压力理论解分别为[20]
 | (16) |
 | (17) |
由式(16) 和式(17) 计算的起裂压力理论解(Pb1和Pb2) 和EPHF模型计算的数值解(Pbs) 及其相对理论解的相对误差(ePb1和ePb2) 如图 9所示.可以看出数值解与不可渗透理论解的相对误差较小，二者结果更为接近.这是由于数值计算中假定岩石材料渗透系数较低，同时水力压裂的开裂过程较快，故算例与不可渗透的理论假设较为符合.
图 9(Fig. 9)
 | 图 9 起裂压力对比Fig.9 Comparison of breakdown pressure |
3 起裂压力影响因素分析起裂压力的影响因素很多，如初始地应力状态、岩石材料性质分布的非均匀性和含有包裹体分布时其与注水点相对位置等[11].本文采用EPHF模型探讨饱和均质多孔岩石材料起裂压力受初始应力状态、水压力、材料性质和注水速率等因素的影响.计算中采用图 2所示的力学模型和有限元网格.当Rm为6，8和10 MPa时，探讨水压力(pw)、垂直有效应力(σ′v)、水平有效应力(σ′h)、有效应力差(Δσ′=σ′v-σ′h) 和有效应力比(ω′=σ′v/σ′h) 等的初始状态对起裂压力的影响.具体计算工况如表 3所示.
表 3(Table 3)
表 3 计算工况详情Table 3 Details of cases?
| 表 3 计算工况详情 Table 3 Details of cases? | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
当初始有效应力σ′v和σ′h分别取值为20 MPa和10 MPa且保持不变，水压力由35 MPa增加至55 MPa时(工况C1)，不同抗拉强度下的起裂压力变化规律如图 10所示.该计算结果表明，当有效应力保持不变时，不同抗拉强度水平下的起裂压力随着初始水压力增加而增加.
图 10(Fig. 10)
 | 图 10 水压力改变对起裂压力的影响Fig.10 Breakdown pressure influenced by the pore pressure |
当满足工况C2或工况C3时，有效应力比改变时对起裂压力的影响如图 11所示.
图 11(Fig. 11)
 | 图 11 有效应力比对起裂压力影响Fig.11 Breakdown pressure influenced by the effective stress ratio (a)—工况C2；(b)—工况C3. |
该计算结果表明，当有效应力差及水压力保持不变(工况C2) 和总应力保持不变(工况C3) 的情况下，不同抗拉强度水平下的起裂压力随着有效应力比ω′的增加而显著降低.当有效应力差Δσ′为5，10和20 MPa(工况C4~C8) 时，有效应力比改变和不同有效应力比下对应的水压力改变对起裂压力的影响如图 12所示.该计算结果表明，不同水压力水平下，起裂压力的改变与有效应力差变化特征一致(图 12a)；不同有效应力比的条件下，起裂压力随水压力增加而增加(图 12b).
图 12(Fig. 12)
 | 图 12 水压力和有效应力比对起裂压力影响Fig.12 Breakdown pressure influenced by the pore pressure and effective stress ratio (a)—pw=35~55 MPa；(b)—ω′=1.5~5.0. |
岩石的材料性质(对数体积模量κ、泊松比v、渗透系数K) 和注水速率qinj对起裂压力的影响如图 13所示.当渗透系数K在7.9×10-12mm/s至7.9×10-8mm/s之间，抗拉强度Rm分别取值为4，6，8和10 MPa时对应的起裂压力发展变化如图 13a所示.该计算结果表明不同抗拉强度水平下，随着渗透系数的增加对应起裂压力有不同程度的降低.当渗透系数较小时(如：K < 7.9×10-10mm/s)，渗透系数的增加对起裂压力影响较小；当渗透系数进一步增大时，起裂压力有不同程度的降低.
图 13(Fig. 13)
 | 图 13 材料性质与注水速率对起裂压力的影响Fig.13 Breakdown pressure influenced by the material properties and injection rate (a)—渗透系数的影响；(b)—κ，ν，K和qinj对起裂压力的影响(Rm=4 MPa). |
以抗拉强度为4 MPa时为例，当对数体积模量κ、泊松比ν和注水速率qinj分别在0.96×10-3~1.76×10-3，0.15~0.35和0.001~0.003mm/s变化时，首先归一化各变量，然后得到归一化后的参量Π改变对起裂压力的影响如图 13b所示.计算结果表明，注水速率qinj改变引起的起裂压力改变量为3.1%，说明其对起裂压力影响较小.而当对数体积模量κ、泊松比ν和渗透系数K增加时，在当前计算的参数范围内对应的起裂压力有不同程度的降低.
4 结论1) 通过与KGD模型计算的裂缝几何特征和水压力随注水时间变化的对比以及与基于多孔弹性理论计算的起裂压力对比，表明EPHF模型能够有效地计算出等效水力几何裂缝特征、水压力发展过程与起裂压力.
2) 采用EPHF模型计算了不同抗拉强度水平下水力裂缝发展过程，得到了注水点位置的水压力随着注水时间的增加而快速增加至起裂压力后降低至传播压力.开裂区域经过位置的水压力先降低后增加至传播压力，而开裂区域未经过的位置水压力变化较小.
3) 采用EPHF模型分析了起裂压力的影响因素.结果表明，在有效应力不变条件下，起裂压力随初始水压力的增加而增加；当总应力保持不变或有效应力差和水压力固定的条件下，起裂压力在不同抗拉强度水平下随有效应力比增加而降低；当初始水压力保持不变时，起裂压力受应力差影响.岩石的材料性质和注水速率等因素对起裂压力有不同程度的影响.
参考文献
| [1] | Rahman M M, Rahman M K. A review of hydraulic fracture models and development of an improved pseudo-3D model for stimulating tight oil/gas sand[J]. Energy Sources, Part A: Recovery, Utilization and Environmental Effects, 2010, 32(1): 1416-1436. |
| [2] | Shel E V, Paderin G V. Analytical solution of the pseudo-3D model for hydraulic fracturing in a storage-dominated regime[J]. International Journal of Rock Mechanics and Mining Sciences, 2019, 114: 92-100. DOI:10.1016/j.ijrmms.2018.12.020 |
| [3] | Lu Y, He M, Wu B, et al. Experimental and theoretical analysis of hydraulic fracturing and gas fracturing of rock under true triaxial compressions[J]. Engineering Fracture Mechanics, 2020, 234: 107100. DOI:10.1016/j.engfracmech.2020.107100 |
| [4] | Li M, Guo P, Stolle D, et al. Development of hydraulic fracture zone in heterogeneous material based on smeared crack method[J]. Journal of Natural Gas Science and Engineering, 2016, 35: 761-774. |
| [5] | Gutierrez R, Sanchez C M, Roehl D, et al. XFEM modeling of stress shadowing in multiple hydraulic fractures in multi-layered formations[J]. Journal of Natural Gas Science and Engineering, 2019, 70: 102950. DOI:10.1016/j.jngse.2019.102950 |
| [6] | Shimizu H, Ito T, Tamagawa T, et al. A study of the effect of brittleness on hydraulic fracture complexity using a flow-coupled discrete element method[J]. Journal of Petroleum Science & Engineering, 2018, 160: 372-382. |
| [7] | Oliaei M N, Pak A, Soga K. A coupled hydro-mechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG mesh-less method[J]. Computers and Geotechnics, 2014, 55: 254-266. DOI:10.1016/j.compgeo.2013.09.001 |
| [8] | Zhuang X Y, Zhou S W, Sheng M, et al. On the hydraulic fracturing in naturally-layered porous media using the phase field method[J]. Engineering Geology, 2020, 266: 105306. DOI:10.1016/j.enggeo.2019.105306 |
| [9] | Bazant Z P, Lin F B. Nonlocal smeared cracking model for concrete fracture[J]. Journal of Structural Engineering, 1988, 114: 2493-2510. DOI:10.1061/(ASCE)0733-9445(1988)114:11(2493) |
| [10] | Hu Y, Chen G L, Cheng W P, et al. Simulation of hydraulic fracturing in rock mass using a smeared crack model[J]. Computers and Structures, 2014, 137: 72-77. DOI:10.1016/j.compstruc.2013.04.011 |
| [11] | Li M, Guo P, Stolle D, et al. Modeling hydraulic fracture in heterogeneous rock materials using permeability-based hydraulic fracture model[J]. Underground Space, 2020, 5(2): 167-183. DOI:10.1016/j.undsp.2018.12.005 |
| [12] | Li M, Guo P, Stolle, D, et al. Heterogeneous rock modeling method and characteristics of multistage hydraulic fracturing based on the PHF-LSM method[J]. Journal of Natural Gas Science and Engineering, 2020, 83: 103518. |
| [13] | Brace W F, Paulding B W, Scholz C. Dilatancy in the fracture of crystalline rocks[J]. Journal of Geophysical Research, 1966, 71(16): 3939-3953. |
| [14] | Nicksiar M, Martin C D. Crack initiation stress in low porosity crystalline and sedimentary rocks[J]. Engineering Geology, 2013, 154: 64-76. |
| [15] | Moghadasi L, Guadagnini G, Inzoli F, et al. Interpretation of two-phase relative permeability curves through multiple formulations and model quality criteria[J]. Journal of Petroleum Science & Engineering, 2015, 135: 738-749. |
| [16] | Corey A T. The interrelation between gas and oil relative permeabilities[J]. Producers Monthly, 1954, 19: 38-41. |
| [17] | Chierici G L. Novel relations for drainage and imbibition relative permeabilities[J]. Society of Petroleum Engineers, 1984, 24(3): 275-276. |
| [18] | Lomeland F, Ebeltoft E, Thomas W H. A new versatile relative permeability correlation[C]//Proceedings of International Symposium of the Society of Core Analysts. Toronto, 2005: 1-12. |
| [19] | Wasantha P L P, Konietzky H, Xu C. Effect of in-situ stress contrast on fracture containment during single-and multi-stage hydraulic fracturing[J]. Engineering Fracture Mechanics, 2019, 205: 175-189. |
| [20] | Yew C H. Mechanics of hydraulic fracturing[M]. Oxford: Gulf Professional Publishing, 1997. |
| [21] | Krief M, Garat J, Stellingwerff J, et al. A petrophysical interpretation using the velocities of p and s waves(full-waveform sonic)[J]. The Log Analyst, 1990, 31(6): 355-369. |
| [22] | 马中高. Biot系数和岩石弹性模量的实验研究[J]. 石油与天然气地质, 2008, 29(1): 135-140. (Ma Zhong-gao. Experimental investigation into Biot's coefficient and rock elastic moduli[J]. Oil & Gas Geology, 2008, 29(1): 135-140.) |
