刚度差异对岩石破裂影响的颗粒流模拟研究

陈军涛，李 昊，褚衍玉，张 毅，唐道增，朱 君，李文昕，王开宇

(1. 山东科技大学 能源与矿业工程学院，山东 青岛 266590 ； 2.煤炭资源高效开采与洁净利用国家重点实验室，北京 100013； 3. 山东能源临沂矿业集团有限责任公司 邱集煤矿，山东 德州 251105； 4. 中国安全生产科学研究院，北京 100012； 5.山东科技大学 电气与自动化工程学院，山东 青岛 266590)

岩土工程的失稳破坏往往与岩体内部裂隙的张开、闭合、扩展及贯通密切相关，研究岩石破裂过程及其规律对岩土工程具有重大意义[1-2]。实际工程中很难采用探测手段观察到岩石的破裂过程[3]，颗粒流模拟则成为研究岩土类材料裂纹扩展破裂的重要手段，而宏观参数与细观参数的转换问题是模拟准确性的关键[4]。国内外学者对颗粒流模拟进行了大量的研究，李守巨等[5]提出一种基于巴西劈裂试验的混凝土细观本构模型参数估计方法；尹小涛等[6]研究了黏结内聚力和摩擦角对岩石各项宏观力学特征的影响；袁康等[7]利用离散元软件分析了加载过程中锚杆和结构面之间的相互作用机理以及锚杆内部轴向和剪切应力情况；邓树新等[8]利用Plackett-Burman试验设计并分析了细观参数对宏观响应的敏感性，建立了宏观力学指标和细观参数之间的线性关系；Yoon[9]采用中心合成设计(central composite design, CCD)标定了接触模型的细观参数，研究了单轴压缩强度和细观参数之间的关系；Coetzee等[10]通过组合的模拟剪切试验和压缩试验研究了颗粒间摩擦系数和黏结刚度；Stavrou等[11]采用离散元方法进行岩体的尺度分析，提出了块石强度与其体积和条件的关系。以上研究大多揭示了平行黏结模型细观参数(如黏结内聚力、摩擦角、黏结强度)与宏观参数之间的关系，而刚度参数与岩石的宏观破裂密不可分[12]，目前鲜有这方面的研究。

图1 线性平行黏结模型[9]Fig. 1 Linear parallel bond mode

1 试验方案

1.1 线性平行黏结模型

假设颗粒有效模量等于黏结有效模量，则法向刚度与切向刚度的比值等于线性平行黏结法向刚度与线性平行黏结切向刚度的比值。此时，模型弹性模量可表示为[16]：

E/Ec=a+bln(kn/ks)。

(1)

式中：E为弹性模量，GPa；Ec为有效模量，GPa；kn/ks为刚度比；a、b为常数。颗粒有效模量与法向刚度的关系如图2所示，表达式为[17]：

图2 线性平行黏结ball-ball接触[9]Fig. 2 Linear parallel bond contact of ball-ball

emod=knL/A，

(2)

(3)

(4)

(5)

式中：R(1)为接触两端某一个颗粒的半径；R(2)为接触两端另一个颗粒的半径；L为接触两端颗粒的中心距；r为较小颗粒的半径；A为接触的截面积；t为接触界面的厚度。

线性平行黏结有效模量与线性平行黏结法向刚度之间的关系[18]为：

(6)

综上分析，当模型刚度比确定时，颗粒、黏结有效模量可以控制颗粒和黏结的法向刚度和切向刚度，并且与宏观弹性模量相关。众所周知，岩石材料是典型的非均质体，而视为均质体的岩石试验模拟获得的结论与实际有一定差距。因此，本研究通过离散元模拟试验，研究细观参数颗粒有效模量、黏结有效模量对岩石裂纹扩展和破裂过程的影响。

图3 模型示意图Fig. 3 Model diagram

1.2 模拟方案设计

建立如图3所示的离散元模型，尺寸为50 mm×100 mm，线性平行黏结模型。模型细观参数颗粒有效模量、黏结有效模量为定值，孔隙率为0.1，颗粒半径0.75×10-3m，密度2 500 kg/m3，颗粒局部阻尼系数0.7。通过墙体建立颗粒生成范围，在范围内随机生成颗粒。

现定义刚度差异系数Ka，反映模型内部颗粒刚度和胶结体刚度的差异程度。

(7)

