李明, 赵岐, 陈昭, 梁力
东北大学 资源与土木学院，辽宁 沈阳 110819
收稿日期：2021-09-30
基金项目：中央高校基本科研业务费专项资金资助项目(N2001015)；国家留学基金资助项目(201906085021)。
作者简介：李明(1980-)，男，辽宁沈阳人，东北大学副教授;
梁力(1955-)，男，辽宁丹东人，东北大学教授，博士生导师。

摘要：天然非均质岩石具有复杂的细观结构和分布特征，本文结合柏林噪声建立考虑非规则边界的颗粒几何模型，实现考虑任意粒径、长短半轴、倾角和平面填充率的非规则颗粒填充.基于完全流-固耦合弹塑性理论，在已有的PHF-LSM(permeability-based hydraulic fracture level set method)模型基础上，通过引入MT (Mori-Tanaka)均匀化方法考虑积分点有效范围内细观尺度的材料分布特征，进而建立了PHF-LSM-MT非均质岩石材料细观特征建模及其水力压裂过程模拟方法.结合误差分析验证了该方法的可行性，在此基础上对非均质岩石材料进行了水力压裂过程模拟，得到了等效开裂区域.模拟结果表明：注水点水压力呈先升高后降低的趋势，起裂压力随包裹体体积分数增加而降低.
关键词：水力压裂细观建模柏林噪声岩石材料界面PHF-LSM
Study on Mesostructure Boundary Modeling and Hydraulic Fracture Propagation for Rock Material Based on the PHF-LSM-MT Method
LI Ming, ZHAO Qi, CHEN Zhao, LIANG Li
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Qi, E-mail: zhaoqi_326326@163.com.

Abstract: Natural heterogeneous rock material has complex mesostructure and distribution characteristics. An approach to model inclusions with irregular boundary by using the Perlin noise, and considering features, e.g., specific radius, semi-axial ratio, orientation angle and plane filling rate, is introduced. In order to include the mesostructure characteristics of rock material at integration point, the permeability-based hydraulic fracture level set method (PHF-LSM) is improved by introducing the Mori-Tanaka (MT) homogenization method based on a coupled fluid-solid elastoplastic theory. The proposed PHF-LSM-MT method is verified by error analysis. Finally, the propagation features of equivalent fracture zone in heterogeneous rock material were simulated by using the proposed method. The results indicate that the water pressure at injection point first increases and rapidly decreases, and the breakdown pressure decreases while the inclusion fraction increases.
Key words: hydraulic fracturemesostructure modelingPerlin noiseinterface of rock materialPHF-LSM(permeability-based hydraulic fracture level set method)
水力压裂技术作为油气增产手段已被广泛应用于实际工程并取得显著效果.目前关于水力压裂扩展规律的研究手段可分为：物理试验、理论研究和数值模拟.其中物理试验是探究岩石水力裂缝扩展特征、分布形态以及影响起裂压力因素的有效方法.例如Zhang等^[1]采用三轴水力压裂试验研究天然裂缝的空间位置和倾角对水力裂缝扩展形态的影响.
在理论研究方面，已有的方法多基于断裂力学的相关研究成果，如PK模型、KGD模型^[2]和PKN模型^[3]等.物理试验虽然可以较真实地模拟原位裂缝扩展形态，但成本较高，理论求解虽然可以得到精确解，但均建立在严格的假设基础上.
数值模拟方法具有低成本高收益的特点，可以弥补物理试验和理论研究的不足.例如：内聚区模型^[4]、弥散裂缝模型^[5]、相场法^[6]等.此外不同方法的结合亦可解决水力压裂的相关问题，如Yan等^[7]采用有限元法耦合离散元法实现了水力压裂过程的模拟，并通过物理试验证明压裂裂隙与地应力的关系.
岩石的宏观力学特性会受到细观结构特征的影响^[8]，因此有必要考虑天然岩石细观边界特征对力学特性和水力压裂过程的影响.如Li等^[9]结合水平集法(level set method, LSM)建立了含有圆形、椭圆形以及空间非均质的岩石材料模型，并结合基于渗透率修正的水力压裂模型(permeability-based hydraulic fracture, PHF)分析了水力裂缝的传播过程.
已有的PHF-LSM法可以考虑颗粒分布对水力压裂过程的影响，但其中未考虑积分点处细观边界特征.本文结合柏林噪声和椭圆水平集方程，建立以颗粒半径为参数的非规则颗粒建模方法，基于重心坐标法^[10]实现非规则颗粒填充，并引入Mori-Tanaka(MT)均匀化方法^[11]考虑高斯积分点处的细观边界特征.最终，建立PHF-LSM-MT法的非均质岩石材料的水力压裂数值计算模型，并分析考虑细观特征的非规则岩石材料的水力裂缝特性.
1 非规则颗粒边界方程在本文的研究中假设非规则颗粒Ω由颗粒包含的椭圆形(Ω₀)和非规则圆环(Ω₁)组成，如图 1所示.
图 1(Fig. 1)

