伯格斯模型参数调试与岩石蠕变特性颗粒流分析

杨振伟，金爱兵，周 喻，严 琼，王 凯，高玉娟

（1.北京科技大学 土木与环境工程学院，北京 100083；2.北京大学 前沿交叉学科研究院，北京 100083）

1 引 言

岩石的力学特性不仅表现出弹性和塑性，还具有与时间相关的性质——流变特性。随着地下工程规模的不断扩大以及深度的逐渐增加，岩石蠕变效应越来越明显，对于蠕变特性的研究也越来越受到重视。Burgers模型是常用的流变模型之一，可以反映岩石的黏弹性，很多学者在Burgers模型的基础上做了改进研究[1－2]。

室内蠕变试验主要集中在单轴蠕变试验[3]和不同围压下的三轴蠕变试验[4－11]。熊良宵等[12]对Burgers模型进行了修正，并将之用于蠕变试验的拟合分析。但室内试验往往具有操作费时费力、载荷加载时间长、试验费用高等缺点。

近年来，许多学者通过数值模拟的方法进行岩石流变特性的研究[13－14]。Sharifzadeh等[15]在室内三轴蠕变试验的基础上，研究了岩体的蠕变特性对于隧道变形的影响，并在有限差分程序FLAC中用Burgers模型进行模拟；张志沛等[16]利用三维有限差分程序(FLAC3D)中的Burgers模型对不同形状的试样进行了蠕变特性模拟分析。

上述模拟和计算均建立在连续介质理论的基础上，而岩石由于长期的地质作用和复杂自然环境的影响，存在着明显的不连续区域和不同程度的缺陷。此时，采用连续介质理论的数值分析方法对岩体进行分析存在一些不足。由Cundall[17]提出的颗粒流法及PFC程序是求解非连续介质力学问题的一个重要数值分析方法。该程序将介质离散为大量的小颗粒单元，并通过颗粒之间的相互作用去描述、研究非连续岩石材料的力学特性。Kang等[18]基于二维颗粒流程序(PFC2D)，采用法向Hertz-Mindlin模型、切向Burgers模型进行了双轴蠕变试验的研究；王涛等[19]在PFC2D中开发出广义Kelvin本构模型，并用于工程计算。结果表明，与线弹性接触模型相比，采用广义Kelvin接触模型得到的结果与由Hoek-Brown强度准则得到的结果更为接近。

本文在对三维颗粒流程序PFC3D中的Burgers本构模型进行深入研究的基础上，分析了Burgers模型中各参数对于瞬时强度特性和流变特性的影响，总结得出参数调试的方法，并结合国内已有的室内试验数据，进行了岩石试样蠕变试验的模拟和研究。

2 Burgers流变模型

2.1 基于连续介质理论的Burgers模型

Burgers模型是由一个马克斯伟尔体和一个开尔文体串联而成的流变模型，其结构如图1所示。根据模型组成特点，基于连续介质理论可以分别求出模型的本构方程、蠕变方程、卸载方程和松弛方程[20]，并绘制蠕变和卸载曲线，如图2所示。

图1 Burgers流变模型Fig.1 Burgers rheology model

图2 Burgers模型流变特性曲线Fig.2 Rheological curve of Burgers model

Burgers模型本构方程为

式中：σ为应力；ε为应变；Em为马克斯伟尔体弹性系数；ηm为马克斯伟尔体黏性系数；Ek为开尔文体弹性系数；ηk为开尔文体黏性系数。

蠕变方程为

式中：σ0为恒定应力；t为加载时间。

卸载方程为

松弛方程为

2.2 基于离散单元法的Burgers模型

连续介质理论中的Burgers模型将试样视为连续体，模型作用的对象是整个连续的试样。而在离散单元法中，模型作用于两个存在接触的实体（颗粒与颗粒或颗粒与墙体），数值积分方案如下[21]：

（1）对于模型中的开尔文体，可得

式中：uk为开尔文体相对位移；F为接触力，±分别表示法向和切向。

利用有限差分法，取 uk和F 的平均值，可得

式中：Δt为时间步长；t 表示当前时步计算结果；t+1表示下一时步计算结果。

对上式进行整理，可以得到 uk的计算公式为

式中：A、B为系数，可由下式求得

（2）对于模型中的马克斯伟尔体，可得

式中：um为马克斯伟尔体相对位移。

利用有限差分法，取F 的平均值，可以得到：