当Ka>1时，颗粒刚度大于胶结物刚度，Ka越大，颗粒与胶结体的刚度差别越大；当0

对模型赋予线性平行黏结后，查看接触Contactpb_state数值，确认接触模型赋值成功。模拟试验采用位移方式加载，轴向加载速率为0.2 mm/s。通过fish函数记录应力-应变变化过程、拉剪裂纹数目及裂纹扩展过程，峰值应力一半处视为弹性阶段，求得试样弹性模量，试验结束后保存模型并记录破裂形态。

2 颗粒有效模量对岩石宏观参数及破裂过程的影响特征

2.1 应力-应变曲线

图4为黏结有效模量一定时，颗粒有效模量分别为1、3、5、7、9和11 GPa时的单轴压缩全应力-应变曲线。可以看出，应力-应变曲线的前半部分基本为直线段，属于弹性变形阶段，峰前曲线明显偏离直线，偏离直线的起点视为屈服点，屈服点到峰值之间为非稳定破裂发展阶段，峰后试样表现出明显的脆性。模拟曲线并未表征试验具有明显的裂隙压密阶段，与室内试验得到的裂隙压密、弹性变形、微弹性裂隙稳定发展、非稳定破裂发展等阶段有差异，这主要是因为软件很难准确模拟岩石复杂的原生孔隙和裂隙网络，而默认为致密模型，故在模拟结果中没有明显的裂隙孔隙压密的现象。

图4 不同颗粒有效模量的应力-应变曲线(单轴压缩)Fig. 4 Stress-strain curves of effective modulus of different particles (uniaxial compression)

从图4中还可以看出：岩石试样破裂前，弹性模量和单轴抗压强度随颗粒有效模量的增大而增大；轴向应变值随颗粒有效模量的增大而减小，呈负相关关系。当颗粒有效模量较大时，模型更容易出现脆性破坏；当颗粒有效模量较小时，模型更容易出现塑性破坏。

PFC2D模拟试验中，颗粒刚度通过法向刚度和切向刚度控制，颗粒刚度越大，颗粒越难以发生弹性变形。由式(2)可知，颗粒有效模量与法向刚度正相关，在黏结有效模量和刚度比确定的情况下，颗粒有效模量越大，法向刚度和切向刚度越大，同样应力条件下模型变形量越小。

黏结有效模量一定时，颗粒有效模量为1、3、5、7、9和11 GPa时，试样拉伸应力峰值、弹性模量和应变峰值几乎没有影响，仅颗粒有效模量为11 GPa的试样在应力峰值后发生了小幅波动。由式(1)可知，颗粒有效模量与颗粒刚度成正相关，拉伸应力主要作用在颗粒与颗粒之间的接触上。因此，在细观参数黏结有效模量确定的情况下，拉伸条件下试样各宏观参数几乎不受颗粒有效模量变化的影响。

2.2 颗粒有效模量对岩石宏观参数的影响分析

根据2.1的研究结果，改变颗粒有效模量时，对岩石压缩过程中弹性模量、峰值应变、单轴抗压强度等宏观参数的变化进行分析。单轴抗压强度和应变随颗粒有效模量的增大呈非线性变化趋势，当颗粒有效模量不超过7 GPa时对模型的单轴抗压强度影响较大，大于7 GPa时模型的单轴抗压强度趋于稳定。当颗粒有效模量不超过5 GPa时对模型的应变影响较大，大于5 GPa时模型的应变趋于稳定。反映颗粒有效模量增大到一定数值(本试验为7 GPa)后对模型单轴抗压强度、应变影响越来越小。

图5～7分别为压缩试验过程中，不同颗粒有效模量的岩石弹性模量、峰值应变、单轴抗压强度的变化曲线，可以看出随着颗粒有效模量的增大，弹性模量和单轴抗压强度均表现出上升趋势，峰值应变则表现出下降趋势。对于拟合结果的R2，弹性模量(0.99)>峰值应变(0.87)>单轴抗压强度(0.84)，表明宏观弹性模量与细观颗粒有效模量呈较好的线性关系；对于拟合结果的斜率，单轴抗压强度(5.20)>弹性模量(0.71)>峰值应变(0.26)，说明颗粒有效模量发生改变时，岩石的单轴抗压强度变化最大。

2.3 刚度差异系数对岩石破裂的影响分析

