李明, 史艺涛, 李鑫, 梁力
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-10-20
基金项目：中央高校基本科研业务费专项资金资助项目(N170104024)；国家自然科学基金资助项目(51474048)。
作者简介：李明(1980-)，男，辽宁沈阳人，东北大学副教授;
梁力(1955-)，男，辽宁丹东人，东北大学教授, 博士生导师。

摘要：采用水平集法(LSM)描述非均质岩石的静态材料边界特性, 探讨分层特性、含有随机分布包裹体特性和具有Weibull分布特性的三种非均质岩石材料有限元建模方法, 并与基于完全流-固耦合弹塑性理论的弥散裂缝模型相结合进行水力压裂过程的模拟.在水力压裂模型中通过Corey关系和相对渗透率计算压裂液与原有液体混合的特性.最后进行符合Weibull分布的非均质岩石材料的水力压裂传播特性的计算, 给出了等效开裂区域的发展过程、注水点处水压力与注水时间变化曲线和有效应力路径.模拟结果表明, 采用水平集法能够有效地描述岩石材料的非均质特性, 同时该法能够与弥散裂缝模型很好地结合, 实现三维非均质岩石的建模与水力压裂过程的模拟.
关键词：水力压裂水平集法弥散裂缝模型非均质岩石流-固耦合
3D Modeling for Heterogeneous Rock Based on LSM and Characteristics of Hydraulic Fracture Propagation
LI Ming, SHI Yi-tao, LI Xin, LIANG Li
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: LI Ming, E-mail: liming@mail.neu.edu.cn
Abstract: The level set method(LSM)was introduced to describe the static material interface in a rock material. The finite element modeling procedures for the three-dimensional rock sample with the layered feature, randomly distributed inclusions and statistical Weibull distribution heterogeneity were given via the LSM. To study the propagation of hydraulic fracture, a smeared crack model based on the fully coupled, elasto-plastic, pore pressure-solid phase interaction theory was coupled to the LSM. The Corey relation and relative permeability were adopted to calculate the viscosity and permeability of mixture in the smeared cracked hydraulic fracture model. Finally, a case study that simulating hydraulic fracture propagation was carried out in a heterogeneous rock sample following Weibull distribution. The development of equivalent crack zone, pore pressure versus injection time and the effective stress path in the injected element were discussed. The simulation results indicate that the LSM coupling with smeared crack model can be well used for modeling the heterogeneity of three-dimensional rock material and simulating the hydraulic fracturing.
Key words: hydraulic fracturinglevel set method(LSM)smeared crack modelheterogeneous rockfluid-solid coupling
大量的天然气和石油存在于裂隙储层中, 水力压裂技术是主要的提取手段之一.水力压裂是用水压将岩石层压裂, 从而释放出其中的天然气或石油.数值计算方法是一种低成本高效率的预测裂缝传播手段.在水力压裂计算的三维分析中, 通常会消耗大量的计算资源, 但可以更好地反映裂缝发展特性, 如, 考虑流固耦合效应、岩石材料的非线性效应以及裂缝扩展的动态效应^[1]; 考虑压裂液的流变性、支撑剂在压裂缝中的运移、压裂液与储层岩石的热交换以及压裂过程中的流固耦合作用^[2]等.
本文重点研究三维条件下非均质岩石的建模方法以及应用弥散裂缝模型(smeared crack model)求解三维条件下的水力压裂过程, 该模型是除了粘结模型(cohesive zone model)^[3]和扩展有限单元法(XFEM)^[4]以外的又一个有效的数值计算手段.弥散裂缝模型是通过修改破坏区域的材料特性(如:刚度、渗透系数等)实现材料损伤过程, 其特点是计算中无需引入真正的裂缝, 但同时也不能得到明确的裂缝形状, 该模型适用于岩石材料的拉伸破坏过程的模拟^[5].此法进行数值计算的一个缺点是计算结果会受到网格形状的影响^[6].避免该问题的方法有引入特征长度和裂缝带宽^[7]或控制单元AR值^[8]等方法, 本文采用后者并通过引入水平集方程描述静态材料边界达到不同的岩石模型均能使用相同的有限元网格^[9].
水平集法(level set method, LSM)是由Osher等^[10]提出的追踪运动界面的方法.此法已经应用到多种数值计算中, 如:扩展有限元法^[11]、SPH法^[12]和接触优化分析^[13]等.高一维方程的使用增大了对计算机硬件的需求^[14].但是此法与弥散裂缝模型结合后的优点是可以通过隐式方法描述不同材料界面的特性, 并保持有限元网格不变.
本文应用水平集法实现三维条件下的非均质材料建模方法, 考虑材料分层、含有随机包裹体分布和符合Weibull分布的非均质特性三种情况, 并与水力压裂计算的弥散裂缝模型结合, 分析单点注水条件下三维非均质岩石材料的水压裂缝传播过程.
1 静态材料边界描述1.1 传统有限元法建模在建立含有材料边界的有限元模型时, 传统的有限元模型中单元边与材料边是必须重合的.以含有随机大小椭球形包裹体分布的三维岩石有限元模型为例, 采用传统的有限元建模方法的有限元网格如图 1所示.可以看到包裹体与基岩的材料界面不同时其网格特征是不同的.当岩石材料含有不同的包裹体分布时网格结构不同.另一方面, 对于具有包裹体等不规则材料边界的岩石模型划分网格时, 控制有限元网格为规则六面体是不可能的.
图 1(Fig. 1)