图 1 非规则颗粒等效分解Fig.1 Equivalent decomposition of irregular particles (a)—非规则颗粒形状；(b)—边界离散化求解.

极坐标系下的非规则颗粒边界的水平集方程，如式(1)所示.

(1)

式中：r₀(θ)=(ab)/和r₁(θ)=AR_PN(θ)分别为椭圆和圆环的极径；a和b分别为椭圆形的长半轴和短半轴的长度；A为圆环厚度的均值；R_PN(θ)为一维柏林噪音(Perlin noise，PN)^[12].当需要考虑颗粒不同角度时，颗粒位置应先旋转后再应用式(1)进行判断.
2 颗粒填充和碰撞检测假设岩石内部的包裹体之间不重叠，在二维区域内生成多个非规则颗粒时，颗粒之间的关系可分为：包含、相交、相切和相离，如图 2所示.本文通过包裹体边界上的点是否在其他包裹体中的方法来判别颗粒间的相对位置关系.
图 2(Fig. 2)

图 2 非规则颗粒位置关系Fig.2 Relative position of irregular particles (a)—包含；(b)—相交；(c)—相切；(d)—相离.

2.1 非规则颗粒接触判定将颗粒离散为三角形单元集合，边界由线性差值简化，如图 3a所示，其面积由式(2)计算.

(2)

图 3(Fig. 3)

图 3 非规则颗粒碰撞检测方法Fig.3 Collision detection method of irregular particles (a)—非规则颗粒离散化；(b)—重心法原理.

式中：S_j为第j个三角形的面积；N为三角形个数，本文建模过程中取N=360.
重心法原理如图 3b所示，直角坐标系下三角形O₁O₂O₃的顶点坐标向量分别为O₁，O₂，O₃，点P′在三角形内部的权重公式为^[13]

(3)

