圆形巷道围岩层裂或板裂化的等效连续介质模型 及侧压系数的影响

王学滨，伍小林，潘一山

（辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000）

1 引 言

巷道围岩的层裂现象是指巷道周围的岩石在压应力的作用下，岩石被劈裂成薄板的一种奇特现 象[1－4]，这种现象在矿柱中[5]也十分常见。层裂现象与岩爆、冲击地压也密切相关。在岩爆和冲击地压的机制分析方面，有些理论模型认为巷道附近的岩石中包含着等间距的裂纹或薄板（薄板位于两裂纹之间）[1－4]，据此可对模型进行某种简化和处理。在岩爆后的巷道围岩中，经常能发现各种形状及尺寸的凹坑，常称这些凹坑为岩爆坑，其中以V 形或三角形的岩爆坑最为常见。文献[1]认为，V 形岩爆坑形成的主要过程包括4 个阶段：损伤阶段、扩容阶段、劈裂成板的剥离阶段及稳定阶段。在第3 阶段，断裂过程区的发展导致了薄板的形成，它们通过剪切、劈裂和屈曲产生，其厚度在1～5 cm 之间。文献[1]还认为，这些薄板之间的凝聚力已损失高达100%，如果不对巷道施加及时的支护，这些像一叠卡片似的薄板易于在重力作用下从开挖边界剥落下来，从而形成呈V 形的岩爆坑。文献[2]认为，紧挨着V 形坑的尖端的岩体中也包含着若干薄板，但凝聚力的损失并不完全，因而它们仍然粘在一起。

在岩石破坏过程的数值模拟方面，通常都采用非均质模型，目前，已经取得了一些重要的研究进展[6－7]，在这些模型中，岩石均被处理成连续介质，只不过单元的力学参数是随机分布的。但岩石不仅是一种非均匀介质，也是一种非连续介质，岩石内部存在各种形式的孔洞和微裂纹等微观缺陷，仅将岩石视为非均匀介质，并不能全面地反映岩石的物理本质。以砂岩为例，砂岩属于典型的颗粒体材料，它们与散体不同，是由圆形或椭圆形的砂粒通过胶结物的胶结而形成的岩石[8]。这类岩石通常都具有较高的孔隙度，以砂岩为例，孔隙度可达30%以上，一般为10%～20%。

本文将岩石视为颗粒堆积体，其中包含众多的颗粒和一定的孔隙。如果颗粒之间发生接触，则在颗粒之间设置界面。另外，将颗粒和界面均视为连续介质，以区别于离散元方法。本文提出的方法在FLAC（连续介质快速拉格朗日分析）中得以实现。本文的特点在于采用处理连续介质的方法处理包含一定孔隙的等效连续介质。

2 方法介绍

本文方法包括6 个步骤：（1）生成颗粒堆积体；（2）切割模型；（3）确定界面的尺寸和位置；（4）剖分单元，修正界面单元；（5）赋予单元本构关系；（6）加载、计算、结果保存。

2.1 生成颗粒堆积体

在本文中，岩石被视为众多圆形颗粒在重力作用下堆积而成的颗粒堆积体。在某一范围A 内随机生成颗粒堆积体模型。颗粒为圆形，其堆积过程从下到上进行。为了模拟单个颗粒在重力作用下的自然堆积过程，需要保证任一颗粒在重力作用下能静力平衡，为此任一颗粒至少要与其下方或两侧的另外两个颗粒接触，此外，拟堆积的任一颗粒的中心应位于已堆积的两个颗粒的中心之间。在颗粒堆积过程中，需要自动记录和修改当前颗粒堆积体的最上一层（称之为表层）的颗粒，例如，当任一新的颗粒被堆积在颗粒堆积体上时，需要对当前表层各个颗粒的记录进行修改（删除旧记录和增加新记录）。在堆积一个新的颗粒时，将优先堆积在颗粒堆积体表层的最低凹处附近，目的是使颗粒堆积体的表层尽可能平坦。

