用DDA方法验证倾倒边坡变形的制动机制

何传永，孙 平，吴永平，段庆伟

（中国水利水电科学研究院 岩土工程研究所，北京 100048）

1 研究背景

岩质边坡的失稳模式主要包括平面滑动、楔形体和倾倒破坏。工程界对边坡倾倒破坏的认识比较晚，进行数值分析也最困难。

倾倒边坡的特征是有一组反倾坡向的岩石层面或高连通率的结构面，破坏形式主要表现为倾倒的岩柱旋转弯曲折断破坏。从坡面上观察，岩柱呈台阶状并伴有裂缝；在地质平洞中，能观察到成组的张性裂隙，岩柱不连续弯曲和折断，在倾倒体与基岩的交界面上有层间错动现象。这些表象说明，边坡曾产生过较大的变形后稳定了下来，Hoek认为倾倒边坡存在制动机制[1]。

倾倒边坡是非连续性大变形问题，连续介质小变形的有限元方法显然不适合求解倾倒边坡问题，刚体极限平衡方法仅适合最简单的倾倒模式，而建立在非连续介质模型基础上应力应变分析的DDA（Discontinuous Deformation Analysis）方法[2]可以模拟岩块的平移、错动、旋转等各种形态的变形，能较真实反映倾倒边坡变形的实际情况，可能是分析岩质边坡倾倒破坏的有效手段[3-5]。非连续变形分析方法DDA是由石根华博士提出的分析块体运动和变形的一种数值方法。它以自然存在的节理面或断层切割岩体形成块体单元，以各块体的位移作为未知量，通过块体间的接触和几何约束形成一个块体系统。在块体运动的过程中，严格满足块体间不侵入条件，将边界条件和接触条件等一同施加到总体平衡方程。总体平衡方程式是由系统的最小势能原理求得，求解方程式组即可得到块体当前时步的位移场、应力场及块体间的相互作用力。反复形成和求解总体平衡方程式，即可得到运动过程中各块体的相对位置及其接触关系，求得各个时步后直至最终平衡时的位移场、应力场等情况。因此，DDA方法适合模拟倾倒边坡岩石块体的移动、旋转及块体间张开、闭合等运动形式和运动过程，并据此判断岩体的破坏程度和破坏范围，根据指定滑面上力的平衡条件，计算出边坡安全系数，从而对倾倒边坡岩体的局部和整体稳定性做出正确的评价。非连续变形分析方法的特有优点是：完全的运动学及其数值可靠性；完全的一阶位移近似严格的平衡要求；正确的能量守恒和高计算效率[6]。因此，DDA方法计算结果更能接近事物的本质，更接近工程实际中的力学规律和破坏现象。

2 倾倒变形的经典例题

Goodman等[7]提出的边坡倾倒稳定分析方法是将边坡用反倾向的结构面切割成n块宽度为ΔX的矩形条块。由于条块所处位置不同，所受作用力不同，可能使条块处于3种状态：（1）稳定；（2）倾倒破坏；（3）滑动。因此将边坡岩体分为稳定区、倾倒区、滑动区3部分。

图1是Hoek书中引用的倾倒边坡的典型结构特征图[1]，最上面2条块为稳定区，最下面2条块为滑动区，中间3个条块是倾倒区。

图1 倾倒边坡的典型结构特征

已知条块右侧作用力的合力，可以根据力的平衡求得左侧合力。对于倾倒块体，将各作用力对条块左下端点取矩：

式中：Pn-1、Pn分别为条块右侧和左侧界面上作用力合力；Wn为岩块重；Yn、Ln、Mn分别为岩块高度、岩块的左、右侧有效接触高度；ΔX为岩块宽度；α为节理倾角；φ为岩块侧面摩擦角。

对滑动块体通过静力平衡可以得：

图2为Hoek和Bray提供的一个理想的实例。这个开挖边坡高92.5m，坡角为56.6°，岩体为倾角60°倾向山体的层状岩体，由16个岩块组成，坐落在阶梯状底面上。几个常数是：该梯状底面为每5个岩块上升1m（角β-α=5.8°）a1=5.0m，a2=5.2m，b=1.0m，ΔX=10.0m，坡顶面仰角为4°（符号参见图1）。

图2 Hoek和Bray经典算例

取岩块底面和侧面的摩擦系数为0.7855，按式（1）和式（2）分别计算各块倾倒力和滑动力，取Pn-1的较大者。从岩块4以上一直到岩块13倾倒力的数值都比较大，而从岩块3以后滑动力的数值比较大。所以，岩块4至13构成了潜在的倾倒区，而岩块1至3构成了滑动区；岩块14、15及16构成稳定区。阻止岩块1滑动所需要的抗滑力趋于零，这说明该边坡非常接近极限平衡状态。

