考虑扩容的围岩局部破坏过程中声发射及能量释放的数值模拟

王学滨，张春野，潘一山，孙秀东

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

在巷道发生岩爆之后，经常能观察到巷道的围岩发生V形坑式破坏［1-5］，V形坑一般有两个以上。在静水压力条件下，V形坑为4个［1-3］。围岩的V形坑式破坏是一种典型的局部化现象，基于一维轴对称弹塑性模型难于对其进行分析。

最近，一些文献采用数值计算方法，研究了V形坑的形成过程及多种因素的影响：文献［6］在静水压力条件下，模拟出了巷道围岩中的4个对称的V形坑，这与岩爆后现场的观测结果及实验结果［1-3］比较吻合，还研究了巷道的尺寸效应，定性地解释了岩爆的级别随着巷道尺寸的增大而增大的调查统计结果［7］；文献［8］考虑了围压和孔隙压力的联合作用，在高孔隙压力条件下，围岩中出现了口小腔大的空洞现象，这与煤和瓦斯突出后形成的梨形空洞［9］比较类似；文献［10］采用了徐林生及Russense判据（环向应力与单轴抗压强度之比）作为岩爆判据，研究了围压对3种岩爆区（严重、中等及轻微岩爆区）分布规律的影响。在上述数值模拟研究中，采用了下列处理手段：（1）将围岩剖分成尺寸相同的矩形网格，这对于模拟局部化现象是有利的，可以避免辐射状网格的不好结果（离巷道越远，剪切带的尺寸越大；或者不能模拟出局部破坏现象，仅能模拟出巷道周围的圆环形的塑性区）；（2）“先加载，后挖洞”的计算顺序，这与在受载的地下岩石中挖掘巷道的工况较为一致。

在上述数值模拟研究中，都采用了带拉伸截断的剪破坏应变软化模型，扩容角均被设置为零。文献［11］研究了扩容角对V形坑形态等方面的影响，发现当扩容角取值小时，V形坑尖且小，这与文献［3］中的观测结果较一致；当扩容角取值大时，V形坑钝且大，这与文献［2］中的观测结果较一致。虽然扩容角的变化范围并不大，例如对于岩土材料和混凝土，其范围一般在0°～20°，但是其影响却是巨大的，吸引了许多研究人员的关注。研究发现，扩容角对地基承载力、边坡稳定性、结构的破坏模式、剪切带的倾角及宽度、失稳的前兆等方面均有不可忽视的影响［12-17］。扩容的机制一般认为是由于微裂纹的张开已经超过了闭合［18］，或者是由于滑动块体在凹凸面上的抬升［19］。在高地应力条件下挖掘巷道时，岩土材料由于破坏而引起的扩容将给巷道的维护和支护带来一定的困难［4，11，15］，因而有必要研究扩容的多种影响，例如对巷道围岩中的声发射及能量释放的影响，这方面的研究目前还未见报道。

本文介绍了FLAC中扩容角的引入方法，阐述了采用数值方法模拟声发射及能量释放的思路，在不同扩容角条件下，研究了巷道开挖之后围岩中的声发射及能量释放的演变特点，解释了未设定抗拉强度与塑性拉伸应变之间的关系，但是实际上却发生拉伸应变软化的原因。

1 基本理论及方法

1.1 FLAC中关于扩容角的基本理论

在FLAC中，一个单元是否发生破坏，取决于其应力状态与初始屈服面之间的关系。初始屈服面包括剪切和拉伸屈服面，它们由初始的内聚力、内摩擦角及抗拉强度所决定，与扩容角并没有关系。

但是，一个单元发生破坏之后，其发生的塑性流动的大小，即塑性应变增量，却与扩容角有关，这由塑性流动法则决定。如果塑性流动法则中的塑性势函数和屈服函数等同，则称之为相关联的流动法则，否则，则称之为非相关联的流动法则。

在FLAC中，塑性势函数也包括两部分，即剪切和拉伸部分，其中塑性剪切势函数gs的表达式为：

式中：σ1和σ3——第1和第3主应力；

Nψ与扩容角ψ有关。

由式（2）可以发现，ψ越大，则Nψ越大。

塑性流动法则建立了塑性势函数对各主应力的偏导数和各塑性主应变增量之间的关系，对于剪切流动而言，有下列关系

