, 张砚琨, 梁力 东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2020-11-10
基金项目：中央高校基本科研业务费专项资金资助项目(N2001015)；国家留学基金委项目(201906085021)；国家自然科学基金资助项目(51474048)。
作者简介：李明(1980-)，男，辽宁沈阳人，东北大学副教授。
摘要：基于PHF(permeability-based hydraulic fracture)模型, 采用水平集方法(level set method, LSM)描述材料边界分布, 并用多尺度方法MT模型, 对材料有效参数进行均匀化, 建立适用于考虑非均质岩石材料细观界面特征分布与水力裂缝动态发展过程的PHF-LSM-MT数值计算模型.该模型解决了已有PHF-LSM模型无法精确描述积分点作用范围内细观尺度材料分布特征的问题.通过比较包裹体几何分布和体积分数影响下的有效弹性参数, 验证了采用MT模型均匀化细观特征的可行性; 在此基础上, 对非均质岩石水力裂缝发展过程进行研究, 得到了考虑细观界面特征情况下的裂缝发展过程, 探讨了等效开裂区域影响范围内关键位置的水压力变化特征和应力路径发展过程.
关键词：水力压裂岩石材料界面PHF-LSM均匀化方法
Numerical Modeling for the Material Interface of Heterogeneous Rock at Mesoscale and Hydraulic Fracturing Based on the PHF-LSM-MT
LI Ming, ZHAO Qi
, ZHANG Yan-kun, LIANG Li School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Qi, E-mail: zhaoqi_326326@163.com.
Abstract: The permeability-based hydraulic fracture(PHF) method is integrated with the level set method(LSM) and Mori-Tanaka(MT) model to study the hydraulic fracture propagation in heterogeneous rock(HR)with the mesoscale material interface being considered. In the proposed method, the PHF method is for simulating the hydraulic fracture propagation, and the LSM is for capturing the material interface with the mesoscale material properties of integration points being homogenized by MT model. To verify the mesoscale homogenization modeling procedure for HR, numerical simulations were carried out to estimate the effective elastic properties of HR with the volume fraction and inclusions' dimension being taken into account, and the results indicate that the proposed method is capable of homogenizing the material properties at mesoscale. Numerical simulations of hydraulic fracture propagation in HR with the mesoscale material interface show the characteristics of hydraulic fracture development, the features of pore pressures and stress paths of the equivalent fracture zone.
Key words: hydraulic fracturematerial interface of rockPHF-LSMhomogenization method
水力压裂技术是目前石油和天然气开采过程中的重要增产措施，自该技术提出以来, 在物理模型试验、理论研究和数值模拟等方面开展了深入和广泛的研究.例如, 基于经典断裂力学理论基础的水力压裂模型包括: 由Perkins, Kern和Nordgren提出的PKN模型, 由Khristianovitch, Geertsma和de Klerk等提出的KGD模型, 准三维水力压裂理论模型P-3D模型, 在此基础上通过引入Carter方程考虑流体滤失的PKN-C模型, KGD-C模型和P-3D-C模型等[1].虽然用于水力压裂计算的理论模型能得出裂缝的主要特征, 但是其理论解均建立在严格的假设基础上, 对于求解复杂地质条件下的水力压裂发展过程有很大限制; 而数值计算方法可以解决这一问题, 同时与试验方法相比有成本低、效率高的优点.
水力压裂计算的数值方法有弥散裂缝模型[2]、内聚力模型[3]、扩展有限元法[4]和相场法[5]等.其中弥散裂缝模型在计算过程中无需引入真正的裂缝和重构网格, 适用于岩石材料拉伸破坏和土体材料剪切破坏过程的模拟.该方法通过修改破坏区域的材料参数(如弹性模量和渗透系数等)实现水力裂缝发展过程模拟[6].同时该模型易于与其他数学方法结合, 例如: 以弥散裂缝模型为基础的PHF(permeability-based hydraulic fracture)水力压裂数值计算模型假设岩石材料在开裂过程中的渗透率为平均有效应力的函数, 基于流固耦合理论, 实现了水力裂缝在均质岩石材料中的传播过程的模拟[7]; 通过引入水平集法(level set method, LSM)[8], 采用隐式方法描述具有不同材料分布特征的非均质岩石材料的静态边界, 避开了采用弥散裂缝模型进行水力压裂计算时材料分布对有限元网格依赖的问题, 进而建立PHF-LSM模型, 该模型实现了非均质岩石的建模及水力压裂过程的模拟.已有研究表明, 采用PHF-LSM模型可建立含有规则和非规则椭圆形硬包裹体分布、具有分层特性且材料空间分布满足韦伯分布的非均质岩石材料模型, 以及建立含有多种类型节理随机分布的岩体模型, 并能够分析水力裂缝特征、起裂压力和多点压裂过程等相关问题[9-11].虽然PHF-LSM模型可以建立多种非均质岩石模型, 但该模型未能精确考虑材料边界的细观特征.地下储层的岩石多具有多孔性质并伴有夹杂物、初始裂隙、节理等, 同时油气开采的工程问题又涉及大尺度的力学模型; 当需要考虑岩石细观结构对水力压裂过程的影响时, 可通过引入多尺度方法来解决这一问题[12-15].
在岩石中通常分布有硬包裹体(如巴拉波石英岩、石灰岩等)和软包裹体(石英云母片岩等), 本文针对具有硬包裹体分布特征的非均质岩石材料, 在已有的PHF-LSM模型基础上, 引入均匀化方法建立能够考虑非均质岩石材料细观界面分布特征的PHF-LSM-MT水力压裂数值计算模型; 采用该模型建立非均质岩石材料力学模型, 进而分析水力裂缝动态传播过程及水压力发展特性等.
1 非均质岩石细观界面模型PHF-LSM法[10-11]在进行有限元建模时, 积分点材料属性的取值只与积分点位置的材料性质有关, 如图 1中包含材料界面单元的积分点i1和i5的材料属性均会被赋值为基岩, 而积分点i2和i6的材料属性会被赋值为包裹体.显然该模型没有考虑积分点周围材料分布特征的影响.本文引入均匀化方法对包含边界单元的积分点进行材料参数均匀化, 进而考虑材料界面细观特征的影响.
图 1(Fig. 1)
 | 图 1 基于LSM建模方法的单元特性示意图Fig.1 Sketch of modeling element characteristics by LSM |