2.2 切割模型

在A 区域内切割一个正方形区域B。选择中心位于正方形区域B 之内的颗粒的集合体作为试样或计算模型。在切割试样时，不允许一个颗粒被切成两半。尽管有些颗粒的中心位于区域B 之内，但如果它们与B 内的其他颗粒仅有一个接触点，则需要删除它们。在本文中，区域B 的尺寸为1 m×1 m，区域B 内的颗粒数目为2 795 个。颗粒的最大半径、最小半径分别为0.012 5 m、0.006 25 m。假设颗粒的半径服从均匀分布，对颗粒半径的统计结果见图1，横坐标为在颗粒半径的范围内将其等分成20 个区间，纵坐标为落在各半径区间内的颗粒数目占颗粒总数的百分比。

图1 颗粒半径分布图 Fig.1 Distribution of granular radius

2.3 确定界面的尺寸和位置

如果颗粒堆积体中的两个颗粒发生接触或者二者之间的距离非常近，则需要在二者之间规定界面。设置界面的目的是模拟颗粒之间的滑动和接触。界面尺寸的确定原则如图2(a)所示。在图2(a)中，圆O1和O2代表2 个圆形颗粒，其半径分别为r1和r2，假定r1＜r2。点O1、O2之间的距离设为d。假定这2个颗粒之间存在一个界面，用矩形ABCD 表示，界面的中心点为O，界面的长和宽分别为l 和b，界面的长度方向是指垂直于O1O2的方向，而界面的宽度方向是指平行于O1O2的方向。界面的位置（由O点确定）及尺寸（由l 和b 确定）的确定需要根据已知的参数，例如点O1的坐标（x1，y1）、点O2的坐标（x2，y2）、r1、r2、d 及将来剖分的正方形单元的尺寸（边长）R。

在本文中，当2 个颗粒之间的距离d0小于一个正方形单元的尺寸R 时，假定这2 个颗粒发生接触，此时需要在二者之间设置界面，即

图2 两个颗粒及它们之间的界面 Fig.2 Two granular and an interface between them

由图2(a)可以发现，界面的宽度b 应大于R。b的值的选取需要凭借经验，可以先给b 一个稍大一些的值，例如取R 的2 倍，在将颗粒和界面剖分成单元之后，如果必要可再对颗粒和界面单元进行适当的修正，以保证界面不致太厚，也不致太薄（致使无法将两个颗粒的单元完全隔开）。

假定界面的2 个角点A、B 位于较大的颗粒的圆周上。设界面的中心点O 的坐标为（x，y）。由于O 点在线段O1O2上，因此，只需确定x 即可。这里假定，x 与O1点的横坐标x1、O2点的横坐标x2有下列简单的关系：

在此，可以对式（2）的依据进行简略的讨论。由式（2）可以发现，当r1→0 时，x→x1。也就是说，当小颗粒很小时，界面ABCD 的中心点O 向小颗粒的中心趋近。这么做的目的是为了确保界面能将2个相差较大的颗粒联接起来（见图2（b））。否则如果r1→0 时，让x→x2（不再遵循式（2）），则完全可能出现图2(c)的情形。此时，界面没有起到联结2 个颗粒的作用，这样，这2 个被判为相互接触的颗粒由于界面位置设置不当将不发生任何作用。

假设线段AB 与线段O1O2相交于E 点，线段EO2的长度为a，线段OO2的长度为c，则c 为

在ΔAEO2中，根据几何关系，可得

这样，利用式（3）、（4）可以得到界面长度l的表达式。界面的位置和尺寸确定之后，需要将界面和颗粒都剖分成正方形单元，其余的部分是孔隙。

2.4 剖分单元，修正界面单元

