基于细观损伤多相耦合的砂砾岩水力压裂裂缝扩展数值模拟

余东合徐康泰车航张登文刘国华马新仿

1. 中国石油华北油田分公司采油工程研究院； 2. 中国石油大学（北京）石油工程学院

基于细观损伤多相耦合的砂砾岩水力压裂裂缝扩展数值模拟

余东合1徐康泰2车航1张登文1刘国华1马新仿2

1. 中国石油华北油田分公司采油工程研究院； 2. 中国石油大学（北京）石油工程学院

低渗透砂砾岩油藏水力压裂裂缝扩展机理及其数值模拟研究，对该类储层压裂改造成功实施具有重要意义。将砂砾岩储层中砾石表征为基质-交界面-砾石的三模态结构，假定砾石分布与几何尺寸及储层物性满足随机分布，结合Moter-Carolo方法，完成砂砾岩储层数学表征；考虑储层渗流场、应力场、水化膨胀湿度场的三相耦合特征，结合损伤力学、断裂力学等原理，利用细观损伤有限元的方法，建立了砂砾岩储层水力压裂裂缝扩展数学模型，并进行数值模拟研究。模拟分析了不同主应力差、基质-砾石交界面强度、砾石强度情况下，水力裂缝遇砾石扩展情况，并最终实现砂砾岩储层水力裂缝动态扩展数值模拟。研究表明，水力裂缝遇砾发生绕砾、穿砾、止裂现象，并以绕砾扩展为主，且裂缝发生明显转向，存在羽状次生裂缝；裂缝转向程度和裂缝延伸长度与主应力差、砾石强度以及交界面强度有关，主要表现有：水平主应力差越小，水力裂缝遇砾转向越明显；基质-砾石交界面强度增加，水力裂缝明显变短，并难以转向；随着砾石强度的增大，裂缝的转向程度增大。

砂砾岩储层；水力裂缝扩展；多相耦合；细观损伤；数值模拟

乌里雅斯太地区低渗透砂砾岩储层油藏开发难度大，具有油层物性差、非均质性强、强水敏性等特征，压裂改造是砂砾岩储层开发的主要技术［1］。砾石的存在使得裂缝扩展形态难以控制，水力裂缝扩展机理不明确，压裂改造难度大。国外对砂砾岩储层水力裂缝扩展问题研究较少，国内学者主要通过室内真三轴物模实验和数值模拟方法对其进行初步研究：赵益忠等［2］对不同岩性岩样进行水力压裂实验，分析了不同岩样对水力压裂裂缝扩展的影响；王昊［3］利用人造岩样添加砾石，研究砂砾岩油藏砾石对压裂裂缝延伸的影响，得出砂砾岩水力压裂裂缝扩展主要以发生绕砾现象为主；李连崇等［4］利用RFPA原理，将其运用到砂砾岩水力压裂裂缝扩展规律研究中，分析砂砾岩水力压裂裂缝形态；李根等［5］利用RFPA方法对砂砾岩储层水力裂缝扩展进行数值模拟分析，基于国内某砂砾岩储层相关参数，分析了水力裂缝扩展情况与规律；罗攀等［6］进行了砂砾岩水力裂缝延伸路径模拟研究，通过编程模拟了砾径、含量以及断裂韧性对砂砾岩储层水力裂缝扩展和压力的影响。以上学者对于砂砾岩储层水力压裂裂缝扩展数值模拟研究均是在RFPA软件的平台上进行，并没有一套专门适用于分析砂砾岩储层裂缝扩展的方法和程序。笔者建立复杂介质条件下的裂缝扩展模型，基于损伤力学、细观有限元等相关方法［7-10］，编制相应程序，探索低渗砂砾岩储层水力压裂裂缝扩展规律和裂缝形态，进一步完善砂砾岩油藏增产改造机理。

1　砂砾岩裂缝扩展数学模型的建立

Building of mathematical model of fracture propagation in glutenite reservoir

1.1砂砾岩储层砾石表征

Gravel characterization of glutenite reservoir将砂砾岩储层中砾石表示为“基质-交界面-

砾石”三模态圆形结构［1，11-12］，考虑砾石分布和几何尺寸的随机分布性［11-13］，完成砂砾岩储层数学表征。（1）砾石中心点坐标（x0，y0）。设中心点位置（x0，