1.1 非均质岩石材料边界的水平集方程以边长为l、含有大小和方向均随机变化的椭圆形包裹体的非均质岩石材料为主要研究对象, 模型受到水平和竖向地应力(σh和σv)作用的同时四周布置初始水压力(pw), 如图 2所示, 采用水平集法(隐式方法)描述有限元模型中该类包裹体分布特征, 其对应的水平集方程见式(1)和式(2)[9].
 | (1) |
 | (2) |
 | 图 2 含随机分布的椭圆形包裹体的岩石材料模型Fig.2 Rock sample with randomly distributed elliptical inclusions |
式中: Li为第i个包裹体, x为求解域Ω内的点向量; 平面问题对应的椭圆水平集方程为?i(x, y), 其中αa和αb分别为半轴ai和bi对应的随机变量; (xi0, yi0)为由满足均匀分布的随机数生成的椭圆包裹体中心点坐标, 对于平面问题还需引入随机取值的倾角βi.
?i < 0和?i> 0分别表示当前材料点处于包裹体内部和包裹体外部(基岩).当材料积分点所在区域(w×h)内包含材料边界时, 已有的简化方法是通过阶跃函数、线性函数或非线性函数对两种材料交界处进行处理[9], 但无法考虑材料分布的影响.本文通过引入均匀化方法考虑有效材料参数与包裹体分布的关系:
 | (3) |
1.2 有效材料参数均匀化为使图 1所示的含有材料界面单元中的积分点(如i1, i5, i2和i6)能够有效反映局部范围材料分布特征的影响, 采用Mori-Tanaka(MT)法[16-17]对该局部区域进行材料弹性参数(E, υ)的均匀化:
 | (4) |
 | (5) |
对于有些问题直接获得多相材料精确的有效参数是非常困难的, 因此可采用近似的材料参数界限.具有代表性的界限模型有Voigt上限(VR+)和Reuss下限(VR-), 分别如式(6)和式(7)所示[17].
 | (6) |
 | (7) |
计算中, 饱和岩石材料中的有效渗透系数κ采用假设包裹体为椭圆形纤维的Halpin-Tsai(HT)模型计算[18], 如式(8)和式(9)所示.
 | (8) |
 | (9) |