图 1 含有包裹体分布的三维岩石传统有限元模型Fig.1 Classical FEM of 3D rock matrix with inclusions

1.2 水平集法描述材料边界水平集法实现运动界面的追踪^[10], 是通过定义高一维的零水平集函数φ(x, t)=0来表示界面的, 如式(1)所示.

(1)

其中:x为计算域内的点坐标向量; t表示时间.
考虑含有包裹体分布的三维岩石模型, 当采用水平集方法建模时, 需记录所有包裹体的形状特征, 通过改变分布函数来实现非均质岩石模型的不同特征, 建模过程中无须重新划分有限元网格.当进行静态材料边界建模时, 可以略去式(1)中与时间相关的参数t.具有n个包裹体的水平集函数可表示为

(2)

其中:L_i为第i个包裹体的水平集方程；Ω_i为第i个计算域；φ_i(x)为第i个包裹体所满足的水平集方程(如, 椭圆形包裹体、圆形包裹体等^[11]).当L_i=1时, 点x处于包裹体内和边界上; 当L_i=0时, 点x位于包裹体外部.
2 三维非均质岩石建模在进行水力压裂特性分析时, 非均质材料需要考虑内部的非均匀特性, 本文探讨三维条件下的三种典型岩石模型的建模方法, 即:岩石中分布有任意大小的包裹体、岩石材料层状分布和非均匀性符合Weibull分布的岩石材料.
2.1 基于水平集法建模流程本文采用有限元软件ABAQUS进行非均质岩石的有限元建模及水力压裂过程分析, 计算中主要使用了SDVINI和USDFLD两个用户子程序.前者用于读取生成计算域内的非均质材料特性分布, 后者依据前者的材料分布用弥散裂缝模型计算水力压裂过程.非均质岩石建模过程及水力压裂分析流程如图 2所示.
图 2(Fig. 2)

图 2 水平集法有限元建模及水力压裂计算流程Fig.2 Schematic of finite element modeling with LSM and hydraulic fracture simulation procedure

2.2 计算域选取本文拟探讨-3 000 m处, 某非均质岩层的有限元模型建立方法及水力压裂过程的模拟.计算域选取边长为1 000 mm的立方体岩石单元进行非均质特性建模方法的研究, 分别考虑具有分层分布、含有随机包裹体分布和符合Weibull分布特性的三种情况模型.计算简图和所用有限元网格分别如图 3和图 4所示.
图 3(Fig. 3)