式中：u，v均为常数，当u=0，v=0时, P′代表O₁点；当u=0, v=1时, O′代表O₂点；当u=1, v=0时, P′代表O₃点.根据图中P′，P₁，P₂，P₃可知，如果判断点在三角形内部，则必须要满足①u≥0；②v≥0；③u+v≤1，取等号时，表示判断点在三角形边上.令V₁=(P′-O₁)，V₂=(O₃-O₁)，V₃=(O₂-O₁)，将式(3)的等式两边分别乘以V₂和V₃即可得到v和u.
2.2 非规则颗粒间的位置关系建立非规则颗粒I₁和I₂并分别包含以下信息：中心坐标(x₀, y₀)；长半轴(a)；短半轴(b)；颗粒倾角(β)；圆环厚度均值(A)；颗粒边界顶点坐标(x_j, y_j).x_j和y_j分别代表第j个顶点的x坐标和y坐标，其中1≤j≤360.顶点坐标的储存顺序按照逆时针方向，当长半轴与水平夹角为0°时记为第一个顶点坐标.
在调用颗粒信息时采用以下书写格式(调用第1个颗粒的中心坐标记为I₁(x₀, y₀))：
① 颗粒I₁包含颗粒I₂：I₂(x₀, y₀)，…，I₂(x₃₆₀, y₃₆₀)均在I₁内部，判别式恒成立.②颗粒I₂包含颗粒I₁：I₁(x₀, y₀)，…，I₁(x₃₆₀, y₃₆₀)均在I₂内部，判别式恒成立.③颗粒I₂与颗粒I₁相交：I₂(x_k, y_k)，…，I₂(x_l, y_l)均在I₁内部，共(l-k+1)个判别式成立；或I₁(x_m, y_m)，…，I₁(x_n, y_n)均在I₂内部，共(n-m+1)个判别式成立，1≤k < l≤360，1≤m < n≤360.④颗粒I₂与颗粒I₁相离：I₂(x₀, y₀)，…，I₂(x₃₆₀, y₃₆₀)均不在I₁内部且I₂(x₀, y₀)，…，I₂(x₃₆₀, y₃₆₀)均在I₁内部，判别式恒不成立.
当颗粒在进行碰撞判定时会出现图 4所示的4种特殊情况.通过细化颗粒边界和控制最小距离D_min避免图 4所示现象，其中D_min与长半轴(a)和厚度均值(A)之间的关系见式(4).

(4)

图 4(Fig. 4)

图 4 非规则颗粒间特殊位置Fig.4 Special positions between irregular particles (a)—类型Ⅰ；(b)—类型Ⅱ；(c)—类型Ⅲ；(d)—类型Ⅳ.

式中：D_gap为颗粒间缝隙；I_ia和I_iA分别为新颗粒I_i的长半轴和厚度均值.
2.3 非规则颗粒间距计算以图 5所示五边形为例，在确定I₂的中心点V₀不在颗粒I₁内部的基础上，计算向量OV₀与水平方向夹角θ_t，在颗粒I₁中一定可以找到2个连续的顶点坐标V_j+1(x_j+1, y_j+1)和V_j(x_j, y_j)，该顶点坐标和I₁中心点连线与水平方向的夹角分别为θ_j+1，θ_j，三者之间的关系满足θ_j+1≥θ_t≥θ_j. OV_t的长度采用式(5)和式(6)计算，其中V_t为OV₀与V_j+1V_j的交点.

(5)

(6)

(7)

图 5(Fig. 5)

图 5 非规则颗粒间距计算Fig.5 Calculation of distance between irregular particles

通过式(4)和式(7)可计算不同颗粒间的最小距离D_gap.采用Processing软件在给定区域内填充不同几何特征的非规则颗粒模型，如图 6所示，其中a=(2.0~30.0)mm，b=(0.0~1.0)×a，β=(0~2π).
图 6(Fig. 6)

图 6 非规则颗粒填充Fig.6 Irregular particle packing (a)—A=10，φ1=20 %；(b)—A=10，φ1=30 %；(c)—A=15，φ1=20 %；(d)—A=15，φ1=30 %.

3 多尺度非规则边界有限元模型已有的基于水平集法建立的有限元模型如图 7所示.由于积分点的有效区域内会被赋值为单一材料属性，颗粒边界会呈现“锯齿状”, 如图 7b所示.
图 7(Fig. 7)

图 7 隐式建模方法Fig.7 Implicit modeling method (a)—背景网格嵌入颗粒前；(b)—背景网格嵌入颗粒后.

3.1 非均质岩石材料边界水平集方程及参数均匀化以图 8所示含有随机分布非规则颗粒(D_gap>0)，边长为l的岩石模型为例.该模型受到初始水压力(p_w)、水平和垂直方向应力(σ_h，σ_v)作用，模型底边受到垂直方向约束，底边中心位置为固定约束.模型中间施加注水荷载.由式(1)建立的水平集方程为

(8)

(9)

图 8(Fig. 8)

图 8 非均质岩石力学模型Fig.8 Mechanical model of heterogeneous rock