图8为颗粒有效模量影响下不同刚度差异系数值时的岩石的压缩破裂形态，图9为颗粒有效模量影响下岩石剪切裂纹数量、拉伸裂纹数量、裂纹总数随刚度差异系数的变化曲线，图10为颗粒有效模量影响下拉伸、剪切裂纹占裂纹总数的比例随刚度差异系数的变化曲线。

从图8中可以看出，总体来说，刚度差异系数越大，岩石破裂产生的裂纹越集中；刚度差异系数越小，岩石破裂产生的裂纹越分散，说明刚度差异系数的大小影响着岩石破裂过程中的裂隙分布形态。从图9中可以看出，当刚度差异系数较小(≤7)时，拉伸裂纹数量和裂纹总数变化不稳定，当刚度差异系数较大(>7)时，拉伸裂纹数量和裂纹总数变化趋于稳定，而剪切裂纹数量较少，在31～63范围内波动。图10表明，随刚度差异系数增大，拉伸裂纹和剪切裂纹占比仅在刚度差异系数为3时有较大波动，刚度差异系数为其他值时拉伸裂纹总数保持在80%左右，剪切裂纹总数保持在20%左右，说明刚度差异系数增大对岩石的剪切拉伸两类裂纹占比几乎没有影响。

图5 弹性模量Fig. 5 Elastic modulus

图6 峰值应变Fig. 6 Peak strain

图7 单轴抗压强度Fig. 7 Uniaxial compressive strength

图8 压缩破裂形态Fig. 8 Compression fracture morphology

图9 裂纹数量Fig. 9 Number of cracks

图10 拉剪裂纹占比Fig. 10 Proportion of tensile and shear cracks

3 黏结有效模量对岩石宏观参数及其破裂过程的影响特征

3.1 应力-应变曲线

图11为颗粒有效模量一定时，黏结有效模量分别为1、3、5、7、9和11 GPa时的单轴压缩全应力-应变曲线。从图11中可以看出：岩石破裂前，曲线斜率随黏结有效模量的增大而增大，即弹性模量随黏结有效模量的增大而增大；单轴抗压强度和轴向应变随黏结有效模量的增大而减小。当黏结有效模量较大时，模型更容易出现脆性破坏，反之更容易出现塑性破坏，但是影响程度较小。

由式(1)可知，黏结有效模量与模型的黏结法向刚度成正相关关系，在颗粒有效模量确定的情况下，增大黏结有效模量相当于增大了接触弹性体的刚度，同样应力条件下弹性体的刚度越小模型变形量越小。

图12为黏结有效模量分别为1、3、5、7、9、11 GPa时的拉伸应力-应变曲线，可以看出随着黏结有效模量增大，岩石的弹性模量增大，相同应力水平下应变量减小，峰值应力几乎保持不变。

图11 不同黏结有效模量单轴压缩应力-应变曲线Fig. 11 Stress-strain curves of uniaxial compression with different bonding effective modulus

图12 不同黏结有效模量拉伸应力-应变曲线Fig. 12 Tensile stress-strain curves of different bonding effective modulus

3.2 黏结有效模量对岩石宏观力学参数的影响分析

图13～15为颗粒有效模量一定时，单轴压缩试验过程中试样压弹性模量、压缩峰值应变和单轴抗压强度随黏结有效模量的变化趋势。图16～18为颗粒有效模量一定时，单轴拉伸试验过程中试样拉弹性模量、拉伸峰值应变、抗拉强度随黏结有效模量的变化趋势。

图13 压弹性模量Fig. 13 Elastic modulus under compression

图14 压缩峰值应变Fig. 14 Compression peak strain

图15 单轴抗压强度Fig. 15 Uniaxial compressive strength

图16 拉弹性模量Fig. 16 Tensile modulus of elasticity

图17 拉伸峰值应变Fig. 17 Tensile peak strain

图18 抗拉强度Fig. 18 Tensile strength

单轴抗压强度、拉应变、压应变随黏结有效模量的增大呈非线性变化趋势，黏结有效模量变化对模型抗拉强度几乎没有影响。当黏结有效模量不超过5 GPa时对模型的拉、压应变影响较大，当黏结有效模量大于5 GPa时模型的拉、压应变趋于稳定；当黏结有效模量不超过3 GPa时对模型的单轴抗压强度影响较大，当黏结有效模量大于3 GPa时模型的单轴抗压强度变趋于稳定。反映黏结有效模量增大到一定数值(本试验为5 GPa)后对模型单轴抗压强度、拉应变、压应变影响越来越小。

