基于矢量和法安全系数的边坡稳定性分析

2013-11-05钟冬望刘建程黄小武

武汉科技大学学报 2013年6期

何理，钟冬望，刘建程，黄小武

（1.武汉科技大学资源与环境工程学院，湖北武汉，430081；2.武汉科技大学理学院，湖北武汉，430065；3.武汉科技大学冶金工业过程系统科学湖北省重点实验室，湖北武汉，430065）

边坡稳定性分析中最重要的概念就是安全系数，但对不同的工况或不同的方法所使用安全系数的定义是不同的。目前求解边坡安全系数最常用的方法主要有极限平衡法和有限元强度折减法。极限平衡法虽经过80多年的发展得到充分完善［1］，但该法只考虑静力平衡条件和 Mohr－Coulmb准则，主要有瑞典法、Bishop法、Morgenstern－Price法、Spencer法和 Sarma法等［2－3］。有限元强度折减法是基于边坡有限元变形分析和强度储备安全系数的稳定性分析方法，其对应的安全系数求解方法有二：一是对岩体的强度参数黏聚力c和内摩擦角φ按同一折减系数F进行折减，此时的折减系数F即为边坡的稳定安全系数；二是安全系数定义为沿滑面各点抗滑力与下滑力积分代数和的比值。极限平衡法与有限元强度折减法均是基于强度储备概念上的稳定性分析方法，但两种方法计算内力的方式不同，相比之下，有限元强度折减法在用于内力计算时较为准确。

通常根据安全系数的大小来评价边坡安全与否，这就涉及安全系数的定义问题。上述所提到的安全系数的定义存在如下问题［4］：①安全系数定义的物理意义不明确；②岩土体强度参数黏聚力c和内摩擦角φ按同一系数F进行折减不合理，如郑宏［5－6］、郭明伟［7］等通过计算表明，在对强度参数进行折减的同时，内摩擦角φ与泊松比μ的取值应满足sinφ≥1－2μ；③计算需要反复迭代，且计算量大。针对安全系数物理意义不明确、人为假定过多、不能真实反映边坡应力状态等问题，葛修润［8］、郭明伟［9］等提出边坡矢量和法，并考虑力的矢量特征，得到抗滑稳定性安全系数应以力的矢量比来定义。但上述两种方法安全系数定义的物理或力学意义模糊的根本原因就在于没有充分考虑滑动的矢量特征。为此，本文针对极限平衡法与有限元强度折减法计算边坡安全系数时存在的问题，采用矢量和法安全系数对某铁矿采场进行边坡稳定性分析，将其计算结果与传统分析法计算结果进行比较，并验证了矢量和法计算安全系数的合理性与适用性，以期为采用矢量和法安全系数用于计算边坡稳定性提供依据。

1 矢量和法安全系数

安全系数最初的定义为：边坡潜在滑动面上总抗滑力∑R与总滑动力∑T的比值，即：

边坡滑面上抗滑力与滑动力均为矢量，故式（1）中对抗滑力与滑动力的求和均应为矢量和。本文应用“矢量和法安全系数”是建立在力的矢量基础上的。

1.1 基本假定

假定条件如下：①若研究对象为二维问题，边坡滑动面l已知，则构成滑动区域为S；②若边坡荷载及岩土体基本物理力学参数已通过勘探得到，则基于有限元计算可得到其应力分布情况；③滑动面上岩土体强度特性服从Mohr－Coulomb屈服准则；④矢量和法安全系数F（θ）定义为：沿计算方向θ，滑动面上提供抗滑力的各力沿此方向投影代数和∑R（θ）与提供滑动力的各力沿此方向投影代数和∑T（θ）的比值：

式中：τf为抗剪强度，Pa；c 为黏聚力，Pa；σ 为法向应力，Pa；φ 为内摩擦角，（°）。

1.2 安全系数的计算

图1为边坡安全系数计算方向示意图。由图1可看出，在边坡荷载、边界条件、滑动面位置、岩土体的基本物理力学参数已知的情况下，采用有限元法计算可得到滑动面上任意一点i在局部坐标系x′o′y′下的应力为σi和τi，滑动面上过点i的切线与总体坐标系x轴正向的夹角为αi。

应力的正负号按岩土力学的规定，角度的正负规定为：从x轴正向出发，沿逆时针方向的角度为正，沿顺时针方向的角度为负，故图1中的αi为负。

图1 边坡安全系数计算方向示意图Fig.1 Sketch for calculating direction of slope safety coefficient