对上式进行整理，可以得到 um的计算公式为

（3）对于全部Burgers模型，有

式中：u为Burgers模型相对位移。

颗粒之间下一时步的接触力Ft+1可以由ut+1、ut、、Ft计算得到：

式中：C、D为系数，可有下式求得

3 PFC3D 中Burgers模型参数调试

在连续介质理论中，Burgers模型的有4个参数，而在离散单元法中，Burgers模型所需的参数有9个，如表1所示。

表1 Burgers模型参数Table 1 Parameters of Burgers model

为了明确各个参数对于瞬时强度特性和流变特性的影响，本文采用控制变量法对各参数进行分类研究。本文通过Fishtank生成直径为50 mm、高为100 mm的圆柱体，最小颗粒半径Rmin=1.38 mm，最大/最小粒径比为1.66。为了方便调试，本文中将各参数法向和切向数值设为相同，即需要调试的参数为5个，分别是：Em、ηm、Ek、ηk和f。

针对上述参数，本文进行5组对照试验，每组试验3次，同一组内的试验之间只变化一个参数，其余参数均取为相同，由此得到各参数对瞬时强度特性和流变特性的影响，每组控制试验所取参数和瞬时强度试验结果如表2所示。

表2 对照试验参数取值和瞬时强度试验结果Table 2 Parameters in controlled experiments and results of instantaneous strength tests

3.1 Burgers模型参数对瞬时强度特性的影响

利用Burgers模型进行蠕变试验或工程计算之前，需要对模型参数进行调试，使得弹性模量、泊松比、单轴抗压强度、黏聚力和内摩擦角等瞬时强度特性与室内试验结果相符合。此时，数值试样可体现真实试样的性质，并用于蠕变试验或工程计算。

按照表2中所列参数取值，进行了5组瞬时强度特性对照试验，并将各组试验的弹性模量、泊松比和单轴抗压强度的变化在表2中列出。根据试验结果，将各参数影响和调试方法总结如下：

（1）由第1组对照试验可知，随着弹性系数Em的增加，弹性模量、泊松比和单轴抗压强度均有所增加，说明瞬时强度特性与弹性系数Em呈正相关。弹性元件的应力-应变与时间无关，可产生瞬时应力和瞬时应变，马克斯伟尔体中的弹性元件与其他部分串联，可独立发挥作用，对弹性模量影响较大。

（2）由第2组和第4组对照试验可知，即使马克斯伟尔体和开尔文体的黏性系数(ηk和ηm)发生较大变化，瞬时强度特性未发生明显变化。说明黏性系数(ηk和ηm)对于瞬时强度特性的影响不大，在匹配瞬时强度特性时无需调整。黏性元件的应力-应变与时间相关，在瞬时强度试验中，由于作用时间较短，应变量较小，对接触力的影响不大，所以黏性元件几乎不发挥作用。

（3）由第3组对照试验可知，当开尔文体的弹性系数Ek发生变化时，弹性模量、泊松比和单轴抗压强度均未发生明显变化。说明开尔文体的弹性系数Ek对于瞬时强度特性的影响不大，在匹配瞬时强度特性时无需调整。开尔文体中的弹性元件由于与黏性体并联，应变受到黏性元件限制，因此，在较短时间内无法像马克斯伟尔弹性体那样发挥作用。

（4）由第5组对照试验可知，随着摩擦因数的增大，弹性模量逐渐增加，泊松比逐渐减小，单轴抗压强度逐渐增加。摩擦因数f 与弹性模量和单轴抗压强度呈正相关，与泊松比呈负相关。摩擦因数通过控制接触实体之间的相对滑动，实现对瞬时强度特性的影响。

综上所述，Burgers模型中参数对瞬时特性有影响的主要是马克斯伟尔体弹性系数Em和摩擦因数f，其余参数均未产生明显影响，在参数调试过程中，只需调整参数Em和f 即可。Burgers模型中各参数对瞬时强度特性的影响关系如表3所示。

表3 Burgers模型参数对瞬时强度特性的影响Table 3 Impact of Burgers model parameters on instantaneous strength

3.2 Burgers模型参数对流变特性的影响

按照表2中所列参数取值，进行了5组流变特性对照试验，分析如下：

（1）马克斯伟尔体弹性系数Em

第1组对照试验主要研究马克斯伟尔体弹性系数Em的影响。Em取不同值时，得到蠕变曲线和卸载曲线如图3所示。由图可知，Em主要影响Burgers模型的瞬时应变量和瞬时回复量。