y0）服从均匀分布，生成砾石的区域为W×H，区域的左下角为坐标原点（0，0），rand（0，1）为（0，1）内的伪随机数，则

式中，W为研究区域长，m；H为研究区域宽，m。

（2）生成粒径l。设砾石粒径服从正态分布，其相应的概率密度函数为

式中，μ为l的均值，m；σ为l的均方差，m。

在假定了砾石上述几何参数的先验概率模型后，采用 Monte-Carlo 模拟实现参数的随机抽样，对于任意给定的分布函数F（x），直接抽样法如下：

式中，XF为由已知分布函数F（x）所产生的简单子样X1，X2，…，XN中的个体；ξ为伪随机序列ξ1，ξ2，…，ξN中的对应伪随机数。

1.2砂砾岩储层非均质性表征

Heterogeneity characterization of glutenite reservoir

由于砾石的存在，砂砾岩储层一般具有强非均质性［13-14］，考虑储层物性参数满足一定的统计规律，采用Weibull分布，完成储层非均质性表征。

式中，α为具有一定统计规律的相关参数；α0为参数α的平均值；m为性质参数，表示参数α的均匀程度；φ（α）为参数α的统计学分布密度。

以弹性模量为例，式（4）写成弹性模量形式

式中，E'为弹性模量平均值，GPa；φ（E）为弹性模量基元分布值。

1.3砂砾岩储层水力裂缝扩展损伤计算

Damage calculation for hydraulic fracture propagation in glutenite reservoir

对于砂砾岩储层水力裂缝的起裂和扩展，考虑其是在剪应力和张应力的共同作用下形成的［14-16］，即单元破坏形式为Mohr-Coulomb 剪切破坏与最大拉伸破坏相结合的方法。当单元的应力状态或者应变状态将满足某个给定的损伤阈值时，单元开始损伤。本文模型采用单元的“Mohr-Coulomb 剪切破坏”和“拉伸破坏”两种损伤阈值确定网格断裂问题。

（1） Mohr-Coulomb 剪切破坏。当单元剪应力达到Mohr-Coulomb 损伤阈值时，

式中，φ为摩擦角，°；fc为单轴抗压强度，MPa；σ1为最大主应力，MPa；σ3为最小主应力，MPa。

（2）最大拉伸破坏。当单元达到抗拉强度ft损伤阈值时

按照上述理论，首先采用载荷连续加载的方式对每一步给定的载荷增量进行应力计算，然后根据破坏准则来检查模型中是否有破坏单元（判断单元是否满足Mohr-Coulomb 剪切破坏或是最大拉伸破坏，并以最大拉伸破坏为主，在不能满足的条件下判断是否满足Mohr-Coulomb 剪切破坏，即将Mohr-Coulomb 剪切破坏作为第二破坏准则）。

当单元应力状态达到临界值时，单元发生破坏。单元的破坏模型主要有损伤判断准则和损伤演化方程。这里只给出弹性模量变化的情况，弹性模量的线性损伤演化方程为

式中，E0为无损伤单元的弹性模量，MPa；E为损伤单元的弹性模量，MPa；B为损伤因子。

1.4砂砾岩储层水力裂缝扩展数学模型

Mathematical model of hydraulic fracture propagation in glutenite reservoir

根据线弹性断裂理论，由水化膨胀湿度场方程、渗流方程和应力平衡方程［16-19］，建立砂砾岩储层水力压裂裂缝扩展的数学模型。

（1）砂砾岩储层水化膨胀湿度场方程。

式中，Cf为水的吸附扩散常数；W为含水率；x，y为裂缝扩展方向。

（2）渗流方程。

式中，k为渗透率，mD；p为流体压力，MPa。

（3）平衡方程。

式中，σij为应力，MPa；bi为体积力，MPa；i，j为张量下标，i=1，2，j=1，2。

（4）几何方程。

式中，εij为应变；u为位移，m。

（5）本构方程。

式中，σij'为有效应力，MPa；σ0为初始地应力，MPa；δij为Kronecker常数；εw为水化膨胀应变；D为弹性矩阵。

（6）边界条件和初始条件。

①水化膨胀湿度场边界条件：研究区域固定边界岩土含水恒定。

②渗流场边界条件：研究区域固定界限上设定为初始地层压力；裂缝扩展端点处的单元应力设定成延伸压力；考虑储层具有对称性，模拟1/2储层，且此对称边界处不考虑渗流速度。