当积分点作用域内含有不同分布的两种或多种材料时, 其材料抗拉强度(Rm)的计算较为复杂.在已有的PHF-LSM方法中, Rm仅与积分点处的材料性质有关.如果Rm直接取值为作用范围内的最小值或最大值，则会相应减小或增大包裹体的边界范围, 因此在本文的数值计算中简化取为抗拉强度的Voigt上限值[19].
1.3 非均质岩石有限元建模及有效弹性参数由于在分析中采用隐式方法描述材料边界, 因此可基于相同有限元网格建立不同的力学模型[11].在进行积分点位置材料参数均匀化时, 为使算法能处理更具有一般性的非规则材料边界, 本文直接采用离散方法首先将积分点作用区域w×h离散为m×n个子网格, 然后采用水平集方程(1)和(2)对子网格进行边界检测, 进而均匀化材料参数.本文计算中取m=n=100, 分别采用MT法、HT法和Voigt法计算该区域的有效弹性参数(弹性模量和泊松比)、有效渗透系数和有效抗拉强度, 如图 3所示.子区域宽度w和高度h取值分别为单元宽度和高度的1/3.
图 3(Fig. 3)
 | 图 3 积分点位置材料参数均匀化Fig.3 Homogenization of material properties at integration point |
当边长为l的岩石模型中包含单个或多个椭圆形包裹体时, 对应的几何参数(包裹体半轴a和b长度、a/b比值、倾角β和体积分数φ)取值范围如表 1所示.采用图 4a所示的有限元网格, 建立如图 4b~图 4e所示的非均质岩石材料的有限元模型, 其中单元高宽比AR=1.61[20].
表 1(Table 1)
 表 1 椭圆形包裹体几何参数取值Table 1 Geometry parameters of elliptical inclusions
| 表 1 椭圆形包裹体几何参数取值 Table 1 Geometry parameters of elliptical inclusions |
图 4(Fig. 4)
 | 图 4 不同参数下椭圆形包裹体有限元模型Fig.4 Finite element model with elliptical inclusions with different modeling parameters (a)—有限元网格(AR=1.61); (b)—单包裹体; (c)—a/b=0.6;(d)—a/b=1.0; (e)—a/b=0.2~1.0. |
以垂直方向y的有效弹性模量Ey为例, 共考虑3种建模方法对计算结果的影响: 传统有限元法模型(显式方法)、基于LSM的隐式有限元模型(隐式方法)和引入均匀化方法后的修正隐式有限元模型(修正方法).计算中基岩和包裹体的弹性模量分别取值20 GPa和80 GPa, 泊松比分别取值0.25和0.20.
当岩石模型中含有单个包裹体, 且体积分数不变时, 包裹体角度β对Ey的影响如图 5a所示.计算结果表明, 对于不同体积分数和不同a/b取值, 当β小于30°时, Ey的增加量为0.23%~0.67%, 受角度β的影响较小; 当角度β大于30°时, Ey的增加量为3.31%~5.19%, 受角度β的影响显著提高; 同时三种方法计算的结果较为接近.当基岩含有多个随机分布的包裹体时, 不同体积分数条件下随机算例的Ey变化如图 5b所示.计算结果表明, 在相同体积分数下, 包裹体的随机分布影响了垂直方向的有效弹性模量; 不同方法计算结果均在VR上下界限内, 在理论解(MT方法)上下波动.隐式方法和修正方法对显式方法的相对误差er采用式(10)计算:
 | (10) |
 | 图 5 不同建模方法垂直方向有效弹性模量对比Fig.5 Comparison on vertical effective elastic modulus with different modeling methods |
式中: Πim为隐式方法或修正方法的计算结果; Πex为显式方法计算结果.
以a/b为固定值和随机值的四种类型多包裹体分布为例, 当体积分数φ1在0~1.0范围变化时, 隐式方法和修正方法计算结果的相对误差如图 6所示.从计算结果可以看出, 修正方法的误差小于隐式方法, 即采用MT法均匀化弹性参数后的修正方法计算结果更接近于显式方法, 进而表明均匀化方法的引入有效地考虑了材料边界的影响.
图 6(Fig. 6)
 | 图 6 多包裹体随机分布模型的有效弹性模量相对误差Fig.6 Relative errors of effective elastic modulus of the rock sample with randomly distributed inclusions (a)—a/b=0.4; (b)—a/b=0.8; (c)—a/b=1.0; (d)—a/b=0.2~1.0. |