式中：x为平面Ω内积分点的位置向量；当φ_i < 0时表示积分点在颗粒内部，反之，当φ_i>0时表示积分点在颗粒外部，当φ_i=0时表示积分点在颗粒边界上.
当积分点的有效区域w_sub×h_sub内含有材料边界时(如图 7中积分点i₁，i₂，i₅，i₉)，已有的简化方法是通过阶跃函数、线性函数和非线性函数对材料进行过度处理^[14]，该三种方法均不能考虑材料边界特征的影响.本文采用均匀化方法对高斯积分点有效区域进行均匀化处理，以消除边界“锯齿状”现象.均匀化后的等效材料参数Λ由式(10)表示，即

(10)

其中：I₀代表基质的材料参数；I₁，…，I_n代表颗粒的材料参数.PHF模型中共考虑了4种材料参数，其中弹性参数(E, v)采用MT^[15]法计算，即：

(11)

(12)

式中：I =(3K, 2G)为等效刚度张量；K和G分别为等效体积模量和等效剪切模量；φ_r和I_r=(3K_r, 2G_r)分别为第r相颗粒的体积分数和有效刚度张量；应变张量H_r可通过式(12)求得，其中T_r为

(13)

式中：I和I₀是四阶单位张量和基质的有效刚度张量；S_r代表第相颗粒的四阶Eshelby张量^[16].
文献[17]的研究表明等效材料弹性参数均在Voigt上限(VR⁺)和Reuss下限(VR^-)之间，即

(14)

式中，和.
等效抗拉强度采用Voigt上限I_Voigt计算方法^[18].本文只研究夹杂单相颗粒的材料，当r=0和r=1时分别对应基质和颗粒材料.饱和岩石材料中等效渗透系数κ采用Halpin-Tsai(HT)方法计算^[19]，即

(15)

式中：M_a=κ₁/κ₀；κ₀和κ₁分别为基质和颗粒的初始渗透系数；φ₁为颗粒的体积分数；ζ是与颗粒形状有关的参数，本文采用的非规则颗粒是以椭圆形为基础，因此计算中ζ取0.5^[19].
3.2 颗粒边界细观模型验证本文依托ABAQUS软件进行数值模拟，通过用户子程序实现图 9所示积分点有效区域参数均匀化.步骤如下：①通过子程序UEXTERNALDB读取非规则颗粒参数信息，根据式(11)建立隐式材料边界模型，见图 9a.②通过子程序SDVINI调用单元信息，根据积分点类型获得坐标信息和有效区域范围w_sub×h_sub，见图 9b.③将有效区域离散成w_min×h_min子网格并计算每个子网格中基质和颗粒占有效面积的体积分数φ₀和φ₁，见图 9c.本文取值w_min=w_sub/10，h_sub=h_min/10.④根据式(11)、式(14)和式(15)对材料参数均匀化，并重新赋值到积分点，见图 9d.⑤重复循环完成所有单元检测.
图 9(Fig. 9)

图 9 积分点参数均匀化Fig.9 Homogenize parameters at integration point

通过建立颗粒边界并对不同区域赋值材料属性得到模型(显式模型)，其结果作为对照组，验证LSM-MT方法的可行性.以模型垂直方向杨氏模量E_y为例，对比LSM法与LSM-MT法建立模型(隐式模型、改进模型)相对于显式模型计算结果的相对误差e_r=|(Ξ_im -Ξ_ex)/Ξ_ex |和绝对误差ε_r=Ξ_im-Ξ_ex，其中Ξ_im为隐式方法计算结果，Ξ_ex为显式方法计算结果，其相对误差和绝对误差记为e_ri和ε_ri，改进模型的相对误差和绝对误差记为e_rm和ε_rm.
三种方法建立的含有非规则颗粒分布的有限元模型如图 10所示，可以看出隐式模型中材料边界呈“锯齿状”，而改进模型中材料边界有明显改善.同时采用表 1所示包裹体的几何参数，建立基质和颗粒的杨氏模量分别取值20 GPa和80 GPa，泊松比分别取值0.25和0.2的模型.通过不同方法计算的等效模量如图 11所示，结果表明，等效弹性参数E_y均落在VR上下界限内，且在MT法的解上下波动.三种模型的E_y均随颗粒体积分数增加而增大.
图 10(Fig. 10)