很显然，这是一个精心设计且恰好接近极限平衡状态的倾倒边坡，当在坡脚岩块1处施加0.5kN支撑力T时，边坡处于稳定状态；反之，则处于不稳定状态。

倾倒例题计算是使用作者在石根华DDA程序的基础上根据工程需要所作的二次开发程序，重新编写了图形处理部分，将计算和前后处理统一在一个程序中，具有Windows界面。计算条件为：岩石弹性模量为2GPa；泊松比0.2；岩块容重γ=25kN/m3；岩块底面和侧面的摩擦系数为0.7855，不考虑凝聚力；块体接触刚度20GPa。

3 用DDA方法计算倾倒变形经典例题

已有学者张国新等[8]和孙东亚等[9]分别用流形元方法和DDA方法对Hoek&Bray的经典倾倒算例进行过计算模拟，取得很好的效果[8-9]，本文则试图通过DDA方法对Hoek&Bray的经典倾倒算例进行计算模拟，探讨倾倒边坡的的制动机制。从DDA计算结果的动画显示可以观察到，在上部倾倒力的作用下，下部的1—3号块首先产生滑动，其后的4—13号岩柱绕底部产生逆时针旋转。随着滑动区块体位移的增加，倾倒区块体的旋转角度也随之增加。由于倾倒区岩块底部旋转时是角对角接触，不容易保持平衡，同时由于底接触面右边台阶限制了岩柱的旋转空间，岩柱在旋转的同时，向下产生错动，这与在实际倾倒边坡的地质平硐中观察到的倾倒体与基岩的交界面上的层间错动现象是吻合的。在变形的过程中，随着岩柱绕底部产生逆时针旋转，岩柱间产生张开裂缝，岩柱间面与面接触转变为棱与面接触，如图3。

图3 计算中间结果，岩柱棱与面接触

图4 摩擦角φ=38.15°，倾倒变形结束后位移云图

随着滑动区1—3号块体向下滑动，倾倒区4—13号岩柱有了更大的逆时针旋转的空间，当岩柱间由棱与面接触逐步重新恢复到面与面接触时，倾倒变形最终停止，边坡重新回到稳定状态。图4是迭代计算结束后位移总量云图，从图中可以看到，最大位移发生在10号岩柱的顶部，总位移达14.75m。

为了分析倾倒位移的产生，发展和收敛过程，在10号块体顶端设置一个测点，在迭代过程中记录测点的水平位移、竖向位移和总位移。从图5测点位移过程线中可以看到，测点位移开始比较缓慢，然后加速，最后趋缓并最终达到平衡状态。

从DDA程序计算的过程中各迭代步的倾倒变形全过程，可以看到，1—3号块体向下滑动，4—13号岩柱倾倒旋转，14—16号岩柱是稳定的，这与理论解的滑动区、倾倒区和稳定区完全一致。

图5 10号岩柱顶端测点位移随迭代步的变化过程线

4 倾倒变形的制动机制探讨

在平面滑动破坏模式中，一旦发生较大位移，会引起整个边坡产生破坏，而对倾倒边坡来说，常有较大的变形或位移发生，但并没有产生严重的破坏后果。从图4中可以看到，计算结束时，尽管滑动块下部仍有位移空间，但摩擦力克服了动力的作用，变形最终还是停了下来，此时安全系数大于1.0，说明倾倒边坡在发生较大变形后，内在的制动机制克服了动力的作用，从而阻止了继续变形，不像平面滑动边坡那样一滑到底。对于这种现象，Hoek认为，处于极限平衡状态的倾倒边坡是一种亚稳定边坡，一旦岩柱产生旋转，岩柱倾角发生变化，随之阻止其进一步翻转所需要的摩擦力就会增加。当岩柱翻转到2（β-α）时，就将使岩块侧面的棱与面接触转换为连续的面接触，同时阻止岩柱进一步翻转所需之摩擦角也将急剧减小，甚至可能小到比初始平衡所需要的摩擦角还小，倾倒变形会随之停止，即倾倒变形的制动机制发挥作用，倾倒边坡由失稳状态重新回到稳态[1]。