式中：λs是一个常数，Δεpis是各塑性主应变增量。由于gs仅和σ1及σ3有关，而与σ2无关，因此可以得到

由式（4）及（5）可以发现，仅第3塑性主应变增量与扩容角ψ有关。ψ越大，第3塑性主应变的值越高，即单元的侧向塑性变形或不可逆剪胀越明显。对于巷道围岩中的一个单元而言，第1塑性主应变增量应位于环向，而第3塑性主应变增量应位于径向，这与σ1位于环向、σ3位于径向是类似的，因此可以期待，当ψ高时，由于破坏引起的扩容的后果将造成巷道断面的极大缩小，实际上也是如此［4，11，15］。

1.2 声发射及能量释放的统计

本文中的声发射率和实验中的声发射数或振铃数具有一定的类似性。如果一个单元刚发生破坏，则被视为一个事件，此后，其不再被视为一个事件。由于将一个刚发生破坏的单元视为一个事件，因此事件的尺度就是一个单元的尺寸，这与实际情况是有差距的，实验中的一个事件可能涉及一定的尺度。

在实验中，在一定的时间间隔之内对事件的数目进行统计，而在本文中，对一定的时步间隔之内的事件的数目进行统计。时步与时间是不同的概念，时步并没有单位，是指一个计算循环，完成从几何方程到本构方程再到运动方程的计算。为了完成一个问题的计算，例如从加载至平衡，往往需要计算大量的时步，才能达到令人满意的结果。对于一个具体的问题，也可以采用一些方法对时步与时间的关系进行校准，即确定多少时步相当于1秒。然而，一般为了研究方便，直接给出不同时步结果，这并不妨碍对物体变形及破坏规律的理解。

FLAC为了加快计算的速度，每10个时步才在窗口中显示一次计算结果，因此统计每10个时步内的事件并将其显示出来是方便的，这样既不会消耗过多的时间用于显示结果，一般也可以保证在10个时步之内有足够的事件用于统计。

如果将每10个时步之内的声发射率加在一起，则可以得到声发射率的累计数，即破坏单元的总数。为了进一步了解由剪切和拉伸破坏所引起的事件的比例，可以将声发射率和声发射累计数细分为剪切和拉伸部分。在声发射率-时步曲线上，如果声发射率的值高且比较密集地分布，则说明事件比较活跃，否则，如果声发射率的值很小或为零，则意味着没有新的事件发生。在声发射累计数-时步曲线上，事件比较活跃体现在声发射累计数的快速上升，没有新事件发生则表现为声发射累计数保持不变。

任一10个时步之内的声发射统计的计算量很小，但是却存在一定的缺陷，例如，不能区分出这些事件的大小，即无论是一个大事件，还是一个小事件，都被计为1，这显然不能反映出事件的全部特征。为此，本文又计算了破坏单元释放的能量，即释放的弹性应变能。在这里，事件不是被定义为一个刚发生破坏的单元，而是定义为释放弹性应变能的破坏单元，如果一个破坏的单元总在释放弹性应变能，则它就总被视为一个事件。也就是说，在每10个时步之内，只要一个破坏单元正在释放弹性应变能，它就被视为一个事件，就将它释放的弹性应变能记录下来，遍历所有的事件，就可得到围岩释放的弹性应变能，在本文中称之为弹性应变能释放率，单位仍然同弹性应变能（J）。对任一10个时步之内的围岩释放的弹性应变能求和即可得到累计的弹性应变能。单元释放的弹性应变能由破坏单元储存的弹性应变能的降低量描述［20］。降低量通过单元储存的弹性应变能前、后作差获得，要求之差大于零。

虽然，围岩释放的弹性应变能的计算比较繁琐，但是和声发射相比，更容易揭示出事件大小上的差别。与声发射随时步的演变规律类似，如果弹性应变能释放率较高且密集，则意味着许多单元同时在释放能量，因此事件非常活跃，在弹性应变能的累计曲线上则表现为快速上升。类似地，为了研究上的需要，也有必要区分出剪切和拉伸破坏释放的能量。

由于任一10个时步结束时，才统计出上述各种结果，因此上述曲线中的时步是指这10个时步结束时当前的总时步。统计声发射及能量的释放规律，可以更深入地了解事件（或局部失稳）的发生、发展及终止规律。上述统计，不局限于平衡过程之中，对于失稳过程同样适用。