图 10 基于不同建模方法有限元模型Fig.10 Finite element model by various modeling methods (a)—显式模型；(b)—隐式模型；(c)—改进模型.

表 1(Table 1)

表 1 非规则颗粒几何参数取值Table 1 Geometry parameters of irregular particles 参数名称 取值范围 参数名称 取值范围 l/mm 100 A/mm 2, 3, 4 ai/mm 2~47 βi/(°) 0~360 bi/mm (0.1~1.0)ai φ1 0.05~0.80 Dgap/mm 0.5

表 1 非规则颗粒几何参数取值 Table 1 Geometry parameters of irregular particles

图 11(Fig. 11)

图 11 体积分数改变对弹性模量的影响Fig.11 Influence of volume fraction on elastic modulus

两种隐式模型在不同体积分数下的相对误差和绝对误差如图 12所示.
图 12(Fig. 12)

图 12 不同方法的有效模量误差对比Fig.12 Error of effective modulus between various methods

结果表明两种模型计算得到的弹性模量均大于理论值，且LSM-MT法的计算误差均小于LSM法.因此本文建立的LSM-MT建模方法可有效考虑材料边界分布特征，并增加计算结果的精度.
4 非均质岩石的水力裂缝特性采用基于Biot固结理论建立的PHF水力压裂数值计算模型，结合LSM-MT法建立含有非规则颗粒的岩石模型实现其水力压裂过程的模拟.
4.1 PHF方法的控制方程在深部岩层中水力压裂过程是渗流与裂隙岩体相互作用的过程，其平衡方程、渗流方程、几何方程、物理方程和有效应力原理如下^[9]：

(16)

(17)

(18)

(19)

(20)

式中：σ′_ij，σ_ij，p_w，ρ，μ，ε_ij和ε_v分别为有效应力张量、总应力张量、流体压力、岩石密度、流体黏度、应变张量和体积应变；u_i和f_i分别为x_i方向的位移和体力；模型的弹性部分采用线弹性本构关系，塑性部分采用Mohr-Coulomb破坏和强化准则；变量α和1/Q分别代表Biot固结系数和储存系数；裂缝开裂区域内渗透系数和平均有效应力的关系函数和压裂液黏度对开裂区渗透系数的影响如式(21)所示.

(21)

式中：参数κ_ave=(κ_max+κ₀)/2，κ_max和κ₀分别为渗透系数上限值和初始值；κ_δ=(κ_max-κ₀)/2；R_m和σ_m分别为岩石材料的抗拉强度和平均有效应力；ξ为损伤系数；μ_mix为混合黏滞系数；k_rw为水相对渗透率，采用Corey^[20]关系计算，即

(22)

式中：k_rw0，S_w，S_wi，S_orw和N_w分别为水相对渗透率端点、饱和度、束缚水饱和度以及剩余油饱和度；本文计算中k_rw0=0.5，N_w=3.66，S_wi=S_orw=0.15^[20].
4.2 PHF-LSM-MT水力压裂数值计算模型PHF-LSM-MT水力压裂数值模型的建立采用ABAQUS软件二次开发功能，流程如图 13所示.子程序SDVINI包含8个状态变量：弹性模量(E)、泊松比(v)、基岩抗拉强度(R_m)、有效渗透系数(κ)、基岩密度(ρ)、饱和度(S)、孔隙比(e)和混合黏滞系数(μ_mix).依据式(11)实现非规则颗粒材料分布，并采用均匀化方法修改材料边界处积分点有效参数.子程序USDFLD用来接受SDVINI传入的状态变量，并通过式(21)修正等效开裂区域渗透系数.
图 13(Fig. 13)

图 13 PHF-LSM-MT模型算法实现Fig.13 Algorithm of PHF-LSM-MT model