图3 马克斯伟尔体弹性系数Em的影响Fig.3 Axial strain with time for impact of Maxwell’s elastic coefficient Em

在应力施加的瞬间，试样产生瞬时应变，Em越大，瞬时应变的数值就越小。在卸载阶段，瞬时回复量与加载时瞬时应变量相同。

（2）马克斯伟尔体黏性系数ηm

第2组对照试验主要研究马克斯伟尔体黏性系数ηm的影响。ηm取不同值时，得到蠕变曲线和卸载曲线如图4所示。由图可知，ηm主要影响Burgers模型的起始蠕变率、起始蠕变量、稳定蠕变率和残余应变。

在应力施加的瞬间，3次试验的瞬时应变相同，但随着ηm的增大，起始蠕变率和稳定蠕变率均减小，蠕变曲线逐渐趋于平缓。在卸载过程中，ηm对于瞬时回复无影响，但对残余应变的影响比较显著，随着ηm的增大，残余应变逐渐减小。

图4 马克斯伟尔体黏性系数ηm的影响Fig.4 Axial strain with time for impact of Maxwell’s viscosity coefficient ηm

（3）开尔文体弹性系数Ek

第3组对照试验主要研究开尔文体弹性系数Ek的影响。Ek取不同值时，得到蠕变曲线和卸载曲线如图5所示。由图可知，Ek主要影响Burgers模型的起始蠕变量和弹性后效回复量。

图5 开尔文体弹性系数Ek的影响Fig.5 Axial strain with time for impact of Kelvin’s elastic coefficient Ek

在应力加载瞬间，瞬时应变均相同，起始蠕变率和稳定蠕变率也几乎相同，但由起始蠕变过渡到稳定蠕变的时间不同，从而导致起始蠕变量不同。在卸载过程中，瞬时弹性回复几乎相同，弹性后效回复量随着Ek的增大而减小，最终残余应变也几乎相同，3次试验曲线逐渐趋于相同。

（4）开尔文体黏性系数ηk

第4组对照试验研究开尔文体黏性系数ηk的影响。ηk取不同值时，得到蠕变曲线和卸载曲线如图6所示。由图可知，ηk主要影响Burgers模型起始蠕变率和弹性后效回复率。

图6 开尔文体黏性系数ηk的影响Fig.6 Axial strain with time for impact of Kelvin’s viscosity coefficient ηk

在蠕变试验过程中，瞬时应变量和稳定蠕变率均相同，起始蠕变率却受到开尔文体黏性系数ηk的影响，ηk的取值越大，起始蠕变率越慢，最终达到稳定蠕变的时间越长。在卸载阶段，弹性后效回复率变化很大，ηk越大，弹性后效蠕变速率越慢，最终达到残余应变的时间越长。

（5）摩擦因数f

第5组对照试验主要研究摩擦因数f 的影响。f取不同值时，得到蠕变曲线和卸载曲线如图7所示。由图可知，f 主要影响Burgers模型的瞬时应变量和瞬时回复量。

图7 摩擦因数f 的影响Fig.7 Axial strain with time for impact of friction factor f

在加载过程中，瞬时应变量随着摩擦因数的增大而减小，在蠕变过程中，3次试验的曲线平行，说明起始蠕变段和稳定蠕变阶段不受f 影响。卸载阶段，瞬时回复量与相应的瞬时应变量相同。

综合以上5组对照试验，可以得到Burgers模型各参数对于试样流变特性的影响，见表4。

表4 Burgers模型参数对流变特性的影响Table 4 Impact of Burgers model parameters on rheological properties

4 Burgers模型应用实例

4.1 工程背景

贵州开磷集团下属的多个矿山开采深度相继进入600 m以下。何涛[22]以深部红页岩为研究对象，进行了瞬时强度特性试验和室内蠕变试验，分析了试样所具有的黏弹性特征，并选用Burgers模型进行拟合分析，所得相关性系数较高，从而验证了Burgers模型用于红页岩等软岩流变特性模拟的合理性。

不足之处在于，每次试验选取不同参数分别拟合，得到的单个试验拟合效果较好，但与实际情况不符，同一岩石试样在不同试验中所选取参数应相同。

本文应用三维颗粒流程序(PFC3D)进行单轴和三轴蠕变试验，在PFC3D中Burgers模型作为用户自定义模型被调用，可以用来模拟软岩的流变特性。对同一岩石试样，每次试验均选取相同细观参数进行计算。