③位移约束：不考虑研究区域内固定边界上位移以及对称边界上法向位移。

④初始条件：压裂施工前，按照σH、σh的值给定原始应力场，固定边界湿度和初始压力均为上述给定的未压裂前湿度（含水率）和压力。

2　水力裂缝遇砾形态数值模拟研究

Numerical simulation on hydraulic fracture shape when encountering gravel

根据目标区块低渗砂砾岩储层特征，考虑渗流场-水化膨胀湿度场-应力场三相耦合问题，建立相应数学方程，依次求解湿度场、渗流场、力位移方程，并进行单元断裂判断与损伤计算，完成数值模拟。

本节研究砂砾岩储层不同主应力差、不同基质-砾石交界面强度、不同砾石强度下，砾石存在对水力裂缝扩展的影响。建立单颗砾石基本模型，模型尺寸1 000 mm×1 000 mm，单元数100×100；设定初始裂缝长度位于x轴，长120 mm；砾石粒径为250 mm，砾石中心点位于x轴，且距初始裂缝左端点为350 mm；设定地层四周围压为30 MPa以控制渗流边界；初始压力设为25 MPa，每步增量0.05 MPa；最大、最小主应力大小分别为35 MPa、25 MPa；地层压力20 MPa；砾石强度平均为30 MPa，基质-砾石交界面强度平均为10 MPa；其余参数见表1。

表1　模型物理力学参数Table 1 Physical mechanical parameters of model

2.1不同主应力差

Different principal stress differences

利用上述建立的基本模型，以表1作为主要输入参数，研究主应力差对水力压裂遇砾形态的影响，选取不同主应力差为0、5、10、15 MPa，观察水力裂缝遇砾扩展形态，模拟结果见图1。

图1　不同主应力差下水力裂缝绕砾形态示意图Fig. 1 Gravel bypassing shape of hydraulic fracture under different principal stress differences

由图1可看出：（1）主应力差0 MPa时，水力裂缝遇砾后发生明显偏转，导致水力裂缝沿着转向方向继续延伸，没有观察到恢复初始延伸方向的趋势；随着主应力差的增大，裂缝遇到砾石后恢复沿最大主应力方向扩展的时间缩短，即转向延伸距离减小；（2）随着主应力差的增大，在相同条件下，裂缝整体延伸的距离变大；但随着应力差进一步增大，由于裂缝延伸压力变大，难度加大，裂缝扩展距离减小。

2.2砾石强度

Matrix-gravel interface strength

利用上述基本模型，选取砾石强度为20、30、40、50 MPa，观察水力裂缝遇砾扩展形态，结果见图2。可以看出：（1）当砾石强度较小时，水力裂缝将全部或者部分穿透砾石继续延伸，如模拟的第1种情况所示，当砾石强度与基质强度一样为20 MPa时，水力裂缝完全穿透砾石延伸，裂缝不发生偏转，此时可看作常规低渗透储层水力压裂情况；（2）随着砾石强度的增加，水力裂缝遇到砾石发生偏转程度增大。

2.3基质-砾石交界面强度

Gravel bypassing shape of hydraulic fracture under different gravel strengths

利用上述基本模型，选取基质与砾石交界面强度为5、10、15、20 MPa，观察水力裂缝遇砾扩展形态，模拟结果见图3。

图2　不同砾石强度下水力裂缝绕砾形态示意图Fig. 2 Gravel bypassing shape of hydraulic fracture under different interface strengths

图3　不同交界面强度下水力裂缝绕砾形态示意图Fig. 3 Gravel bypassing shape of hydraulic fracture under different gravel strengths

由图3可以看出：（1）随着基质与交界面强度的增加，水力裂缝明显变短，并难以转向；交界面强度5 MPa时，裂缝的长度最长，裂缝转向程度最弱；（2）随着基质与交界面强度的增加，当其强度达到基质强度时，水力裂缝发生止裂的现象，即为20 MPa时的模拟结果，这是由于交界面强度足够大，水力裂缝难以穿透或是绕过砾石继续延伸；（3）整体来看，水力裂缝遇砾转向及延伸主要受基质与交界面强度的影响，即基质-交界面强度是水力裂缝遇砾转向和扩展的主导因素。