4.3 均质岩石材料水力裂缝扩展以边长为100 mm的二维均质岩石模型为例，将PHF模型计算得到的水力裂缝扩展形态与KGD模型的理论解进行对比.力学模型如图 8所示，材料参数采用某地下3 088 m的实测数据，如表 2所示.
表 2(Table 2)

表 2 基岩材料力学参数Table 2 Material properties of rock matrix 参数名称 数值 弹性模量E/GPa 22.719 泊松比v 0.283 抗拉强度Rm/MPa 4.172 摩擦角Ψ/(°) 30 膨胀角ψ/(°) 10 密度ρ/(kg·m-3) 2 538.32 渗透系数κ/(mm·s-1) 1.0×10-5 孔隙率e 0.072

表 2 基岩材料力学参数 Table 2 Material properties of rock matrix

模型受到初始水平和垂直方向地应力为σ_h= 61.47 MPa和σ_v= 71.62 MPa，初始水压力p_w= 52 MPa，注水速度q_inj=0.001 mm/s，注水时间为600 s.由于PHF模型是通过修改破坏区渗透系数表示裂缝扩展过程的一种弥散裂缝模型，因此模拟过程中没有出现真实的裂缝，计算得到的等效开裂区域发展过程如图 14所示.采用PHF法计算的裂缝等效高度h_f为开裂区域的垂直方向范围，裂缝等效宽度w_f为注水点位置裂缝开裂区域的水平方向位移差(u_r-u_l).
图 14(Fig. 14)

图 14 裂缝等效开裂发展过程Fig.14 Development of equivalent fracture zone (a)—t=12 s；(b)—t=421 s；(c)—t=600 s.

基于KGD理论模型的水力裂缝高度h_f、裂缝开度w_f和注水压力p_w的理论解如式(23)~式(25)^[9]所示：

(23)

(24)

(25)

式中：G为剪切模量；当水力裂缝为双向扩展时，参数a，b，c分别为0.48，1.32和1.91^[9].
通过PHF模型计算的水力裂缝等效高度、等效开度和注水点水压力随注水时间变化过程与KGD理论解对比如图 15所示.
图 15(Fig. 15)

图 15 水力裂缝特征验证Fig.15 Verification of hydraulic fracture characteristics (a)—裂缝高度扩展；(b)—裂缝宽度扩展；(c)—注水点水压力发展.

该计算结果表明PHF模型计算的数值解与解析解基本吻合.由于KGD模型假设裂缝处于无限半空间中，而在注水过程中随着裂缝不断扩展并逐渐接近模型边界，导致两者有一定偏差.
4.4 非均质岩石水力裂缝发展特性建立图 8所示的非均质岩石的数值计算模型，边界条件与均质岩石情况相同.计算中假设包裹体不易开裂，所用材料参数如表 3所示.每种包裹体体积分数包括3个随机算例.
表 3(Table 3)

表 3 包裹体材料力学参数Table 3 Material properties of inclusion sample 参数名称 数值 弹性模量E/GPa 80 泊松比v 0.20 抗拉强度Rm/MPa 40 渗透系数κ/(mm·s-1) 1.0×10-18 密度ρ/(kg·m-3) 2 538.32 厚度均值A/mm 2 体积分数φ1/% 10, 20, 30, 40

表 3 包裹体材料力学参数 Table 3 Material properties of inclusion sample

当注水时间为600 s时，基于PHF-LSM-MT法计算的非均质岩石模型等效开裂区域如图 16所示.该结果表明：当传播路径存在包裹体时，水力裂缝在材料交汇处发生变化.当包裹体与裂缝发展方向的有效面积和倾角均较大，裂缝沿包裹体表面扩展，如图 16a和图 16c所示；当包裹体在裂缝发展方向的有效面积较大但倾角较小时，包裹体会阻止裂缝扩展，随注水时间增加，裂缝沿包裹体表面发展，最终绕过包裹体，如图 16c和图 16d所示；当包裹体在裂缝发展方向的有效面积较小时，裂缝直接绕过包裹体继续发展，如图 16b和图 16d所示.
图 16(Fig. 16)