为了采用FLAC 模拟等效连续介质，颗粒O1、O2和界面ABCD 都将被离散为正方形单元。因此，实际上，图2 中的矩形ABCD 和颗粒O1、O2仅是离散后的界面和颗粒的轮廓。位于矩形ABCD 内的正方形单元的集合构成了真实的界面，这些单元称之为界面单元，而位于圆内的单元称之为颗粒单元。

如果界面单元不能将2 个颗粒的单元完全隔开，则2 个颗粒的一些单元可能共用一些节点，由于在计算中节点无法分离，这对于模拟2 个颗粒沿界面的相对运动是不利的，因为这将使颗粒无法做相对运动或相对运动受到阻碍。因此，适当修正是必要的。如果界面单元显得过厚，则要将那些多余的界面变成颗粒单元或孔隙；如果界面显得过薄，则需要适当将颗粒单元或孔隙当成界面单元。

将颗粒和界面剖分成单元，并经过适当修正之后，经过统计可以得到下列数据：计算模型水平及垂直方向上的最大单元数约为500，颗粒单元的总数为174 937，界面单元的总数为29 353。据此，计算模型的孔隙度大约为0.182 8（（500×500－ 174 937－29 353）/（500×500））。对于Berea 砂岩，其孔隙度在17%～25%之间[9]，本文的孔隙度与其下限更接近。

2.5 赋予单元本构关系

颗粒单元被视为各向同性的弹性材料，弹性模量取为26.52 GPa，泊松比取为0.21，而界面单元在破坏之后被视为摩尔-库仑材料，界面单元的凝聚力取为2 MPa，内摩擦角取为30°。在界面单元破坏之前，其本构关系与颗粒单元一致。

2.6 加载、计算、结果保存

加载步骤如下：首先，给定本构关系、边界条件及模型四周的围压，对未开挖的颗粒堆积体模型进行计算，阻尼由FLAC 自行施加，直到达到静力平衡状态。在本文中，若最大失衡力小于0.05 N 时则停止计算，认为模型已经达到了静力平衡状态。由于模型尺寸较小，因此，在计算中没有考虑重力。

然后，为了模拟开挖，将巷道内部的单元删除。以模型的中心做1 个直径为模型尺寸1/3 的圆。将那些位于圆内的颗粒和界面单元删除，以模拟巷道的一次性开挖。为了使巷道的表面看起来更圆一些，这里允许将颗粒和界面一分为二，仅删除那些位于圆内的单元。

图3 颗粒堆积体模型 Fig.3 Model of circle granular material

开挖后的模型见图3，从中可以看到两种单元和孔隙。巷道开挖之后，模型的四周仍然受到压应力作用。由于开挖卸荷，因而在巷道附近的节点上将存在失衡力，模型将不再处于静力平衡状态。

最后，对开挖后的模型重新进行计算，阻尼由FLAC 自行施加。对于每个方案，从开挖开始之时开始计算，一直计算到4×104步。从失衡力-时步的演变规律来看，各个方案均达到了静力平衡状态。

本文采用了4 个计算方案（见表1）。方案1～4的模型四周的压应力不同，而其余参数全相同。

表1 计算方案 Table 1 Computational schemes

3 计算结果及分析

3.1 剪切应变增量的分布规律

图4 给出了不同方案的剪切应变增量高值区分布的等值线图，是从巷道开挖之后再计算4×104步的结果。在图中，颜色越深代表单元的剪切应变增量越大。本文中随后的其他图片也是这一时步时的结果，不再赘述。此时，开挖后的模型已经重新达到平衡。由图还可以发现，当水平方向的应力小于垂直方向的应力时，在巷道两侧的围岩中出现了较高的剪切应变增量，见图4(a)、(b)。与方案2 相比，方案1 的剪切应变增量的高值区严格位于巷道的两侧，而在巷道的顶部及底部的围岩中剪切应变增量较低。当模型的四周的压应力相等时（见图4(c)）围岩中的剪切应变增量高值区的等值线图呈圆环状，位于巷道附近。图4(d)的结果与图4(a)刚好相反。在图4(d)中，围岩中的剪切应变增量的高值区位于巷道顶部和底部，长条形的剪切应变增量等值线基本上都在水平方向上。