图 3 岩石材料计算域力学模型Fig.3 Mechanical model of solution domain

图 4(Fig. 4)

图 4 有限元网格Fig.4 Finite element meshes

2.3 含有包裹体分布特性的岩石材料建模在含有随机分布包裹体的岩石材料分析中, 以包含多个任意大小椭球形包裹体的岩石材料为例, 建立其有限元模型.其中第i个椭球形包裹体与基岩材料边界所满足的水平集方程为

(3)

其中:x_ij0为第i个包裹体的中心点坐标; a_j为椭球各轴长度最大值; α_j为用于改变轴长的随机数.当考虑椭球型包裹体与坐标轴间随机变化角度θ和λ时, 需进行坐标变换后再应用式(3)进行判断.包裹体生成时通过控制各包裹体中心点之间的距离实现两种模式, 一种是可以重叠, 另外一种是不可重叠.
当随机分布包裹体数量为4, 9, 10, 27, 150和300时有限元模型如图 5所示.可以看出, 虽然包裹体分布在不断地随机变化, 但有限元网格均为同一个.精确的材料边界可以通过增加网格密度的方式实现.
图 5(Fig. 5)

图 5 相同有限元网格下含包裹体分布的三维岩石模型Fig.5 3D rock model with randomly distributed inclusions under the same finite element meshes (a)—N=4; (b)—N=9; (c)—N=10; (d)—N=27; (e)—N=150; (f)—N=300.

2.4 具有分层特性岩石材料建模对于分层材料, 可以认为是符合特殊条件的水平集方程, 当只考虑材料沿着垂直方向(y)变化时, y处材料参数φ_i的水平集方程可表示为

(4)

其中, 位置y处的材料参数依据现场资料计算得到.
当所考虑的岩石材料弹性模量具有图 6所示的分层分布特性时, 采用水平集法建立分层有限元模型的仍然可以使用图 4所示的有限元网格.此时, 具有分层特性的岩石模型如图 6所示.
图 6(Fig. 6)

图 6 三维分层特性岩石建模Fig.6 3D rock modeling with layered features

2.5 非均质特性符合Weibull分布的岩石材料建模当岩石的细观物理属性满足Weibull分布时, 任意位置的材料属性是由满足此概率分布的水平集方程即式(5)计算得到, 此时材料边界为一点.

(5)

其中：λ>0为分布的比例系数；变量u为0到1符合均匀分布的随机数；k为形状参数.
当岩石试件的抗拉强度平均值为5 MPa, 即比例系数λ=5和反应材料参数分布向平均值的集中程度的形状参数k=5时, 岩石材料的抗拉强度概率密度曲线和有限元模型分别如图 7和图 8所示.可以看出, 材料的随机分布特性在相同的有限元网格中得以实现.
图 7(Fig. 7)

图 7 抗拉强度概率密度分布Fig.7 Probability density function for tensile strength

图 8(Fig. 8)

图 8 三维非均质岩石模型抗拉强度分布Fig.8 Tensile strength distribution of 3D heterogeneous rock

3 水力压裂计算模型本文采用的弥散裂缝模型^[15]是基于完全流-固耦合理论, 并假设开裂后的岩石单元渗透率增大.由于渗透率κ与渗透系数K呈线性关系, 故而在软件ABAQUS中通过修改渗透系数即可.开裂后渗透系数的改变与平均有效应力σ′_m和岩石抗拉强度R_m的关系, 如式(6)和式(7)所示:

(6)

(7)

式中：Kⁱⁿⁱ和K^up分别为计算域内初始渗透系数和开裂后最大渗透系数, 本文取K^up=0.1 m/s；ξ为影响渗透系数曲线变化的损伤系数.
当考虑两相流(压裂液和原液体)混合情况时, 采用相对渗透系数的计算方法考虑水力裂缝内混合后流体的渗透系数K_mix,