2 计算模型、计算参数及计算过程

矩形平面应变模型见图1（a），模型长及高均为1m，未开挖的模型被划分为40000个面积相同的矩形单元。模型四周受到静水压力p=7.5MPa作用，即模型受到的水平及垂直方向的压缩应力相同。数值计算在小变形模式下进行。

采用3个计算方案，各方案的扩容角ψ不同（方案1至3的ψ分别为0°、10°及20°），其余参数相同。各方案的孔径为0.31m。在弹性阶段，岩石的本构关系取为各向同性的线弹性模型，弹性模量取为26.6GPa，泊松比取为0.21。峰值强度后岩石的破坏准则取为带拉伸截断的剪破坏莫尔-库仑模型，抗拉强度取为2MPa。在峰后，内摩擦角及内聚力与塑性剪切应变的关系先是降低（应变软化），然而分别达到残余内摩擦角及残余内聚力，当塑性剪切应变为零时的内摩擦角及内聚力分别称之为初始的内摩擦角及内聚力。内摩擦角-塑性剪切应变关系与文献［21］中密实岩石的软化规律相同。内聚力-塑性剪切应变关系与文献［21］中密实岩石的软化规律类似，但是初始的内聚力取为4MPa。

图1 模型的几何特征及边界条件Fig.1 Model geometry and boundary conditions

计算步骤如下：首先，建立未开挖的模型（图1（a）），给定本构关系、边界条件及加载条件（模型四周的围压），为了使模型尽快达到平衡条件，采用了与文献［6，8，10-11］不同的处理方式，按平面应变条件初始化模型内部单元的各个方向的应力值。在本文中若最大失衡力小于1.5×10-1N则认为已经达到了静力平衡条件。

然后，利用编写的 FISH 函数［6，8，10-11］开挖巷道，即将模型中巷道直径范围内的单元一次性删除（图1（b））。不管巷道开挖与否，模型始终受到围压的作用，由于开挖卸荷而引起了最大失衡力的突增，模型将不再处于静力平衡状态。

最后，对开挖后的模型进行计算，利用编写的各种FISH函数［20］计算声发射率及累计数（与破坏单元总数等同）、能量释放率及能量释放的累计（总量），直到达到一个新的平衡状态。

3 结果分析及讨论

图2给出了当时步=6万时3种破坏区的分布规律，其中巷道附近的黑色区域表示发生破坏的单元。图2（a-c）、图2（d-f）及图2（g-i）分别是ψ=0°、10°及 20°时的结果。图 2（a，d，g）、图 2（b，e，h）及图2（c，f，i）分别是破坏单元、剪破坏单元和拉破坏单元的分布规律。图3（a-c）分别给出了3种破坏单元的数目（声发射数的累计）随时步的演变规律；图3（d-e）分别给出了两种释放的能量的累计随时步的演变规律；图4（a-d）分别给出了4种量的（变化）率随时步的演变规律。

3.1 破坏区的分布规律

由图3及图4可以发现，当时步=6万时，各方案均已达到了静力平衡状态，因此，图2给出的结果是巷道开挖之后围岩再次平衡后的结果。由图2可见，当围岩再次平衡时，在围岩中产生了4个对称的V形坑，这与静水压力条件下有关的实验及现场观测结果［1-3］吻合；剪切破坏单元形成的V形坑较深，而拉伸破坏单元形成的V形坑较浅。在剪切破坏单元形成的V形坑尖端的位置，并无拉伸破坏单元出现，这与这一位置具有强烈的压应力集中（σ1的值较高）有关，不易发生拉伸破坏。

图2 当6万时步时三种破坏区的分布Fig.2 Distribution of three kinds of failed zones at 60000 time steps

由图2（b，e，h）可以发现一个重要的现象：随着ψ的增加，V形坑变大、变钝、变深，下面对其原因进行解释。图5示意地显示了当V形尖坑和钝坑出现后，剪切带方向与径向应力（σ3）之间夹角的差异，即α1＜α2。这说明，随着ψ的增加，剪切带方向与σ3之间的夹角增加，这可由著名的Roscoe和Arthur倾角公式［22］进行解释，而Coulomb倾角公式由于不包含ψ，并不能解释上述现象。Roscoe和Arthur倾角公式都表明，剪切带方向与σ3之间的夹角随着ψ的增加而增加。