计算表明，方案1～4 的最大剪切应变增量分别为 2.412 1×10-3、1.891 3×10-3、1.837 6×10-3及 2.456 1×10-3。可以发现，方案1、4 的最大剪切应变增量的值较高，而方案3 的值较小。另外，在图4 中，为了能更清晰地显示出长条形的剪切应变增量等值线图，未显示大于2.5×10-4的剪切应变增量。

图4 不同方案的剪切应变增量（当时步=4×104） Fig.4 Contour maps of shear strain increment at different schemes (timesteps=4×104)

3.2 垂直方向应力的分布规律

图5 给出了4 方案的垂直方向应力的高值区分布的等值线图，该图仅显示了大于0 的垂直方向的应力分布。在FLAC 中，以压应力为负，而以拉应力为正。为了能更清晰地显示出垂直方向的应力为拉应力的单元的拉应力分布规律，在图中未显示大于1.5 MPa 的垂直方向应力。在图中，只有单元的垂直方向应力是拉应力才被显示出来。可以发现，在图5(a)～(c)中，垂直方向应力为拉应力的单元的分布规律具有一定类似性。随着水平方向应力的增加，垂直方向应力为拉应力单元的数目有增多的趋势；垂直方向应力的高值区在巷道的顶部及底部形成的V 形坑越来越明显。图5(d)的计算结果与图5(a)～(c)差别较大，如果显示方案4 的水平方向的应力分布（若只显示大于0 部分，且旋转90°后再比较），就不应与图5(a)有大的差别。

图5 不同方案的垂直方向应力（当时步=4×104） Fig.5 Vertical stress distributions at different schemes (timesteps=4×104)

由计算结果可以发现，方案1～4 大于0 的垂直方向应力的最大值分别为2.89、2.62、3.03、3.03 MPa。这些值都与所施加的最大压应力值（3.00 MPa）相近。这说明，尽管在模型的四周施加的是压应力，但由于岩石被处理为等效连续介质，由于颗粒之间的相对运动规律较为复杂，因此，在模型中诱发出了与所施加的压应力接近的拉应力。众所周知，通常岩石的抗拉强度是抗压强度的1/10～1/4，因此，岩石在压应力作用下产生的如此高的拉应力是易于将岩石拉坏的。笔者的数值计算经验表明，采用均质或非均质模型，在岩石内部不可能产生如此大的拉应力。尽管在压缩载荷作用下，在岩石内部将诱发出拉破坏在岩石力学界被广泛承认，但试图模拟出拉应力或者超过一定值的拉应力一直是困难的。本文的方法从等效连续介质的角度出发，得到了不同于均质和非均质模型的计算结果。

3.3 最小主应力的分布规律

图6 给出了不同方案的最小主应力的高值区分布的等值线图。应当指出，为了使各方案的最小主应力分布规律显示得更清晰，在图6 中仅显示了那些小于0 的最小主应力，而且未显示超过10 MPa的压应力。在FLAC 中，最小主应力一般为绝对值较大的负值，代表很高的压应力。可以发现，图6的计算结果与图4 中的剪切应变增量分布的计算结果比较类似，图6(a)、(d)的结果也刚好相反。计算表明，方案1～4 的最小主应力的最小值（很高的压应力）分别为-17.01、-15.36、-18.79、-24.29 MPa，而最小主应力的最大值（较高的拉应力）分别为3.34、3.02、2.80、2.84 MPa。由此可以发现，模型中的最大压应力远超过模型边界上施加的压应力（最大值为3 MPa），可以高达8 倍，例如，在方案4 中，最小主应力的最小值为-24.29 MPa，其大小是3.00 MPa 的8.1 倍。在如此高的压应力作用下，模型内部将不可避免地发生损伤。模型中的最大压应力与边界上所施加的压应力的比值必然与材料的微观结构有关。如果将本来是非连续介质视为连续介质，纵然在设计中考虑了一定的安全系数，也可能无法避免存在的安全隐患。