根据摩擦理论，边坡潜在滑动面上任意一点i的滑动方向应沿着该点滑动面的切向方向，由修正的粘着理论可知，该点应力在滑裂面切向方向的分量τi与Δli的乘积为该点处静滑动摩擦力的大小，该点的滑动趋势方向与其静滑动摩擦力方向相反。对于整个滑体而言，其静滑动摩擦力就是滑动面各点处静滑动摩擦力的合力，合力方向的反方向为整个滑体的潜在滑动趋势方向。而矢量和法安全系数的计算方向就是滑体的整体滑动趋势方向，也就是滑动面各点处静滑动摩擦力合力方向的反方向，将每点静摩擦力的合力分别投影到x和y轴后的比值，即得计算方向与水平线的夹角为

式中：Fxi＝τiΔlicosαi、Fyi＝τiΔlisinαi分别为τiΔli在x轴与y轴上的投影。

1.3 矢量和法安全系数的计算

图2为抗滑稳定性分析时安全系数求解示意图。由图2可看出，以滑动面为研究对象，滑面上一点处的岩土体黏聚力为ci，内摩擦角为φi，滑体作用于基岩上的荷载为正应力σi和剪应力τi，基岩作用在滑体上的反力分别为σ′i和τ′i，它们是一对大小相等、方向相反的作用力与反作用力。矢量和法安全系数表达式的计算如下：

图2 安全系数求解示意图Fig.2 Sketch for calculating safety coefficient

滑动面上的滑动力由自重及外荷载引起。在微弧段Δli内，将引起滑动力的σi和τi分别沿安全系数计算方向的反方向的投影有：

滑动面上抗滑力由岩土体的黏聚力、摩擦力和基岩对滑动面的法向反力提供。在微弧段Δli内，黏聚力与摩擦力之和为滑动面上岩土体的抗剪强度，其值由摩尔－库仑强度准则计算，基岩对滑动面的法向反力为σ′i（＝σi）。将抗滑力分别沿安全系数的计算方向投影，即：

1.4 矢量和法安全系数求解步骤

求解步骤如下：①根据地质勘察资料建立有限元数值计算模型，在模型中反映地质勘察的滑面；②对模型进行常规有限元弹（塑）性分析计算，得到各点的应力状态；③根据有限元弹（塑）性计算结果，计算滑面上各处剪力的矢量和，该剪力矢量和的方向即为安全系数的计算方向。通常按有限元网格划分节点，将滑面分为若干个小段进行计算；④安全系数的计算先按式（7）和式（11）计算滑面上各处滑动力与抗滑力沿此方向投影的代数和，再通过式（12）计算就可得边坡整体的安全系数。

2 工程应用

2.1 边坡概况

某铁矿采场由狮子山、尖林山、象鼻山3个矿体组成。采场总体走向为NW300°，总体长度约为2200m，坑底总面积为8150m2。东露天采场地质结构复杂，出现断层、节理发育，边坡稳定性差。Ⅵ区位于区域断层F9上盘，其西侧与Ⅴ区有区域性节理J2，整体坡度约为48°。区内有两组区域性节理J3、J4，这两组主节理将该区划分为上下两个部分，上部形似锲形体，下部形似正三角形体。另外还有两组节理发育，一组节理近于平行，产状320°∠75°，密集发育；另一组节理顺坡面发育，两组节理将岩体切割成碎裂结构。Ⅵ区岩体结构与Ⅲ－Ⅲ′剖面工程地质图分别如图3、图4所示。

图3 Ⅵ区岩体结构Fig.3 Rock structure of RegionⅥ

图4 Ⅲ－Ⅲ′剖面工程地质图Fig.4 Engineering geological map of SectionⅢ－Ⅲ′

2.2 岩体力学参数取值

本文参考中国科学院岩土力学研究所1995年编写的《大冶铁矿狮子山北邦A区尖F9断层上盘－60～－96m边坡稳定性和滑坡防治研究》报告提供的岩体参数，如表1所示。

表1 岩体力学参数Table1 Rock mass mechanical parameters

2.3 计算结果分析

2.3.1 极限平衡法

基于极限平衡原理，采用边坡工程中广泛运用的边坡稳定性计算软件Slope／w，包括 Morgenster－Price法、瑞典条分法、Bishop法和Janbu法等对边坡稳定性进行计算，同时考虑孔隙水压力、程序自动搜索最危险滑动面，并给出最小安全系数，如图5所示（边坡岩性见图4）。