2 非均质岩石水力压裂计算模型在基于流固耦合理论并假设岩石材料开裂区域内渗透率为平均有效应力函数的PHF水力压裂数值计算模型基础上, 采用LSM法描述椭圆形硬包裹体几何分布，采用MT等均匀化方法计算考虑材料积分点影响范围内边界细观特征的有效材料参数, 建立用于分析非均质岩石材料水力裂缝传播过程的数值计算模型.
2.1 PHF水力压裂数值计算模型PHF模型[7]中饱和岩石材料的固体骨架和孔隙流体的平衡方程分别如式(11)和式(12)所示.其中应变的定义如式(13)所示.计算中采用线弹性理论描述岩石材料开裂前的力学行为, 见式(14);岩石材料开裂后采用摩尔-库伦塑性理论进行描述.计算中的有效应力原理由式(15)给出.渗透系数和渗透率的关系如式(16)所示.
 | (11) |
 | (12) |
 | (13) |
 | (14) |
 | (15) |
 | (16) |
式(17)和式(18)分别为开裂区域渗透系数由平均有效应力引起的改变和压裂液黏度特性对开裂区域渗透系数的影响.
 | (17) |
 | (18) |
2.2 PHF-LSM-MT水力压裂计算模型结合ABAQUS软件提供的二次开发功能实现PHF-LSM-MT水力压裂数值计算模型, 流程图如图 7所示.计算过程中共设置8个场变量, 分别为有效弹性模量(E)、有效泊松比(υ)、有效抗拉强度(Rm)、有效渗透系数(κ)、密度(ρ)、饱和度(Sw)、孔隙比(e)和混合黏滞系数(μmix); 各参量的初值存储在对应的8个状态变量中.涉及的三个主要用户子程序为UEXTERNALDB, SDVINI和USDFLD.在ABAQUS求解之前, 采用PYTHON程序实现椭圆填充并计算包裹体的随机位置分布，以及进行几何特征参数的前处理, 并存储于数据文件中; 子程序UEXTERNALDB用于实现求解之前的岩石材料细观特征模型参数读取; 子程序SDVINI通过LSM法实现积分点所处位置坐标的判断, 并采用均匀化方法(式(4)~式(9))计算有效材料参数, 以及初始化8个状态变量; 子程序USDFLD实现PHF水力压裂数值计算模型(式(17)和式(18)).
图 7(Fig. 7)
 | 图 7 PHF-LSM-MT方法流程图Fig.7 Flow chart for the PHF-LSM-MT method |
3 水力裂缝动态发展特性针对图 2所示力学模型, 建立考虑细观特征分布的非均质岩石材料模型.考虑包裹体体积分数为10%, 30%, 50%的三种工况, 每种工况包括三个随机算例.所建立的数值计算模型如图 8所示.
图 8(Fig. 8)
 | 图 8 基于均匀化方法的非均质岩石有限元模型Fig.8 Finite element models of heterogeneous rock sample based on homogenization method (a)(b)(c)—φ1=10%; (d)(e)(f)—φ1=30%; (g)(h)(i)—φ1=50%. |
从图 8可以看出, 当采用隐式方法描述材料边界时, 由于只考虑积分点所在位置的材料属性, 故会出现锯齿形状; 而均匀化方法的引入修正了原有隐式模型对材料边界描述的不足, 如图 8a所示.
针对采用修正方法建立的非均质岩石模型进行水力压裂数值模拟试验, 模型受到的垂直和水平方向初始应力分别为70和57 MPa, 初始水压力为49 MPa.注水时间和注水速率分别为1 200 s和0.001 mm/s.计算中假设包裹体不易开裂, 并赋予较高的抗拉强度, 相应的材料参数取值如表 2所示.
表 2(Table 2)
表 2 非均质岩石材料力学参数Table 2 Material parameters of heterogeneous rock sample
| 表 2 非均质岩石材料力学参数 Table 2 Material parameters of heterogeneous rock sample |
当注水时间为1 200 s时, 不同包裹体含量下的水力裂缝最终形态如图 9所示, 可以看出裂缝发展方向均沿着大主应力方向(垂直方向).同时水力裂缝的发展也受到了包裹体分布的影响: ①当包裹体与裂缝发展方向接触面较小时, 等效开裂区域会绕过包裹体继续沿大主应力方向传播, 见图 9b; ②当接触面较大, 且与裂缝传播方向夹角较大时, 裂缝会沿着包裹体边界发展, 见图 9a和图 9c; ③当夹角较小时, 裂缝会被阻止进一步延伸, 同时沿着横向发展直到可绕过包裹体, 见图 9d、图 9e和图 9i.以图 9f为例, 其裂缝动态发展过程如图 10b~图 10f所示.
图 9(Fig. 9)
 | 图 9 不同岩石模型的水力裂缝(t=1 200 s)Fig.9 Hydraulic fractures of different rock samples (a)(b)(c)—φ1=10%; (d)(e)(f)—φ1=30%; (g)(h)(i)—φ1=50%. |