从DDA验算Hoek例题的变形运动计算过程中，可以清楚的看到先是倾倒体底部块体产生滑动，中部岩柱产生旋转，块体之间产生裂缝并逐渐增大，相互间由面与面接触变成为棱与面接触；随着滑动区块体下滑，变形岩柱翻转角度的增加，岩柱间裂缝又逐渐变小，重新恢复面与面接触时，最终块体运动停止，倾倒边坡重新回到稳定状态。DDA计算模拟倾倒变形的过程验证了Hoek关于倾倒变形制动机制的论述。

DDA验算倾倒变形最终停止时，岩柱的倾角为42.5°，与初始的倾角60°倾角相比，岩柱翻转角为17.5°，而Hoek认为倾倒变形停止时岩柱翻转角为2（β-α）=11.6°，两者有5.9°的差别，DDA计算结果的岩柱翻转角相对大一些。分析产生差别的原因，应该是Hoek的理论计算是静态的，没有考虑块体的动力加速度，而DDA方法是模拟块体运动的过程，在迭代计算中，当前时步保持前一时步岩块的速度，总刚矩阵中加入了重要的惯性力项，因此可以认为，DDA计算与Hoek理论解的岩柱翻转角有5.9°的误差是动力计算与静力计算所产生的差异，倾倒变形停止时DDA计算的岩柱翻转角略大也是合理的。

在许多实际倾倒变形边坡中，可以观察到大量的表面位移和张裂缝的形成，在勘探平洞中可以看到倾倒变形体底部与基岩有较大位移而产生的层间错动，这些现象与DDA数值模拟的结果是吻合的。尽管倾倒边坡产生了较大的位移，但并没有像平面滑动边坡那样会瞬间失稳，释放巨大的动能，岩体高速运动并造成很大的破坏，而是随着岩柱倾角变化，阻止其进一步翻转所需要的摩擦力增加，在变形过程中自动停了下来，重新回到稳态，边坡仍然是稳定的，岩体分离量也比较有限，这是倾倒变形过程中制动机制作用的结果。

倾倒边坡能够产生较大的变形，但变形量并不能作为倾倒边坡可能失稳的依据，边坡安全系数的选择取决于允许边坡有多大的变形量。

5 结论

通过对Hoek经典倾倒例题的验算表明，DDA能够模拟倾倒边坡下部滑动区岩块的滑动，倾倒区岩柱旋转、错动，直到倾倒边坡重新回到稳态的动态过程，与理论方法的滑动区、倾倒区和稳定区的计算结果完全一致。

文中对倾倒边坡的制动机制进行了初步探讨，在DDA程序开始计算时，倾倒区岩柱绕底部产生旋转，岩柱间产生张开裂缝，随着岩柱的继续翻转，岩柱间产生裂缝又逐渐变小，当岩块侧面的棱与面接触重新转换为连续的面与面接触时，倾倒变形会随之停止，倾倒边坡由失稳状态重新回到稳定状态。验证了Hoek关于倾倒区岩柱旋转使阻止岩块进一步翻转所需之摩擦角急剧减小，是倾倒边坡变形的重要制动机制的论述。

[1]Hoek E，Bray J W.岩石边坡工程[M].卢世宗译，北京：冶金工业出版社，1981：183-193.

[2]Shi，G-H.Discontinuous Deformation Analysis：A New Numerical Model for the Statics and Dynamics of Block Systems[D].Ph.D.Thesis，Department of Civil Engineering，University of California，Berkeley，1988：378.

[3]Ishida T C，Chigira M，Hibino S.Application of the distinct element method for analysis of toppling observed on a fissured slop[J].Rock Mechanics and Rock Engineering，1987，20（4）：277-283.

[4]Jiang Y，Esaki T，Nagatomi M.Studies on toppling failure mechanism of slope in discontinuous rock mass[M].Mechanics of Jointed and Faulted Rock.Rotterdam：Balkema，1995：605-610.

[5]Prichard M A，Savigny K W.Numerical modelling of toppling[J].Canada Geotech，1991，27：823-833.

[6]石根华.数值流形方法与非连续变形分析[M].裴觉民译，北京：清华大学出版社，1997：92-93.

[7]Goodman R E，Bray J W.Toppling of rock slopes.Proceedings of Specialty Conference on Rock Engineering for Foundations and Slopes ASCE/Boulder[J].Colorado，ASCE，1976，2：201-234.

[8]张国新，赵妍，石根华，等.模拟岩石边坡倾倒破坏的数值流形法[J].岩土工程学报，2007（6）：11-16.

[9]孙东亚，彭一江，王兴珍.DDA数值方法在岩质边坡倾倒破坏分析中的应用[J].岩石力学与工程学报，2002，21（1）：39-42.