由图2可以发现，当ψ由0°增至10°时，破坏区的形态变化较大，而当ψ由10°增至20°时，破坏区的形态变化不大。图3（a，c-d）也反映了类似的规律，即当ψ=10°及20°时，上述3种曲线的最终值（当时步=6万时的结果）相差不大。然而，图3（b）的结果并不是这样，ψ=20°时的结果反常地位于ψ=0°及10°的结果之间，观察图2（f，i）可以发现，当ψ=20°时，拉伸破坏区尚未形成完整的V形坑，而当ψ=10°时，拉伸破坏区形成的V形坑已较清晰，因此，ψ=20°时的拉伸破坏单元数会少于ψ=10°时的结果。在图3（e）中，ψ不同时的结果并没有太大的差别。

3.2 声发射及能量释放的演变规律

由图3中各种曲线图可以发现，随着ψ的增加，当时步=6万时各种量的结果并不是严格呈单调的变化规律。也就是说，并不是ψ高时的结果都在最上方，ψ=10°的结果在最上方有3次（图3（b-c，e）），ψ=20°的结果在最上方有两次（图3（a，d））。尽管如此，在通常的情况下，ψ=0°时的结果都在最下方，仅图3（e）的结果是个反例。这意味着ψ高时各种量通常都较高。并且，由图3及图4还可以发现，随着ψ的增加，达到静力平衡时的时步要多一些，即平衡要困难一些；事件发生快速突增要稍晚一些。

图3 不同扩容角时，5种累计量随时步的演变规律Fig.3 Evolution of five kinds of accumulated quantities with time step at different dilation angles

图4 不同扩容角时，4种率随时步的演变规律Fig.4 Evolution of four kinds of rates with time step at different dilation angles

由图3（e）和图4（d）可以发现，有一定的拉伸应变能被释放出来。然而，在本文中，并没有规定抗拉强度的软化规律，即抗拉强度的值被设置为一个定值。这似乎是矛盾的，下面对此现象进行解释。

图5 不同V形坑时的剪切带与第3主应力之间的夹角Fig.5 Angles between the maximum principal stress and the shear band for different shapes of V-shaped notches

图6 规定的抗拉强度和由残余内聚力及内摩擦角决定的最大抗拉强度Fig.6 The prescribed tensile strength and the maximum tensile strength determined by the residual cohesionand internal frictional angle

图6给出了σ1-σ3平面内对称线σ1=σ3上方的初始屈服面（粗实线），包括剪切屈服部分和拉伸截断部分，当岩石发生破坏之后，其抗剪强度参数（内聚力和内摩擦角）将发生相应地降低，直至达到二者的残余值（粗虚线是残余屈服面）。如果由残余的内聚力和内摩擦角决定的最大抗拉强度小于设定的抗拉强度，FLAC程序就会将岩石的抗拉强度调整为由残余的内聚力和内摩擦角所决定的最大抗拉强度，以保证屈服函数的连续性要求，这通过检查C++源代码获知。也就是说，尽管没有规定抗拉强度的软化规则，如果条件满足，FLAC程序会自动地将抗拉强度按脆性模式降低。因此，在此过程中，有一定的拉伸应变能被释放出来的现象可以得到解释。

文献［17］的研究表明，当ψ高时，矩形平面应变试样破坏的前兆变得明显。由于当ψ高时，会引起试样侧向变形的明显增加，因此这一结论显然是正确的。然而，在图3（a-c）及图4（a-b）中并不能看到类似的规律，这似乎也说明只统计事件的数目，而不计事件的大小存在一定的弊端。确实在图4（c-d）及图3（e）中都能发现当ψ不同时，这些量的演变规律有所差异。例如，在图3（e）中，当ψ高时，拉伸应变能累计的突增相对平缓，而当ψ=0°时，其突增非常突然；在图4（c-d）中，当ψ高时，两种应变能释放率在第1次突增时的范围比较宽，因而，V形坑的形成过程要缓慢一些，这对于破坏的预警是有利的；而当ψ=0°，这两种应变能释放率在第1次突增时的范围较窄，而且值都较高，因而即使能觉察到有异常，V形坑会迅速形成，进而可能发生岩爆，因而前兆并不明显。

3.3 ψ变化的深层次原因及引起的复杂现象