3　砂砾岩储层水力裂缝扩展数值模拟实现

Numerical simulation realization of hydraulic fracture propagation in glutenite reservoir

在本文建立的砂砾岩储层水力压裂裂缝扩展数学模型的基础上，进行数值模拟实现。建立砂砾岩储层基本模型，模型尺寸：1 000 mm×1 000 mm；单元数：200×200；设定砾石粒径为15 mm，砾石个数为300个，砾石分布及几何特征满足本文1.1节随机分布函数，模型其余参数同第2节单砾模型。

根据建立的模型，得到如图4所示的裂缝扩展形态，可以看出：（1）由于储层的非均质性及砾石的存在，裂缝出现“羽状”结构，这种羽状次生裂缝的生成，加大了砂砾岩储层裂缝的复杂性，同时也验证了砂砾岩储层水力压裂施工过程中压力变化不稳定的情况；（2）水力裂缝存在较多分支缝，且水力裂缝遇砾后主要以绕砾扩展为主，存在穿砾与止裂的现象，水力裂缝发生明显转向。

图4　砂砾岩储层水力压裂裂缝扩展形态Fig. 4 Shape of hydraulic fracture propagation in glutenite reservoir

4　结论

Conclusions

（1）应用细观损伤有限元的方法，结合砂砾岩储层特征，包括砾石表征、非均质性表征、多相耦合效应等，建立了砂砾岩储层水力压裂裂缝扩展数学模型，该模型可用于分析砂砾岩储层水力裂缝动态扩展规律。

（2）裂缝转向程度和裂缝延伸长度与主应力差、砾石强度以及交界面强度有关，主要表现有：主应力差越大，裂缝遇到砾石后恢复沿最大主应力方向扩展的时间越短，即转向延伸距离越小，水力裂缝遇砾转向越不明显；随着砾石强度的增加，水力裂缝遇到砾石发生转向程度增大；当砾石强度较小时，水力裂缝将全部或者部分穿透砾石继续延伸；随着基质与交界面强度的增加，水力裂缝明显变短，并难以转向，当其强度达到基质强度时，水力裂缝发生止裂的现象，整体来看，基质-砾石交界面强度是水力裂缝遇砾转向与扩展的主导因素。

（3）砂砾岩储层水力裂缝不同于常规低渗储层，水力裂缝不是沿最大主应力方向双翼延伸，存在明显转向和羽状次生裂缝；遇砾发生绕砾、穿砾、止裂现象，其中以绕砾扩展为主。

References:

［1］ 昝灵，王顺华，张枝焕，张立生，邢辉.砂砾岩储层研究现状［J］.长江大学学报：自然科学版，2011，8（3）：63-66. ZAN Ling， WANG Shunhua， ZHANG Zhihuan， ZHANG Lisheng， XING Hui. Research status of sandy conglomerates reservoir［J］. Journal of Yangtze University: Natural Science Edition， 2011， 8（3）: 63-66.

［2］ 赵益忠，曲连忠，王幸尊.不同岩性地层水力压裂裂缝扩展规律的模拟实验［J］.中国石油大学学报：自然科学版，2007，31（3）： 63-66. ZHAO Yizhong， QU Lianzhong， WANG Xinzun. Simulation experiment on prolongation law of hydraulic fracture for different lithologic formations［J］. Journal of China University of Petroleum， 2007， 31（3）: 63-66.

［3］ 王昊.砂砾岩油藏砾石对压裂裂缝延伸的影响研究［D］.山东青岛：中国石油大学（华东）， 2011. WANG Hao. The study of hydraulic fracture propagation influenced by gravel in sand-gravel reservoir ［D］. China University of Petroleum （East China）， 2011.

［4］ 孟庆民，张士诚，郭先敏，陈新浩，张云.砂砾岩水力裂缝扩展规律初探［J］.石油天然气学报，2010，32（4）：119-123. MENG Qingmin， ZHANG Shicheng， GUO Xianmin，CHEN Xinhao， ZHANG Yun. A primary investigation on propagation mechanism for hydraulic fracture in glutenite formation［J］. Journal of Oil and Gas Technology， 2010，32（4）: 119-123.