图 16 非均质岩石水力裂缝模型等效开裂区域Fig.16 Equivalent fracture zone of heterogeneous rocks (a)—φ1=10 %；(b)—φ1=20 %；(c)—φ1=30 %；(d)—φ1=40 %.

注水点位置的水压力变化曲线如图 17所示，两种建模方法得到的结果均呈现先增加至起裂压力后降低至传播压力的趋势，与已有的计算结果一致^[14].
图 17(Fig. 17)

图 17 注水点单元水压力发展历程Fig.17 Development of pore pressure at injection element (a)—φ1=10 %；(b)—φ1=20 %；(c)—φ1=30 %；(d)—φ1=40 %.

注水位置应力路径变化如图 18所示.结果表明包裹体分布影响了应力路径初始应力状态. 当注水时间增加时，平均有效应力p由初始状态降低至零(受压状态)，随后开始增加(受拉状态).剪应力q在开裂过程处于降低状态，最终为零.
图 18(Fig. 18)

图 18 注水点单元应力路径发展历程Fig.18 Development of stress paths at injection element (a)—φ1=10 %；(b)—φ1=20 %；(c)—φ1=30 %；(d)—φ1=40 %.

图 19给出了不同体积分数，两种方法计算的起裂压力.当包裹体体积分数相同时，三种随机工况均表明改进模型计算得到的起裂压力高于隐式模型，绝对误差最大值为0.101 MPa.基于线弹性理论的起裂压力理论解为p_b=3σ_h-σ_v+R_m= 116.962 MPa^[9]，与本文中隐式模型和改进模型计算结果的相对误差分别在8.98 % ~12.02 % 和8.94 % ~11.98 % 之间.
图 19(Fig. 19)

图 19 不同方法计算的起裂压力Fig.19 Break pressure calculated by different methods

5 结论1) 建立了基于柏林噪声的非规则几何边界的水平集方程，采用重心坐标法检测颗粒碰撞，实现以粒径、长短半轴比、倾角和填充率为参数的非规则颗粒填充算法.
2) 在已有LSM建模方法的基础上，采用均匀化方法建立考虑非规则颗粒细观边界特征的LSM-MT方法.通过与显式方法的计算结果对比，相对误差为0.03 % ~2.58 %，与LSM模型计算结果相比，改进幅度为0.12 % ~2.30 %，结果表明LSM-MT法可以有效考虑非规则颗粒细观边界特征.
3) 所建立的PHF-LSM-MT方法实现了非均质岩石材料建模及水力裂缝传播过程模拟，通过与线弹性理论的起裂压力对比，得到隐式模型和改进模型的相对误差分别在8.98 % ~12.02 % 和8.94 % ~11.98 % 之间，且随包裹体含量的增加，模型起裂压力降低.
参考文献