图 10(Fig. 10)
 | 图 10 水力裂缝动态发展过程Fig.10 Dynamic development process of hydraulic fracture (a)—关键点位置; (b)—t=98 s; (c)—t=226 s; (d)—t=708 s; (e)—t=948 s; (f)—t=1 200 s. |
各随机算例(如图 9所示)注水点位置水压力发展过程和应力路径分别如图 11和图 12所示.由图 11可以看出, 注水点位置的水压力先增加至起裂压力后降低至传播压力, 该发展特征与已有计算结果一致[10].包裹体分布影响了水压力的传播过程, 由水压力发展曲线可得各算例的起裂压力Pb在100.1~134.5 MPa之间, 与基于线弹性理论的起裂压力解(3σmin-σmax+Rm =108 MPa)[22]的相对误差在3.3%~24.5%之间.包裹体的分布也同时影响了注水点位置的应力路径发展过程, 如图 12所示.结果表明, 应力路径起点(初应力状态)受包裹体分布的影响而不同; 随着注水时间增加, 当岩石处于受压状态时, 平均有效应力p逐渐降低至0;随着注水时间的进一步增加, 注水点位置由受压状态转变为受拉状态后, 平均有效应力p开始增加并达到岩石材料的抗拉强度, 产生水力裂缝; 在整个注水过程中, 等效剪应力q先增加后沿着屈服面降低至0.
图 11(Fig. 11)
 | 图 11 注水点水压力发展过程Fig.11 Pore pressure development of injection points 图中a~i对应图 9a~图 9i |
图 12(Fig. 12)
 | 图 12 注水点应力路径发展过程Fig.12 Stress path development of injection points 图中a~i对应图 9a~图 9i |
除注水点以外, 水力裂缝发展路径上的点可分为三类: 裂缝经过的基岩内点A2~A10, 包裹体边界点B1~B5, 以及包裹体内部点C1~C3, 如图 10a所示.三类点对应的水压力随注水时间发展过程和应力路径变化分别如图 13和图 14所示.由于水力压裂过程为张拉破坏, 故开裂区域经过点的水压力为先降低再增加至传播压力.而包裹体内部点无裂缝经过(未开裂), 其水压力增加较小; 距离开裂区域较近的包裹体(点C1和C2)水压力增加高于较远处包裹体(点C3); 点C4距离裂缝较远, 故水压力维持在初始状态.
图 13(Fig. 13)
 | 图 13 关键点水压力发展Fig.13 Pore pressure development of key points |
图 14(Fig. 14)
 | 图 14 关键点应力路径Fig.14 Stress path development of key points |