由图5可看出，由于上述计算方法均是基于极限平衡原理，故计算得到最危险滑动面基本重合，但各方法严格程度及考虑因素不尽相同，以致最小安全性系数各不相同。极限平衡法研究边坡稳定性的合理性与适应性已得到工程验证。4种情况下边坡安全性系数如表2所示。

图5 基于极限平衡法的边坡最危险滑动面Fig.5 The most dangerous sliding surface of slope based on limit equilibrium method

表2 边坡安全性系数Table2 Slope safety coefficients

2.3.2 有限元强度折减法

借助有限元分析软件Ansys模拟边坡应力－应变关系，边坡按照平面应变问题处理，单元类型选用PLANE42平面单元，岩土体材料本构模型采用非线性弹塑性德鲁克－普拉格（D－P）模型。DP模型不仅计入静水压力作用，而且还考虑中间主应力的影响，克服了 Mohr－Coulomb模型的主要弱点，被广泛应用于国内外岩体力学计算与数值分析中［10］。

Ⅲ－Ⅲ′剖面的实际坡高为431.12m，宽为566.5m。为了消除计算时边界条件对坡体应力－应变求解的影响，在进行求解时将模型尺寸取高为700m，宽为1500m。选取位移边界条件，左右边界位移x方向固定，y方向自由，底边x、y方向固定。模型经有限元网格划分，共得到节点4132个，单元4007个。边坡有限元计算模型如图6所示。

图6 边坡有限元模型Fig.6 Finite element model of the slope

对岩体的强度参数黏聚力c和内摩擦角φ按同一折减系数F进行折减，利用Ansys默认的计算收敛准则（力和位移不收敛），得到有限元强度折减安全系数为1.26。极限状态下边坡应力、应变分布如图7所示。由图7可看出，当岩体结构面强度参数黏聚力c和内摩擦角φ同时折减1.26倍后，最大剪应力区域集中出现在边坡Ⅵ区Ⅲ－Ⅲ′剖面断层区域，此时xy塑性剪应变在断层及节理面上已基本贯通，形成塑性带，边坡处于极限平衡状态，坡体具有失稳下滑倾向。

2.3.3 矢量和法

图7 边坡剪应力、剪应变云图Fig.7 Shear stress and shear strain nephogram of the slope

选用由Slope／W搜索出的最危险滑动面作为临界滑动面，滑面为圆弧形（见图5），利用Ansys用户可编程特性APDL语言的二次开发功能，自定义算法，根据上述矢量和法安全系数的定义与求解步骤编写计算程序，得到该边坡给定滑动面上各处剪力的矢量和方向，即安全系数的计算方向为32.84°，安全系数F（θ）值为1.34。

3 结果对比分析

采用不同方法得到边坡安全系数计算结果的比较如表3所示。由表3可看出，矢量和法计算结果与极限平衡法和有限元强度折减法之间最大相对误差为9.7%，误差范围为4.2%～9.7%，由此表明矢量和法计算安全系数的适用性。极限平衡法主要是基于重度增加方法来计算边坡安全系数，随着边坡体的重度增加，正应力的增大比例较大于剪应力的增大比例，故边坡滑体抗剪强度得到加强，从而导致极限平衡法中的Morgenster－Price法、瑞典条分法、Bishop法较矢量和法计算结果大。

表3 边坡安全系数计算结果比较Table3 Comparison of calculated slope safety coefficient

矢量和法同有限元强度折减法计算结果相对误差较小，仅为5.9%，主要由于本文选取圆弧滑面为最危险滑动面，按照强度安全系数的定义是滑动面上各处的抗滑力对滑弧圆心的力矩与各处滑动力对滑弧圆心的力矩之比，如同矢量和法安全系数，具有明确的物理意义。

有限元强度折减法计算结果与极限平衡法安全系数相对误差较大，分析其原因主要在于：①有限元计算的数值收敛性受多种因素的影响，由此得到的安全系数的合理性与惟一性受到质疑；②抗剪强度с和内摩擦角φ在边坡失稳时发挥程度以及衰减速度与程度并非完全相等，两者具有各自的安全储备。二者采用同一折减系数会造成塑性区失真；③强度折减后，岩体残余强度考虑不足。工程实践表明，自重情况下，某铁矿Ⅵ区边坡整体稳定性良好，表明采用矢量和法安全系数用以衡量边坡稳定性切实可行。