扩容的微观机理之一是微裂纹的大量出现。基于连续介质力学的基本原理，仅Δεps3与ψ有关（式（5）），这是从唯象的角度描述物理上的微裂纹的大量出现。

若以一个单元作为分析对象，ψ以十分简单的方式对Δεps3产生影响。然而，在围岩或岩石结构中，情况则要复杂得多。模型或结构之中的单元通过节点相连，一些破坏单元某一方向上的塑性应变的增加，必然会对其附近的其他单元的变形及应力产生影响。随着ψ的增加，这种作用将增强。本文的数值计算结果表明，ψ的增加导致了围岩的破坏花样、声发射、能量释放及前兆的明显变化。这些复杂的变化从不同的侧面刻画ψ的影响，但是，都是基于简单的规则（式（5））。显然，对于复杂的结构，ψ的影响将更加复杂，以致于难于事先预估，因而开展一定的数值模拟研究十分必要。

4 结论

（1）高扩容角时，开挖巷道之后围岩再次平衡时发生的剪切破坏单元较多，释放的剪切应变能也较多，后者可由前者得到解释；V形坑大且钝，这可由第3主应力与剪切带之间的夹角受扩容角控制得到解释。即使没有规定抗拉强度的软化规则，也有一定的拉伸应变能被释放出来，这是由于当设定的抗拉强度大于由残余的内聚力和内摩擦角决定的最大抗拉强度时，FLAC程序会自动按脆性规则降低设定的抗拉强度，这通过检查C++源代码获知。

（2）高扩容角时，事件集中出现的时间有所推迟，从与能量释放有关的量上一般能观察到巷道围岩发生V形坑式破坏的征兆，而从声发射数上一般并不能观察到，这说明了统计事件的能量比数目更具优越性。从节省计算时间的角度讲，将扩容角置零，会使模型的静力平衡容易达到，但是却可能低估破坏区的尺寸及释放的能量，也可能观察不到破坏的前兆，这一方面给研究带来了不便，另一方面会给岩石工程带来潜在的风险。

［1］Guenot A.Borehole breakout and stress fields［J］.International Journal of Rock Mechanics and Mining Sciences，1989，26：185-195.

［2］Ortlepp W D，Stacey T R.Rockburst mechanisms in tunnels and shafts［J］.Tunnelling and Underground Space Technology，1994，9（1）：59-65.

［3］陆家佑，王昌明.根据岩爆反分析岩体应力研究［J］.长江科学院学报，1994，11（3）：27-30.LU Jiayou，WANG Changming.Study on back analysis for stress of rock mass from information of rock bursts［J］.Journal of Yangtze River Scientific Research Institute，1994，11（3）：27-30.

［4］Ewy R T，Cook N GW.Deformation and fracture around cylindrical openings in rock——II.Initiation，growth and interaction［J］.International Journal of Rock Mechanics and Mining Sciences，1990，27：409-427.

［5］Martin CD.The effect of cohesion loss and stress path on brittle rock strength［J］.Canadian Geotechnical Journal，1997，34（5）：698-725.

［6］王学滨，潘一山，陶帅.不同尺寸的圆形隧洞剪切应变局部化过程模拟［J］.中国地质灾害与防治学报，2009，20（4）：101-108.WANG Xuebin，PAN Yishan，TAO Shuai.Numerical simulation of shear strain localization processes along circular tunnels with different diameters［J］.The Chinese Journal of Geological Hazard and Control，2009，20（4）：101-108.

［7］孟陆波，李天斌，周春宏，等.某水电站洞室岩爆的尺寸效应［J］.水力发电，2008，34（1）：29-31.MENG Lubo，LITianbin，ZHOU Chunhong，et al.Rock burst's dimensional effect on the tunnels of a hydropower station［J］.Water Power，2008，34（1）：29-31.

［8］王学滨，王玮，潘一山.孔隙压力条件下圆形巷道围岩的应变局部化数值模拟［J］.煤炭学报，2010，35（5）：723-728.WANG Xuebin，WANG Wei，PAN Yishan.Numerical simulation of the strain localization of the surrounding rock of a circular tunnel at different pore pressures［J］.Journal of China Coal Society，2010，35（5）：723-728.