4.2 颗粒流试样生成

在PFC3D中，利用墙体的概念来模拟介质的边界。试样生成主要包括以下4个阶段：

（1）生成指定数量和大小的颗粒；

（2）通过循环消除试样内部的非均匀应力；

（3）消除悬浮颗粒。本文中，将与周围颗粒接触数量小于3个的颗粒定义为悬浮颗粒。这些颗粒的存在导致试样内部产生“空区”，并影响模型的受力状态，造成计算结果的失真，因此，需要消除；

（4）赋予黏结属性，最终生成试样，本文中颗粒之间赋予平行黏结。

4.3 瞬时强度试验

试样生成后，调用Burgers模型，代替默认的线性模型，并赋予参数，通过调用Fishcall函数，在计算过程中，如果有新的接触产生，则将其赋予Burgers模型及参数值。根据本文第3节中的参数调试方法，经调试后，所采用的颗粒细观力学参数如表5所示。颗粒流试样和单轴压缩试验所得应力-应变曲线如图8所示。

表5 数值试样细观力学参数Table 5 Mesomechanical parameters of numerical samples

图8 颗粒流试样和应力-应变曲线Fig.8 Particle flow specimen and stress-strain curve

在进行三轴压缩试验时，分别将围压设置为1、2、5、10 MPa，压缩至试件破坏。根据三轴试验结果，绘制摩尔圆和强度包络线，最终确定试样的内摩擦角和黏聚力。

最终所得数值试样的瞬时强度参数与室内试验数据对比如表6所示。通过与室内试验数据对比可以看出，数值试样与真实试样试验结果较为吻合，可以用来进行蠕变试验。

表6 红页岩室内试验强度特性Table 6 Instantaneous properties of red shale in laboratory

4.4 单轴蠕变试验

（1）分级增量加载试验

通过删去试样侧墙，实现单轴压缩应力状态，采用伺服机制将轴向压力分别设置为7.5、10.0、12.5、17.5 MPa。分级增量加载试验轴向应变呈台阶状，蠕变过程中的轴向应变如图9所示。由图可知，各级应力水平下蠕变速率和蠕变量随着应力增大而增大。

将数值试验与室内蠕变试验结果进行对比分析，数值试样在压缩试验的前期，应力-应变表现出线性关系，如图8所示。但实际岩石试样往往存在一定的微裂隙，应力-应变关系是非线性的，如图10所示。应力-应变曲线可划分为5个阶段[20]：①微裂隙压密阶段(OA 段)；②弹性变形阶段(AB 段)；③裂隙发生和拓展阶段(BC 段)；④裂隙不稳定发展直到破裂阶段(CD 段)；⑤破裂后阶段(DE 段)。

图9 分级增量加载试验轴向应变Fig.9 Axial strains with time at staged increment loading tests

图10 室内单轴压缩试验应力-应变曲线Fig.10 Stress-strain curve of uniaxial compression test in laboratory

针对微裂隙压密阶段(OA 段)，赵东宁等[23]提出，变形可以分为由微裂隙被压密产生的变形和试样的压缩弹性变形两部分，其中前者占主导，并给出轴向应变ε3计算公式为

式中：ε1为微裂隙被压密产生的变形；ε2为试样本身受压产生的弹性压缩变形，可由下式计算：

式中：σ1为试样被压密时的压密强度；E为试样弹性模量。

在微裂隙压密阶段，可通过消除微裂隙压密应变ε1的方法，得到等效线弹性岩石试样的应力-应变曲线。将岩石试样的线弹性阶段AB 反向延长与坐标轴交于O′点，曲线O′AB 可认为是对压密后等效线弹性试样的压缩。将数值试验结果与消除裂隙压密应变ε1之后的结果进行拟合，如图11所示。由图可知，数值试验结果与等效线弹性阶段拟合效果良好。（2）分级增量循环加、卸载试验

图11 数值试验结果与室内试验数据对比图Fig.11 Axial strains with time from results of numerical tests and creep tests in laboratory

在分级增量循环加、卸载试验中，轴向压力分别设置为7.5、10.0、12.5、17.5 MPa。每级应力水平加载一定时间后，将轴向压力卸载至0，持续一段时间，再进行下一级加载，试验中的轴向应变如图12所示。由图可知，轴向应变呈锯齿状，且蠕变量随着应力水平的增大而增大。从起始蠕变过渡到稳定蠕变的时间逐渐增长。