(8)

其中：k_r为相对渗透率, 计算中采用描述相对渗透率与饱和度S_w的Corey关系得到^[16]；μ_mix为开裂区域混合流体的黏度.
4 三维水力裂缝传播特性以符合Weibull分布的非均质岩石材料为例进行水力裂缝过程计算.算例中取地下-3 000 m处, 长度为1 m的立方体为计算单元, 初始应力状态为垂直方向应力σ_v=70.0 MPa, 水平方向两个应力分别为σ_hmax= 65.0 MPa和σ_hmin= 56.8 MPa, 水压力为49.0 MPa.计算简图和所用有限元网格分别如图 3和图 4所示.拟注水点位置为模型中心位置, 注水速率为0.005 mm/s.
当Weibull分布中弹性模量(E)、泊松比(μ)、抗拉强度(R_m)和渗透系数(K)的平均值分别取值为30 GPa, 0.3, 5.0 MPa和6×10^-9 mm/s时, 形状参数k取值分别为4, 15, 5和0.8时, 材料属性的三维分布情况如图 9所示.
图 9(Fig. 9)

图 9 三维非均质材料特性分布Fig.9 Characteristics distribution of 3Dheterogeneous rock (a)—弹性模量；(b)—泊松比；(c)—抗拉强度；(d)—渗透系数.

4.1 水力压裂区域开展过程当注水时间达到5 400 s时, 不同观测方向的水力压裂等效开裂区域如图 10所示.水力裂缝从注水开始到注水结束, 水力裂缝从1 005 s到5 400 s的动态发展过程如图 11所示.可以看出水压裂缝的传播方向为垂直于小主应力方向所在平面内.
图 10(Fig. 10)

图 10 水力压裂等效开裂区域(t=5 400 s)Fig.10 Equivalent hydraulic fracture zone(t=5 400 s) (a)—开裂区域透视图；(b)—开裂区域x-y方向视图；(c)—开裂区域z-y方向视图；(d)—开裂区域x-z方向视图.

图 11(Fig. 11)

图 11 水力压裂等效开裂区域发展过程Fig.11 Progress of equivalent hydraulic fracture zone

4.2 注水点特性注水点处的水压力变化曲线如图 12所示, 随着注水时间的增加水压力达到峰值(产生水力裂缝, 此时水压力为125 MPa)后降低到稳定状态(在注水条件保持不变的情况下裂缝形成后保持稳定状态, 此时水压力约为102 MPa), 该模拟结果反映了岩石破坏过程中注水压力的变化特点.
图 12(Fig. 12)

图 12 注水点处水压力变化曲线Fig.12 Pore pressure changing of injection point

注水点的应力路径变化特点与二维分析结果基本一致, 即:平均有效应力p由压应力转为拉应力达到岩石抗拉强度后出现水力压裂破坏，而有效剪应力q逐渐降低为零，如图 13所示.
图 13(Fig. 13)

图 13 注水点所在单元有效应力路径Fig.13 Stress path at injection element

5 结论1) 水平集法可以在保持有限元网格不发生改变的情况下, 有效地描述三维条件下岩石材料的非均质特性.
2) 详细探讨了采用相同的有限元网格实现具有随机包裹体分布特性、分层分布特性和非均质特性符合Weibull分布的三维非均质岩石材料的建模过程.
3) 通过对非均质特性符合Weibull分布的三维岩石模型的水力压裂过程分析, 可知水平集法与弥散裂缝能够有效地结合, 并得到符合工程实际的水压裂缝动态发展过程.
参考文献