图6 不同方案的最小主应力（当时步=4×104） Fig.6 The minor principal stress distributions at different schemes (timesteps=4×104)

另外，方案1～4 的最小主应力的最大值在3 MPa附近，也近似等于在模型的边界上所施加的最大的压应力的值。

3.4 最大主应力的分布规律

图7 给出了不同方案的最大主应力的高值区分布的等值线图。为了使各方案的最大主应力分布规律显示得更加清晰，图中仅显示了那些大于0 的最大主应力，而且未显示大于1.5 MPa 的最大主应力。在FLAC 中，最大主应力通常为正值，代表较高的拉应力。可以发现，与最小主应力的高值区呈环向分布不同（见图6），最大主应力的正高值区基本上都位于径向。计算表明，方案1～4 的最大主应力的最大值分别为4.69、4.25、3.54、4.35 MPa，而最小值分别为-4.15、-3.90、-3.59、-3.78 MPa。还可以发现，对于方案2、3，最大主应力的正高值区都位于巷道附近，而对于方案1、4，最大主应力的正高值区已经深入到围岩内部。这是由于模型边界上所受的2 个方向上的压应力之差较大的缘故。

图7 不同方案的最大主应力（当时步=4×104） Fig.7 The maximum principal stress distributions at different schemes(timesteps=4×104)

3.5 一行单元应力的分布规律

图8给出了方案1～4的一行单元的水平和垂直方向应力的分布规律。这一行单元选在巷道左侧的围岩中，而且取在通过巷道直径的水平方向上。这一行单元仅是实体单元，即包括界面和颗粒单元，不包括孔隙。图8 的横坐标为从左到右的单元数，第1 个单元表示离模型左边界最近的单元，而第139个单元表示离巷道表面最远的单元。因此，水平方向应力相当于径向应力rσ ，而垂直方向应力相当于环向应力θσ 。由图可以发现，rσ 、θσ 的分布规律和以往的均质模型得到的解析解或数值解有重大的差别。在一些位置的θσ 的值非常高，因而，二者的分布曲线具有多个峰值，既不光滑，也不连续。

图8 不同方案的巷道左侧一行单元应力分布 （当时步=4×104） Fig.8 Stress distributions of a row of elements in surrounding rock at left side of a circular tunnel at different schemes (timesteps=4×104)

由图8 可以发现，在巷道的表面，rσ 为0，而在模型的左边界上，rσ 与所施加的应力相等。这些结果的正确性是显然的。在方案1～3 中，随着Ph的增加，小于0 的θσ （表示压）的峰值的绝对值下降，而小于0 的rσ 的峰值有上升的趋势。方案4 的rσ 、θσ 的分布与方案1～3 差别较大，这是由于对于方案4 而言，在模型左、右边界上施加的压应力的值大于在模型上、下边界上施加的压应力值。如果显示方案4 的位于拱顶的一列单元的水平和垂直方向应力的分布，则不会与方案1 相比有大的差别。另外，由图可以发现，在离巷道附近的围岩中，有许多单元的rσ 为正值，这表明这些单元受到的径向应力为拉应力，而且具有多个峰值；与方案2、3相比，方案1 中的拉应力的范围最大，而且值最高。拉应力的分布具有多个峰值的特征应该与围岩发生层裂有关。

3.6 讨论

（1）图9(a)、(b)给出了2 个方向不等压条件下，巷道围岩的破坏形式[1－2，5]。这与本文方案1、2 的计算结果较为一致。与图4(a)、(b)、(d)中的剪切应变增量的高值区的分布相比，图6(a)、(b)、(d)中的最小主应力高值区的分布的结果与这些文献给出的结果更为一致。环向的压应力高值区和剪切应变增量高值区似乎具有一定的间距，该间距应该与颗粒的平均尺寸有关，这一点需要进一步研究。