[1]	Zhang J, Li Y W, Pan Y S, et al. Experiments and analysis on the influence of multiple closed cemented natural fractures on hydraulic fracture propagation in a tight sandstone reservoir[J]. Engineering Geology, 2021, 281: 105981. DOI:10.1016/j.enggeo.2020.105981
[2]	Zhou Y, Yang D, Zhang X, et al. Numerical investigation of the interaction between hydraulic fractures and natural fractures in porous media based on an enriched FEM[J]. Engineering Fracture Mechanics, 2020, 235: 107175. DOI:10.1016/j.engfracmech.2020.107175
[3]	Linkov A M, Markov N S. Improved pseudo three-dimensional model for hydraulic fractures under stress contrast[J]. International Journal of Rock Mechanics and Mining Sciences, 2020, 130: 104316. DOI:10.1016/j.ijrmms.2020.104316
[4]	Liu C, Zhang J N, Yu H, et al. New insights of natural fractures growth and stimulation optimization based three-dimensional cohesive zone model[J]. Journal of Natural Gas Science and Engineering, 2020, 76: 103165. DOI:10.1016/j.jngse.2020.103165
[5]	Roostaei M, Taghipoor S, Nouri A, et al. Smeared modeling of hydraulic fracture using partially coupled reservoir and geomechanics simulators[J]. International Journal of Rock Mechanics and Mining Sciences, 2019, 113: 99-111. DOI:10.1016/j.ijrmms.2018.11.008
[6]	Sun Y, Liu Z Y, Tang X H. A hybrid FEMM-phase field method for fluid-driven fracture propagation in three dimension[J]. Engineering Analysis with Boundary Elements, 2020, 113: 40-54. DOI:10.1016/j.enganabound.2019.12.018
[7]	Yan C Z, Zheng H, Sun G H, et al. Combined finite-discrete element method for simulation of hydraulic fracturing[J]. Rock Mechanics and Rock Engineering, 2016, 49(4): 1389-1410. DOI:10.1007/s00603-015-0816-9
[8]	张国凯, 李海波, 夏祥, 等. 岩石细观结构及参数对宏观力学特性及破坏演化的影响[J]. 岩石力学与工程学报, 2016, 35(7): 1341-1352. (Zhang Guo-kai, Li Hai-bo, Xia Xiang, et al. Effects of microstructure and micro parameters on macro mechanical properties and failure of rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(7): 1341-1352. DOI:10.13722/j.cnki.jrme.2015.1154)
[9]	Li M, Guo P J, Dieter S, 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. DOI:10.1016/j.jngse.2020.103518
[10]	Skala V. Barycentric coordinates computation in homogeneous coordinates[J]. Computer & Graphics, 2008, 32: 120-127.
[11]	Meng Q X, Wang H L, Xu W Y, et al. Numerical homogenization study on the effects of columnar jointed structure on the mechanical properties of rock mass[J]. International Journal of Rock Mechanics and Mining Sciences, 2019, 124: 104127. DOI:10.1016/j.ijrmms.2019.104127
[12]	Perlin K. An image synthesizer[J]. Computer & Graphics, 1985, 19(3): 287-296.
[13]	Zirwes T, Zhang F, Wang Y Q, et al. In-situ flame particle tracking based on barycentric coordinates for studying local flame dynamics in pulsating Bunsen flames[J]. Proceedings of the Combustion Institute, 2021, 38: 2057-2066. DOI:10.1016/j.proci.2020.07.033
[14]	李明, 郭培军, 李鑫, 等. 基于水平集法的非均质岩石建模及水力压裂传播特性研究[J]. 岩土力学, 2016, 37(12): 3591-3607. (Li Ming, Guo Pei-jun, Li Xin, et al. Modelling method of heterogeneous rock material based on level set method and hydraulic fracture propagation features[J]. Rock and Soil Mechanics, 2016, 37(12): 3591-3607.)
[15]	Barral M, Chatzigeorgiou G, Meraghni F, et al. Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: theory and numerical validation[J]. International Journal of Plasticity, 2020, 127: 102632. DOI:10.1016/j.ijplas.2019.11.011
[16]	Nguyen B N, Hou Z S, Last G V, et al. Three-dimensional analysis of a faulted CO2 reservoir using an Eshelby-Mori-Tanaka approach to rock elastic properties and fault permeability[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2016, 8: 828-845. DOI:10.1016/j.jrmge.2016.06.007
[17]	Strazzeri V, Karrech A, Elchalakani M. Micromechanics modelling of cement stabilized rammed earth[J]. Mechanics of Materials, 2020, 148: 103540. DOI:10.1016/j.mechmat.2020.103540
[18]	Zhang H T, Zhu J M, Liu Y H, et al. Strength properties of jointed rock masses based on the homogenization method[J]. Acta Mechanica Solida Sinica, 2012, 25(2): 177-185. DOI:10.1016/S0894-9166(12)60018-4
[19]	Halpin J C, Kardos J L. Halpin-Tsai equations: areview[J]. Polymer Engineering and Science, 1976, 16(5): 344-352. DOI:10.1002/pen.760160512
[20]	Corey A T. The interrelation between gas and oil relative permeabilities[J]. Producers Monthly, 1954, 19: 38-41.