［5］ 李根，唐春安，李连崇，梁正召.水压致裂过程的三维数值模拟研究［J］.岩土工程学报，2010，32（12）：1875-1881. LI Gen， TANG Chunan， LI Lianchong， LIAN Zhengzhao. Numerical simulation of 3D hydraulic fracturing process ［J］. Chinese Journal of Geotechnical Engineering， 2010，32（12）: 1875-1881.

［6］ 罗攀，李勇明，江有适，李源.砂砾岩水力裂缝延伸路径模拟研究［J］.油气地质与采收率，2013，12（5）： 103-106. LUO Pan， LI Yongming， JIANG Youshi， LI Yuan. Research on micro-morphology characteristics of hydraulic fractures for conglomerate reservoir.［J］. Petroleum Geology and Recovery Efficiency， 2013， 12 （5）: 103-106.

［7］ 李林地，张士诚，张劲，潘林华.缝洞型碳酸盐岩储层水力裂缝扩展机理［J］.石油学报，2009，30（4）：570-573. LI Lindi， ZHANG Shicheng， ZHANG Jin， PAN Linhua. Mechanism of hydraulic fracture propagation in fracturecavity carbonate reservoirs［J］. Acta Petrolei Sinica，2009， 30（4）: 570-573.

［8］ 李林地，张士诚，庚勐.煤层气藏水力裂缝扩展规律［J］.天然气工业，2010，30（2）：72-74. LI Lindi， ZHANG Shicheng， Geng Men. A study of the propagation law of hydraulic fractures in coalbed gas reservoirs［J］. Natural Gas Industry， 2010， 30（2）: 72-74.

［9］ Norbeck J H， Mcclure M W， Lo J W， Horne R N. An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation［J］. Computational Geosciences， 2016， 20（1）: 1-18.

［10］ ZHU Haiyan， ZHAO Xing， GUO Jianchun， JIN Xiaochun， AN Fengchen， WANG Yonghui， LAI Xiangdong. Coupled flow-stress-damage simulation of deviated-wellbore fracturing in hard-rock［J］. Journal of Natural Gas Science & Engineering， 2015， 26: 711-724.

［11］ 杜修力，金浏.考虑过渡区界面影响的混凝土宏观力学性质研究［J］.工程力学，2012，29（12）：72-79. DU Xiuli， JIN Liu. Research on the influence of interfacial transition zone on the macro-mechanical properties of concrete［J］. Engineering Mechanics，2012， 29（12）: 72-79.

［12］ 杜修力，金浏.混凝土材料宏观力学特性分析的细观单元等效化模型［J］.计算力学学报，2012，29（5）：654-661. DU Xiuli， JIN Liu. Meso-element equivalent model for macro-scopic mechanical properties analysis of concrete materials［J］. Journal of Computational Mechanics，2012， 29（5）: 654-661.

［13］ 温韬，唐辉明，刘佑荣，邹宗兴，王康，林成远.影响因子修正的新型岩石损伤统计本构模型［J］.中国矿业大学学报，2016，45（1）：141-149. WEN Tao， TANG Huiming， LIU Yourong， ZHOU Zongxin， WANG Kang， LIN Chenyuan. Newly modified damage statistical constitutive model of rock based on impact factor［J］. Journal of China university of mining， 2016， 45 （1）: 141-149.

［14］ 唐春安，刘红元，秦四清，杨志法.非均匀性对岩石介质中裂纹扩展模式的影响［J］.地球物理学报，2000，43（1）：116-121 . TANG Chunan， LIU Hongyuan， QIN Siqing， YANG Zhifa. Influence of heterogeneity on crack propagation modes in brittle rock［J］ . Journal of Geophysics， 2000 ，43（1）: 116-121.

［15］ BARTELMANN M， FABIS F， BERG D， KOZLIKIN E，LILOW R， VIERMANN C. Non-equilibrium statistical field theory for classical particles: Non-linear structure evolution with first-order interaction［J］. Angewandte Chemie International Edition in English， 2014， 3（7）: 519-519.

［16］ 于超云，唐世斌，唐春安.湿度场下矩形巷道围岩破坏规律的数值研究［J］.济南大学学报：自然科学版，2015，29（3）：204-209 . YU Chaoyun， TANG Shibin， TANG Chunan. Numerical analysis on the failure law of rectangular roadways under humidity conditions［J］. Journal of University of JiNan: Natural Science Edition， 2015， 29（3）:204-209.