图12 分级增量循环加卸载轴向应变Fig.12 Axial strains with time at staged increment cyclic loading and unloading experiments

4.5 三轴压缩蠕变试验

采用自编伺服机制，控制轴向和侧向应力水平，进行三轴蠕变试验，为深入研究三轴蠕变试验特性，分别进行了固定围压和固定轴压两组试验，加载方式为分级增量加载。

（1）固定围压

将围压设置为5 MPa，轴压取不同应力值时，轴向应变-径向应变曲线如图13所示。

图13 围压5 MPa轴向应变和径向应变Fig.13 Axial and lateral strains with time at confining stress 5 MPa

由图可知，三轴蠕变试验中轴向应变也呈台阶状，但各应力水平应变值均小于单轴压缩蠕变试验值。径向应变也有相同蠕变趋势，每个应力水平加载瞬间产生瞬时应变，试样径向扩大，然后在蠕变作用下逐渐收缩。

（2）固定轴压

将轴压固定为20 MPa，围压取不同值时，轴向和径向蠕变曲线如图14所示。从图中可以看出，每级应力加载瞬间，试样轴向伸长，径向收缩。而随着时间的推移，轴向逐渐收缩，径向却一直收缩。在每一级应力水平加载瞬间轴向伸长，径向收缩。围压达到20 MPa时，径向应变已经超过轴向应变。

图14 固定轴压20 MPa的轴向和径向应变Fig.14 Axial and lateral strains with time at axial stress 20 MPa

5 结 论

（1）数值试样的瞬时强度主要影响因素为马克斯伟尔体弹性系数Em和摩擦因数f：弹性模量和单轴抗压强度与Em和f 均呈正相关；泊松比与Em呈正相关，与f 呈负相关；其余参数未产生明显影响。

（2）岩体流变特性受Burgers模型的弹性系数(Em和Ek)、黏性系数(ηk和ηk)及摩擦因素f 共同影响，岩体流变特性与各参数均呈负相关。

（3）模拟结果与室内单轴蠕变试验结果对比表明，PFC3D中通过Fishtank生成的Burgers模型对单轴蠕变试验中的等效线弹性段模拟效果良好。

（4）通过固定轴压和固定围压两种方式，开展了三轴压缩蠕变试验，得到的蠕变曲线与理论曲线具有相同趋势，验证了Burgers模型在三轴蠕变试验模拟中的适用性。

[1]康永刚，张秀娥.岩石蠕变的非定常分数伯格斯模型[J].岩土力学，2011，32(11):3237－3241.KANG Yong-gang，ZHANG Xiu-e.Nonstationary parameter fractional Burgers model of rock creep[J].Rock and Soil Mechanics，2011，32(11):3237－3241.

[2]何利军，孔令伟，吴文军，等.采用分数阶导数描述软黏土蠕变的模型[J].岩土力学，2011，32(增刊2):239－243.HE Li-jun，KONG Ling-wei，WU Wen-jun，et al.A description of creep model for soft soil with fractional derivative[J].Rock and Soil Mechanics，2011，32(Supp.2):239－243.

[3]LI Yong-sheng，XIA Cai-chu.Time-dependent tests on intact rocks in uniaxial compression[J].International Journal of Rock Mechanics and Mining Sciences，2000，37(3):467－475.

[4]ZHANG HUABIN，WANG ZHIYIN，ZHENG YALI，et al.Study on triaxial creep experiment and constitutive relation of different rock salts[J].Safety Science，2012，50(4):801－805.

[5]徐进，张家生，赵同顺，等.软弱路基土体三轴蠕变试验及蠕变模型研究[J].中南大学学报(自然科学版)，2011，42(10):3136－3142.XU Jin，ZHANG Jia-sheng，ZHAO Tong-shun，et al.Creep triaxial tests and constitutive model of embankment foundation soil[J].Journal of Central South University(Science and Technology)，2011，42(10):3136－3142.

[6]ZHANG Z L，XU W Y，WANG W，et al.Triaxial creep tests of rock from the compressive zone of dam foundation in Xiangjiaba hydropower station[J].International Journal of Rock Mechanics and Mining Sciences，2012，50(1):133－139.