（2）众所周知，在单向压缩条件下，岩石易于劈裂，图8 的结果表明，在巷道表面附近，一些单元的径向应力居然还是拉应力，这更加剧了单元发生沿上述两种高值区的劈裂。这应该是层裂现象的根源，图9(c)给出了一个岩石试件在一向压应力及另一向拉应力作用下发生拉破坏的示意图。

（3）在本文中，由于将颗粒和界面都剖分成了尺寸相同的单元，因此，采用的模拟等效连续介质的方法计算量较大。一般认为，这种方法不适于工程实际问题的分析，更适于破坏机制的分析。然而，以离散元为基础的颗粒流方法也常见被用于工程实际问题的分析之中。一些作者所选取的颗粒尺寸都普遍大于岩石颗粒的真实尺寸。以砂岩为例，其平均颗粒半径可以达到0.25 mm[9]。文献[10]在数值计算中选取的颗粒半径为0.2 m；文献[11]在数值计算中选取的颗粒半径为0.01 m，这与本文接近；文献

[12]在数值计算中选取的颗粒半径约为1.5 cm。在计算中，采用大颗粒进行模拟，一方面是不得已而为之，另一方面或许可以在某种程度上近似岩石的非连续性，或者认为，大颗粒与节理岩体中的岩块具有某种等效性。从这一方面讲，利用比真实颗粒尺寸大得多的颗粒进行模拟，是可以用于工程尺度问题的定性分析的。此时，颗粒仅是岩体中完整部分的某种近似，而颗粒之间的接触、错动以及孔隙在某种程度上代表了岩体中的非连续性的某些侧面。由于本文中颗粒的尺寸大于砂岩的实际颗粒的尺寸，因此，似乎应该认为本文模型是节理岩体的某种等效或近似。尽管本文中颗粒的尺寸大于砂岩颗粒的真实尺寸，但从图6(a)、(b)、(d)给出的最小主应力高值区的分布和图9(a)、(b)给出的破坏区域可以发现，二者的范围大体相当。

图9 围岩层裂的一些结果[1－2, 5]及岩样受拉压应力 发生拉破坏的示意图 Fig.9 Some results about the exfoliation phenomenon of the surrounding rock [1－2, 5] and schematic about the tensile failure of a rock specimen under combination of compressive and tensile stresses

（4）通过计算可以发现，模型中个别单元的应力或应变的值会很高，这与这些单元附近的微结构的细节有密切的关系，例如孔隙度及颗粒的排列方式等，但其余的绝大部分单元的应力都在某一范围之内，这应该得到重点的关注。如果将所有的单元的应力都显示在一张图片上，个别的高值就将主流的与高值相比的低值掩盖了，因此，在图4～7 种舍弃了个别的高值。

（5）限于本文篇幅，未能研究下列问题：界面单元强度的影响，孔隙存在及填充后的计算结果的差别，颗粒尺寸及排列方式的影响等。在文献[13]中比较详细地介绍了颗粒的堆积算法，研究了圆形洞室尺寸效应问题，关于其他形状洞室的研究也有必要开展。

4 结 论

（1）将岩石视为等效连续介质后，所得的数值计算结果与均质和非均质模型的计算结果相比，具有重大的差异。本文的计算结果显示，在模型中可以诱发出与所施加的压应力大小相当的拉应力，模型内部的压应力远高于在模型边界上所施加的压应力的值。上述结果无论从哪一方面讲，过去的均质和非均质模型的计算结果都可能存在潜在的风险。因此，非连续性有必要在设计中予以考虑。