［17］ BAI Yilong， WANG Haiying， XIA Mengfen， KE Fujiu. Statistical mesomechanics of solid linking coupled multiple space and time scales［J］. Applied Mechanics Reviews， 2005， 58（6）: 372-388.

［18］ NOIRIEL C， RENARD F， DOAN M L， GRATIER J P. Intense fracturing and fracture sealing induced by mineral growth in porous rocks［J］. Chemical Geology， 2010，269（3-4）: 197-209.

［19］ 李根，唐春安，李连崇.水岩耦合变形破坏过程及机理研究进展［J］.力学进展，2012，42（5）：593-619. LI Gen， TANG Chunan， LI Lianchong. Advances in rock deformation and failure processunder water-rock coupling［J］. Advances in Mechanics， 2012， 42（5）: 593-619.

（修改稿收到日期 2016-03-27）

〔编辑 朱 伟〕

Numerical simulation on hydraulic fracture propagation in glutenite reservoir based on
microscopic damage multiphase coupling

YU Donghe1， XU Kangtai2， CHE Hang1， ZHANG Dengwen1， LIU Guohua1， MA Xinfang2
1. Oil Production Engineering Research Institute of Huabei Oilfield Company， CNPC， Renqiu， Hebei 062552， China；2. College of Petroleum Engineering of China Uniνersity of Petroleum （Beijing）， Beijing 102249， China

The numerical simulation on hydraulic fracture propagation in low permeability glutenite oil reservoir is of great significance to the successful fracturing stimulation of this reservoir. By characterizing the gravel in the glutenite reservoir as a matrixinterface-gravel three modal textures， and assuming that the gravel distribution， geometry and reservoir property meet the random distribution， the mathematical characterization of glutenite reservoir was completed with the Moter-Carolo method. Taking into account the tri-phase coupling of reservoir seepage field， stress field and hydration propagation moisture field， and according to the principles of damage mechanics and fracture mechanics， the microscopic damage finite element method was used to establish a mathematical model

glutenite reservoir； hydraulic fracture propagation； multiphase coupling； microscopic damage； numerical simulation

马新仿（1972-），副教授，主要从事采油工程理论与技术研究和教学工作。通讯地址：（102249）北京市昌平区府学路18号。电话：010-89734593。E-mail： maxinfang@cup.edu.cnof hydraulic fracture propagation in glutenite reservoir. This mathematical model was used to simulate and analyze the hydraulic fracture propagation when the fractures encounter gravels under different principal stress differences， matrix-gravel interface strengths and gravel strengths. Finally， numerical simulation of dynamic propagation of hydraulic fractures in glutenite reservoir was completed. The study results show that some phenomena such as bypassing， passing through and arresting occur when hydraulic fracture encounters gravels. Especially， propagation through bypassing gravel predominates. The fracture apparently diverts， and pinnate secondary fractures occur. The fracture diverting level and penetration are related to the principal stress difference， gravel strength and interface strength. Essentially， the smaller the horizontal principal stress， the more apparent the diversion of hydraulic fracture when it encounters gravels. As the matrix-gravel interface strength increases， the hydraulic fracture apparently shortens and is difficult to divert. As the gravel strength increases， the diverting level of fractures increases.

TE357.1

A

1000 - 7393（ 2016 ） 03 - 0352- 07

10.13639/j.odpt.2016.03.015

YU Donghe， XU Kangtai， CHE Hang， ZHANG Dengwen， LIU Guohua， MA Xinfang. Numerical simulation on hydraulic fracture propagation in glutenite reservoir based on microscopic damage multiphase coupling［J］. Oil Drilling & Production Technology， 2016， 38（3）: 352-358.

中国石油天然气股份有限公司重大科技专项“华北油田上产稳产800万吨关键技术研究与应用”（编号：2014E-35-08-03）。

余东合（1966-），1987年毕业于西南石油学院石油地质专业，一直从事油田开发技术研究和管理工作，高级工程师。通讯地址：（062552）河北省任丘市华北油田公司采油工程研究院。电话：0317-2728443。E-mail： pjb_ydh@petrochina.com.cn

引用格式：余东合，徐康泰，车航，张登文，刘国华，马新仿.基于细观损伤多相耦合的砂砾岩水力压裂裂缝扩展数值模拟［J］.石油钻采工艺，2016，38（3）：352-358.