[1]	朱君, 叶鹏, 王素玲, 等. 低渗透储层水力压裂三维裂缝动态扩展数值模拟[J].石油学报, 2010, 31(1): 119–123. ( Zhu Jun, Ye Peng, Wang Su-ling, et al. 3D numerical simulation of fracture dynamic propagation in hydraulic fracturing of low-permeability reservoir[J].Acta Petrolei Sinica, 2010, 31(1): 119–123.)
[2]	刘建军, 冯夏庭, 裴桂红. 水力压裂三维数学模型研究[J].岩石力学与工程学报, 2003, 22(12): 2042–2046. ( Liu Jian-jun, Feng Xia-ting, Pei Gui-hong. Study on mathematical model of three dimensional hydraulic fracturing[J].Chinese Journal of Rock Mechanics and Engineering, 2003, 22(12): 2042–2046.DOI:10.3321/j.issn:1000-6915.2003.12.017)
[3]	Priester C, Morton L C, Kinsey S T, et al. Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[J].Engineering Fracture Mechanics, 2012, 79(1): 312–328.
[4]	Ren Q, Dong Y, Yu T. Numerical modeling of concrete hydraulic fracturing with extended finite element method[J].Science in China Series E:Technological Sciences, 2009, 52(3): 559–565.DOI:10.1007/s11431-009-0058-8
[5]	Pak A.Numerical modeling of hydraulic fracturing[D].Edmonton: University of Alberta, 1997.
[6]	Bazant Z P, Lin F B. Nonlocal smeared cracking model for concrete fracture[J].Journal of Structure Engineer, 1988, 114(11): 2493–2510.DOI:10.1061/(ASCE)0733-9445(1988)114:11(2493)
[7]	Hu Y, Chen G, Cheng W, et al. Simulation of hydraulic fracturing in rock mass using a smeared crack model[J].Computers and Structures, 2014, 137(6): 72–77.
[8]	李明, 梁力, GuoPei-jun, 等. 弥散裂缝模型水力压裂数值方法的网格敏感性分析[J].东北大学学报(自然科学版), 2015, 36(9): 1337–1341. ( Li Ming, Liang Li, Guo Pei-jun, et al. Mesh sensitivity analysis of solution to hydraulic fracture problems based on a smeared crack model[J].Journal of Northeastern University(Nature Science), 2015, 36(9): 1337–1341.DOI:10.3969/j.issn.1005-3026.2015.09.026)
[9]	Li M, Guo P J, Stolle D, et al. Modeling method for a rock matrix with inclusions distributed and hydraulic fracturing characteristics[J].Journal of Petroleum Science & Engineering, 2017, 157: 409–421.
[10]	Osher S, Sethian J A. Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics, 1988, 79(1): 12–49.
[11]	Sukumar N, Chopp D L, Moes N, et al. Modeling holes and inclusions by level sets in the extended finite-element method[J].Computer Methods in Applied Mechanics and Engineering, 2001, 190(46): 6183–6200.
[12]	Marrone S, Colagrossi A, Touzé D L, et al. Fast free-surface detection and level-set function definition in SPH solvers[J].Journal of Computational Physics, 2010, 229(10): 3652–3663.DOI:10.1016/j.jcp.2010.01.019
[13]	Myslinski A. Level set method for optimization of contact problems[J].Engineering Analysis with Boundary Elements, 2008, 32(11): 986–994.DOI:10.1016/j.enganabound.2007.12.008
[14]	Stolarska M, Chopp D L. Modeling thermal fatigue cracking in integrated circuits by level sets and the extended finite element method[J].International Journal of Engineering Science, 2003, 41(20): 2381–2410.DOI:10.1016/S0020-7225(03)00217-9
[15]	Li M, Guo P J, Stolle D, et al. Development of hydraulic fracture zone in heterogeneous material based on smeared crack method[J].Journal of Natural Gas Science & Engineering, 2016, 35: 761–774.
[16]	Corey A T. The interrelation between gas and oil relative permeabilities[J].Producers Monthly, 1954, 19: 38–41.