（2）模型中的最小主应力和剪切应变增量的分布具有一定的类似性。当水平方向和垂直方向的压应力不相等时，二者的高值区位于巷道的两侧或顶部和底部，而在静水压力条件下，最小主应力和剪切应变增量的高值区的分布近似呈环形，主要位于巷道的周围。最小主应力和剪切应变增量的高值区分布基本平行，这些位置是将来裂纹出现的位置。围岩层裂现象的原因被认为是环向的高压应力和径向的高拉应力的共同作用。

[1] MARTIN C D. The effect of cohesion loss and stress path on brittle rock strength[J]. Canadian Geotechnical Journal， 1997， 34(5)： 698－725.

[2] BAŽANT Z P， LIN F B， LIPPMANN H. Fracture energy release and size effect in borehole breakout[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 1993， 17(1)： 1－14.

[3] 冯涛， 潘长良. 洞室岩爆机理的层裂屈曲模型[J]. 中国有色金属学报， 2000， 10(2)： 287－290. FENG Tao， PAN Chang-liang. Lamination spallation buckling model for formation mechanism of rock burst[J]. The Chinese Journal of Nonferrous Metals， 2000， 10(2)： 287－290.

[4] 王敏强， 侯发亮. 板状破坏的岩体岩爆判别的一种方法[J]. 岩土力学， 1993， 14(3)： 53－59. WANG Min-qiang， HOU Fa-liang. A method to determine rock burst in rock mass with plate-form damage[J]. Rock and Soil Mechanics， 1993， 14(3)： 53－59.

[5] KLERCK P A， SELLERS E J， OWEN D R J. Discrete fracture in quasi-brittle materials under compressive and tensile stress states[J]. Computer Methods in Applied Mechanics Engineering， 2004， 193： 3035－3056.

[6] TANG C A， THAM L G， LEE P K K， et al. Numerical studies of the influence of microstructure on rock failure in uniaxial compression-Part II： Constraint， slenderness and size effect[J]. International Journal of Rock Mechanics and Mining Sciences， 2000， 37(4)： 571－583.

[7] FANG Z， HARRISON J P. Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks[J]. International Journal of Rock Mechanics and Mining Sciences， 2002， 39(4)： 443－457.

[8] 张凯， 周辉， 胡大伟， 等. 弹塑性条件下岩土孔隙介质有效应力系数理论模型[J]. 岩土力学， 2010， 31(4)： 1035－1041. ZHANG Kai， ZHOU Hui， HU Da-wei， et al. Theoretical model of effective stress coefficient for rock/soil-like porous media under elastoplastic condition[J]. Rock and Soil Mechanic， 2010， 31(4)： 1035－1041.

[9] HAIMSON B C. Borehole breakouts in Berea sandstone reveal a new fracture mechanism[J]. Pure and Applied Geophysics， 2003， 160(5－6)： 813－831.

[10] 王涛， 吕庆， 李杨， 等. 颗粒离散元方法中接触模型的开发[J]. 岩石力学与工程学报， 2009， 28(增刊2)： 4040－4045. WANG Tao， LÜ 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.

[11] 周健， 邓益兵， 贾敏才， 等. 基于颗粒单元接触的二维离散-连续耦合分析方法[J]. 岩土工程学报， 2010， 32(10)： 1479－1484. ZHOU Jian， DENG Yi-bing， JIA Min-cai， et al. Coupling method of two-dimensional discontinuum-continuum based on contact between particle and element[J]. Chinese Journal of Geotechnical Engineering， 2010， 32(10)： 1479－1484.

[12] HAZZARD J F， YOUNG R P. Moment tensors and micromechanical models[J]. Tectonophysics， 2002， 356： 181－197.

[13] 伍小林， 王学滨， 潘一山. 含不同半径孔洞的颗粒体模型的力学行为数值模拟[J]. 中国地质灾害与防治学报， 2011， 22(1)： 107－113. WU Xiao-lin， WANG Xue-bin， PAN Yi-shan. Numerical simulation of mechanical behavior of the granular model with a hole with different radii[J]. The Chinese Journal of Geological Hazard and Control， 2011， 22(1)： 107－113.