[7]ZHANG YU，XU WEIYA，GU JINJIAN，et al.Triaxial creep tests of weak sandstone from fracture zone of high dam foundation[J].Journal of Central South University，2013，20(9):2528－2536.

[8]YANG SHENGQI，JIANG YUZHOU.Triaxial mechanical creep behavior of sandstone[J].Mining Science and Technology，2010，20(3):339－349.

[9]ZHAO YANLIN，CAO PING，WANG WEIJUN，et al.Viscoelasto-plastic rheological experiment under circular increment step load and unload and nonlinear creep model of soft rocks[J].Journal of Central South University of Technology，2009，16(3):484－491.

[10]MARANINI E，BRIGNOLI M.Creep behaviour of a weak rock:Experimental characterization[J].International Journal of Rock Mechanics and Mining Sciences，1999，36(1):127－138.

[11]李亚丽，于怀昌，刘汉东.三轴压缩下粉砂质泥岩蠕变本构模型研究[J].岩土力学，2012，33(7):2035－2040.LI Ya-li，YU Huai-chang，LIU Han-dong.Study of creep constitutive model of silty mudstone under triaxial compression[J].Rock and Soil Mechanics，2012，33(7):2035－2040.

[12]熊良宵，杨林德，张尧.岩石的非定常 Burgers模型[J].中南大学学报(自然科学版)，2010，41(2):679－684.XIONG Liang-xiao，YANG Lin-de，ZHANG Yao.Non-stationary Burgers model for rock[J].Journal of Central South University (Science and Technology)，2010，41(2):679－684.

[13]徐平，李云鹏，丁秀丽，等.FLAC3D黏弹性模型的二次开发及其应用[J].长江科学院院报，2004，21(2):10－13.XU Ping，LI Yun-peng，DING Xiu-li，et al.Secondary development and application of visco-elastic constitutive model in FLAC3Dsoftware[J].Journal of Yangtze River Scientific Research Institute，2004，21(2):10－13.

[14]褚卫江，徐卫亚，杨圣奇，等.基于FLAC3D岩石黏弹塑性流变模型的二次开发研究[J].岩土力学，2006，27(11):2005－2010.CHU Wei-jiang，XU Wei-ya，YANG Sheng-qi，et al.Secondary development of a viscoelasto-plastic rheological constitutive model of rock based on FLAC3D[J].Rock and Soil Mechanics，2006，27(11):2005－2010.

[15]SHARIFZADEH M，TARIFARD A，MORIDI M A.Time-dependent behavior of tunnel lining in weak rock mass based on displacement back analysis method[J].Tunnelling and Underground Space Technology，2013，38:348－356.

[16]张志沛，王芝银，彭惠.陕南泥岩三轴压缩蠕变试验及其数值模拟研究[J].水文地质工程地质，2011，38(1):53－58.ZHANG Zhi-pei，WANG Zhi-yin，PENG Hui.Study on triaxial compression creep test and numerical simulation about mud rock in southern Shaanxi[J].Hydrogeology &Engineering Geology，2011，38(1):53－58.

[17]CUNDALL P A，STRACK O D L.A discrete numerical model for granular assemblies[J].Geotechnique，1979，29(1):47－65.

[18]KANG D H，YUN T S，LAU Y M，et al.DEM simulation on soil creep and associated evolution of pore characteristics[J].Computers and Geotechnics，2012，39(1):98－106.

[19]王涛，吕庆，李杨，等.颗粒离散元方法中接触模型的开发[J].岩石力学与工程学报，2009，28(增刊2):4040－4045.WANG Tao，LU Qing，LI Yang，et al.Development of contact model in particle discrete element method[J].Chinese Journal of Rock Mechanics and Engineering，2009，28(Supp.2):4040－4045.

[20]赵文，曹平，章光.岩石力学[M].长沙:中南大学出版社，2006:12－14.

[21]Itasca Consulting Group.PFC3D(Particle flow code in 3 dimensions) theory and background[R].Minneapolis，USA:Itasca Consulting Group，2008:33－39.

[22]何涛.沙坝土矿红页岩蠕变特性室内试验与数值模拟研究[D].长沙:中南大学，2012.

[23]赵东宁，李明霞，张艳霞，等.微裂隙对灰质泥岩强度的影响分析[J].人民长江，2013，44(22):86－89.ZHAO Dong-ning，LI Ming-xia，ZHANG Yan-xia，et al.Influence study of micro-fissure on strength of lime mudstone[J].Yangtze River，2013，44(22):86－89.