通过观察拟合结果的R2，拉弹性模量(1)>压弹性模量(0.995)>单轴抗压强度(0.75)>单轴抗拉强度(0.73)>压缩峰值应变(0.72)>拉伸峰值应变(0.66)，发现宏观参数拉、压弹性模量与细观参数黏结有效模量表现出良好的线性关系。从拟合结果斜率来看，拉弹性模量(1.61)>压弹性模量(1.15)>单轴抗压强度(0.97)>压缩峰值应变(0.33)>拉伸峰值应变(0.22)>单轴抗拉强度(0.04)，说明黏结有效模量增大对拉、压弹性模量的影响最大。

3.3 刚度差异系数对岩石破裂的影响分析

图19为黏结有效模量变化影响下不同刚度差异系数值时试件的压缩破裂形态，图20为黏结有效模量变化影响下剪切裂纹数量、拉伸裂纹数量和裂纹总数随刚度差异系数变化曲线，图21为黏结有效模量变化影响下拉伸、剪切裂纹数量占比随刚度差异系数变化曲线。结果表明，刚度差异系数较小时(≤1/5)，产生的裂纹数量较多，且裂纹数量变化较大；刚度差异系数较大时(>1/5)裂纹数目较稳定。随着刚度差异系数的减小，拉、剪裂纹的占比上下波动，剪切裂纹保持在16%～31%，拉伸裂纹占比保持在69%～84%，拉伸裂纹是剪切裂纹的2～4倍，说明刚度差异系数减小时对岩石的裂纹演化及破裂有一定影响，但总体影响不大。

图19 压缩破裂形态Fig. 19 Compression fracture morphology

图20 裂纹数量变化Fig. 20 Change of crack number

图21 拉剪裂纹占比Fig. 21 Proportion of tensile and shear cracks

4 实例验证

以文献[19]模拟结果为例，参考上述本研究结果，对单轴压缩试验进行参数优化研究。图22(a)为文献[19]的室内试验结果，图22(b)为文献[19]的模拟试验结果(Ka=1)，可以发现试件破坏模式以剪切为主，形成了与试件轴向呈一定夹角的宏观剪切裂纹，而模拟试验结果中裂隙分布散乱，岩石裂隙没有相互贯通形成宏观的剪切裂隙，与室内试验获得的结果相差较大。

图22 岩石破裂特征Fig. 22 Rock fracture characteristics

由前面研究可知，当颗粒有效模量或黏结有效模量增大时，试样越容易出现剪切破坏，且颗粒有效模量影响更大；刚度差异系数值增大时，裂纹数目减小。因此，保持其他参数不变，调整颗粒有效模量，增大刚度差异系数值，使得岩石裂纹集中在实际破裂区域，最终出现明显的宏观剪切裂隙，破坏形态与室内岩石破裂吻合性较好，调整前后的模拟试验结果见图22，细观力学参数见表1，验证了前面研究结论的准确性。

表1 调整前后的细观力学参数Tab. 1 Meso mechanical parameters before and after adjustment

5 结论

1) 模型黏结有效模量恒定时，颗粒有效模量越大，宏观参数压弹性模量和单轴抗压强度越大，压缩应变反而减小；颗粒有效模量对拉伸条件下的宏观参数几乎没有影响；当颗粒有效模量较大时，模型易出现脆性破坏，反之易出现塑性破坏；颗粒有效模量与压弹性模量具有较好的线性关系。

2) 模型颗粒有效模量恒定时，黏结有效模量越大，宏观参数压弹性模量和拉弹性模量越大，单轴抗压抗拉强度、压缩和拉伸应变越小；当黏结有效模量较大时，模型易出现脆性破坏，反之易出现塑性破坏；黏结有效模量与压、拉弹性模量都具有良好的线性关系。模拟试验时，应优先标定黏结有效模量，以提高模拟结果的准确性。

3) 刚度差异系数反映了模型内部颗粒和胶结体刚度的差异程度。刚度差异系数越大，模型越难产生裂隙，稳定性好；刚度差异系数越小，模型越易产生裂隙，稳定性差。

4) 以室内岩石试验为例，对细观参数进行了优化，参数优化后的试验结果与室内岩石破裂形态吻合较好，验证了本研究结论的准确性。