等效开裂区域内点(A5, A10, B4和B5)的平均有效应力和剪应力在裂尖到达前逐渐增加, 开裂过程中逐渐降低, 由受压状态转至受拉状态直到开裂.因包裹体边界点(A10, B4和B5)所在单元包含两种材料, 故均匀化的抗拉强度高于基岩(如点A5处Rm=7 MPa)且低于包裹体(如点C1处Rm=40 MPa), 因此该类点应力路径的屈服面受抗拉强度影响而不同, 其发展过程如图 14所示.该计算结果表明: 当水力裂缝到达该类点(如点A5和A10)之前, 平均有效应力p和剪应力q均有不同程度的增加; 当裂缝经过该点时, p逐渐降低至0后由受压状态转为受拉状态, p开始进一步增加至抗拉强度后产生裂缝; q先增加后沿着屈服面降低至0.由于假设包裹体内部点不易开裂, 并赋予较高的抗拉强度(40 MPa), 故该类点(如点C1和C3)在水力裂缝的传播过程中处于受压状态未开裂, 其平均有效应力p在裂缝传播过程中维持不变, 而剪应力q增加.
4 结论1) 在已有的基于LSM的隐式有限元建模方法基础上, 通过引入均匀化方法考虑材料边界所在单元积分点位置的有效材料参数, 得到非均质岩石材料有限元建模的改进方法.
2) 将改进的建模方法与PHF水力压裂数值计算模型结合, 建立了可用于考虑岩石材料细观边界特征与水力压裂过程模拟的PHF-LSM-MT法.
3) 采用PHF-LSM-MT模型计算了不同体积分数下的非均质岩石材料的水力裂缝传播过程, 得到了注水点、等效开裂区域内部点、包裹体边界和内部点的水压力和应力路径发展特征.
参考文献
| [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] | Hu Y J, Chen G L, Cheng W P, et al. Simulation of hydraulic fracturing in rock mass using a smeared crack model[J]. Computer and Structure, 2014, 137: 72-77. DOI:10.1016/j.compstruc.2013.04.011 |
| [3] | Zhai L H, Zhang H, Pan D B, et al. Optimisation of hydraulic fracturing parameters based on cohesive zone method in oil shale reservoir with random distribution of weak planes[J]. Journal of Natural Gas Science and Engineering, 2020, 75: 103130. DOI:10.1016/j.jngse.2019.103130 |
| [4] | Saberhosseini S E, Ahangari K, Mohammadrezaei H. Optimization of the horizontal-well multiple hydraulic fracturing operation in a low-permeability carbonate reservoir using fully coupled XFEM model[J]. International Journal of Rock Mechanics and Mining Sciences, 2019, 114: 33-45. DOI:10.1016/j.ijrmms.2018.09.007 |
| [5] | 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 |
| [6] | 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 and Engineering, 2016, 35: 761-774. DOI:10.1016/j.jngse.2016.09.018 |
| [7] | Li M, Guo P J, 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 |
| [8] | 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. DOI:10.1016/0021-9991(88)90002-2 |
| [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 and Engineering, 2017, 157: 409-421. DOI:10.1016/j.petrol.2017.07.017 |
| [10] | Li M, Guo P J, Stolle D, et al. Modeling method and hydraulic fracture propagation for jointed rock mass using PHF-LSM method[J]. Journal of Natural Gas Science and Engineering, 2019, 68: 102910. DOI:10.1016/j.jngse.2019.102910 |
| [11] | Li M, Guo P J, Stolle D, et al. Heterogeneous rock modelling 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 |
| [12] | Klusemann B, B?hm H J, Svendsen B. Homogenization methods for multi-phase elastic composites with non-elliptical reinforcements: comparisons and benchmarks[J]. European Journal of Mechanics A: Solids, 2012, 34: 21-37. DOI:10.1016/j.euromechsol.2011.12.002 |
| [13] | Will J, Eckardt S, Ranjan A. Numerical simulation of hydraulic fracturing process in an enhanced geothermal reservoir using a continuum homogenized approach[J]. Procedia Engineering, 2017, 191: 821-828. DOI:10.1016/j.proeng.2017.05.249 |
| [14] | 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 |
| [15] | Cao Y J, Shen W Q, Shao J F, et al. Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale[J]. European Journal of Mechanics A: Solids, 2020, 81: 103958. DOI:10.1016/j.euromechsol.2020.103958 |
| [16] | Mori T, Tanaka K. Averages stress in matrix and average elastic energy of materials with misfitting inclusions[J]. Acta Metallurgica et Materialia, 1973, 21(5): 571-574. DOI:10.1016/0001-6160(73)90064-3 |
| [17] | Qu J M, Cherkaoui M C. Inclusions and inhomogeneities[M]// Fundamentals of Micromechanics of Solids. [S. l. ]: John Wiley & Sons, 2006. |
| [18] | Halpin J C, Kardos J L. Halpin-Tsai equations: a review[J]. Polymer Engineering and Science, 1976, 16(5): 344-352. DOI:10.1002/pen.760160512 |
| [19] | 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 |
| [20] | 李明, 梁力, Guo Pei-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) |
| [21] | Corey A T. The interrelation between gas and oil relative permeabilities[J]. Producers Monthly, 1954, 19: 38-41. |
| [22] | Yew C H, Weng X W. Mechanics of hydraulic fracturing[M]. [S. l. ]: Gulf Professional Publishing, 1997. |