［9］蒋承林，俞启香.煤和瓦斯突出机理的球壳失稳假说［J］.煤矿安全，1995（2）：17-25.JIANG Chenglin，YU Qixiang.A hypothesis of spherical shell losing stability about outburst mechanism［J］.Safety in Coal Mines，1995（2）：17-25.

［10］王学滨，陶帅，潘一山，等.基于环向应力岩爆判据的不同级别岩爆区分布的数值模拟研究［J］.四川大学学报（工程科学版），2011，43（5）：7-12.WANG Xuebin，TAO Shuai，PAN Yishan，etal.Numerical simulation of the distribution of rockburst regions with different magnitudes based on loop stress criteria［J］.Journal of Sichuan University（Engineering Science Edition），2011，43（5）：7-12.

［11］王学滨，潘一山，张智慧.扩容角对圆形巷道岩爆过程的影响［J］.中国工程科学，2010，12（2）：40-46.WANG Xuebin，PAN Yishan，ZHANG Zhihui.Numerical simulation of the rockburst processes of a circular tunnel at different dilation angles［J］.Engineering Sciences，2010，12（2）：40-46.

［12］郭保华.岩样尺度、孔道及端部摩擦效应的数值分析［J］.岩石力学与工程学报，2009，28（S2）：3391-3401.GUO Baohua.Numericalanalysis of size scale，inner hole and end restraint effects of rock samples［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（S2）：3391-3401.

［13］张培文，陈祖煜.剪胀角对地基承载力的影响［J］.岩土力学，2006，27（3）：463-465.ZHANG Peiwen，CHEN Zuyu.Significance of soil dilatancy in bearing capacity of ground［J］.Rock and Soil Mechanics，2006，27（3）：463-465.

［14］邓楚键，郑颖人，王凯，等.有关岩土材料剪胀的讨论［J］.岩土工程学报，2009，31（7）：1110-1114.DENG Chujian，ZHENG Yingren，WANG Kai，et al.Some discussion on the dilatancy of geotechnical materials［J］.Chinese Journal of Geotechnical Engineering，2009，31（7）：1110-1114.

［15］张春，茅献彪，倪海敏，等.岩石剪胀效应的数值模拟研究［J］.西安科技大学学报，2011，31（2）：181-187，229.ZHANG Chun，MAO Xianbiao，NIHaimin，et al.Numerical simulation of rock dilatancy effect［J］.Journal of Xi’an University of Science and Technology，2011，31（2）：181-187，229.

［16］赵金勇，蒋刚，杨磊，等.剪胀角对地基极限承载力影响的有限元极限分析［J］.南京工业大学学报（自然科学版），2008，30（5）：28-32.ZHAO Jinyong，JIANGGang，YANG Lei，etal.Effect of shear dilatancy angle on ultimate bearing capacity of foundation by FEM limit analysis［J］.Journal of Nanjing University of Technology（Natural Science Edition），2008，30（5）：28-32.

［17］王学滨.剪胀对岩样全部变形特征的影响［J］.四川大学学报（工程科学版），2005，37（5）：25-30.WANG Xuebin.Effects of shear dilatancy on entire deformational characteristics of rock specimen［J］.Journal of Sichuan University（Engineering Science Edition），2005，37（5）：25-30.

［18］Cherry J T，Schock R N，Sweet J.Theoretical model of dilatant behavior of a brittle rock［J］.Pure and Applied Geophysics，1975，113：183-196.

［19］GerbaultM，Poliakov A N B，Daignieres M.Prediction of faulting from the theories of elasticity and plasticity：what are the limits［J］.Journal of Structural Geology，1998，20：301-320.

［20］WANG Xuebin，MA Jin，LIU Liqiang.Numerical simulation of failed zone propagation process and anomalies related to the released energy during a compressive jog intersection［J］.Journal of Mechanics of Materials and Structures，2010，5（6）：1007-1022.

［21］王学滨.扩容角对初始随机材料缺陷岩石试样破坏过程及力学行为的影响［J］.中国有色金属学报，2007，17（7）：1063-1071.WANG Xuebin.Effects of dilation angle on failure process and mechanical behavior for rock specimen with random material imperfections［J］.The Chinese Journal of Nonferrous Metals，2007，17（7）：1063-1071.

［22］Vardoulakis I.Shear band inclination and shear modulus of sand in biaxial tests［J］.International Journal for Numerical and Analytical Methods in Geo mechanics，1980，4：103-119.