层状岩体边坡整体安全系数数值求解方法对比分析

汪良峰 吴小刚 徐鼎平,3,4 周玉新

(1.中钢集团马鞍山矿山研究院有限公司；2.金属矿山安全与健康国家重点实验室；3.中国科学院武汉岩土力学研究所；4.岩土力学与工程国家重点实验室)

岩土工程中的边坡稳定性评价是相对“古老”而又十分重要的研究课题。随着工程建设持续而快速地开展，越来越多的边坡稳定性问题受到学者们的广泛关注和研究。边坡稳定性分析需要研究2个主要问题：一是潜在主滑动面，二是整体安全系数。层状岩体作为人类工程活动最密切的地质体之一，在地壳浅部占比超过70%，因此，层状岩体边坡稳定性分析也同样面临这2个急需解决的问题。层状岩质边坡的一个突出特点是岩体的强度既受层面和层间岩石力学性质的协同控制，也受到边坡临空面与层面产状组合关系的制约。倘若忽略层状岩体的这一突出特点，而将其假定或等效为各向同性岩体进行处理，很可能错估边坡的潜在破坏模式，给出与实际相差较大的边坡整体安全系数，进而误导边坡工程的开挖、加固设计与失稳治理。

随着计算机技术的发展和弹塑性力学理论的逐步完善，各种反映层状岩体横观各向同性力学特征的力学模型和模拟手段被提出并应用于层状岩体边坡的稳定性分析中。就数值分析方法的经验积累和工程实用性而言，连续介质分析方法(有限元法、有限差分法等)仍是当前边坡稳定数值分析最主要的手段，因而采用它们进行层状岩体边坡稳定性数值分析方法也日益增多。赵尚毅等[1]采用低强度弹塑性薄层单元模拟岩体中的软弱夹层，采用接触面单元模拟层面，进而对层状岩体边坡采用强度折减法进行稳定性分析。林杭等[2]、张社荣等[3]采用低强度弹塑性薄层单元来模拟软硬互层岩体中的软弱夹层，然后采用强度折减法研究了软弱夹层产状、边坡角与边坡整体安全系数、潜在破坏模式间的关系。吴顺川等[4]则基于遍布节理模型，通过对薄层状岩质边坡层间岩石和层面折减相同系数进行了稳定性问题研究。总体而言，当前这些层状岩体边坡稳定性数值模拟方法是在客观认识层状岩体结构组成对其力学性质影响，将层状岩体视为层间岩石和层面的复合材料的前提下，采用不同模拟手段来反映其力学性质的层面依赖性特征及其对边坡稳定的影响。在这些研究的基础上，进而形成了4种层状岩体边坡稳定性的安全系数数值求解方法，即基于层面接触单元的强度折减法、基于层面薄层单元的强度折减法、基于隐式节理等效连续模型的强度折减法和基于层状岩体点安全系数的面积加权平均法。其中，基于层状岩体点安全系数的面积加权平均法是作者在层状岩体横观各向同性力学模型基础上推导得到其点安全度表达式之后，借鉴杨涛等[5]的各向同性岩体边坡整体安全系数面积加权法思想在本文中首次提出。这些研究对加深层状岩体边坡变形破坏的认识和促进层状边坡稳定定量评价的发展做出了重大贡献。但是，当前对层状边坡稳定性分析往往是孤立地采用上述某一方法具体开展，计算结果的合理性也多采用层状岩体各向同性化后的经典条分法计算结果来验证。事实上，边坡稳定性数值分析方法由于考虑了岩土体的本构关系及变形对计算结果的影响，它实质上是与基于刚体极限平衡理论的经典条分法并行的一种分析方法。而且，层状岩体边坡的经典条分法在固有的条间力假设情况下又附加了岩体力学性质各向同性这一假设。显然，对于同一分析对象，这两者计算结果间的相互校核不仅在理论上缺乏严密性，在工程上也缺乏说服力。因此，对现有顺层岩质边坡稳定安全系数数值求解方法进行系统、全面地横向对比研究，对于客观认知和准确评价层状岩体边坡稳定十分必要。

本文从这4种层状岩体边坡整体稳定性数值分析方法的基本原理出发，同时，基于相同的边界条件、计算收敛标准和相同的力学参数，采用这4种方法对已有文献中的经典算例进行精细分析对比验证，系统总结各种求解方法的优势、局限性和适用范围，以期为使用者求解层状岩体边坡整体稳定性提供借鉴和参考。

1 整体稳定性数值分析方法原理

1.1 方法一

(1)

(2)

1.2 方法二

基于层面薄层单元的强度折减法将层间岩石处理为各向同性岩体，层面采用薄层实体单元模拟，层间岩石和薄层均采用各向同性的理想弹塑性力学模型。求解安全系数时，同时对层间岩石和薄层抗剪强度参数折减相同的系数并计算至最终状态，不断试算直至计算不收敛，以此时的折减系数作为边坡整体安全系数，以此时等效塑性云图贯通所确定的滑动面为边坡潜在滑动面。层状岩体力学模型中，局部坐标系由层面倾向d、走向s和法向n构成；全局坐标系由正东向x,正北向y和竖直向z组成，α和β分别为走向角和倾向角。见图1。

图1 层状岩体力学模型坐标系

1.3 方法三

基于隐式节理等效连续模型的强度折减法在将层状岩体视为层间岩石和层面的复合材料的基础上，分别基于全局坐标系和局部坐标系(图1)建立基于线性Mohr-Coulomb准则和最大拉应力准则的层间岩石和层面力学模型，而后合成得到层状岩体的横观各向同性理想弹塑性模型，即Ubiquitous-Joint模型[6]；以此模型为基础，再分别对层间岩石和层面抗剪强度参数折减相同的系数，不断试算直至不收敛，以此时的折减系数作为边坡整体安全系数，以此时等效塑性云图贯通所确定的滑动面为边坡潜在滑动面。Ubiquitous-Joint模型本质上是考虑层状岩体强度各向异性的一类隐式节理等效连续模型，该模型通过特定的几何参数(层面倾向、倾角)和强度参数(层面黏聚力和摩擦角)来反映定向节理/层面的力学性质，因而无需在网格模型里给出特定的接触面单元或实体单元来模拟层面。

1.4 方法四

1.4.1 整体安全系数定义

基于层状岩体点安全系数的面积加权平均法同样是基于Ubiquitous-Joint模型，不同的是基于该模型在未进行强度参数折减的情况下进行弹塑性计算获得边坡体内的真实应力场，而后插值获得临界滑动面每一个网格质心上的应力分量，再基于其上的应力分量，求得每一个网格单元的点安全系数，最后基于点安全系数的面积加权法求得应力-应变分析获得的边坡整体安全系数，即

(3)

式中，Fsi为指定网格单元的点安全度;Ai为指定网格单元的点面积。

基于式(3)这一定义，要确定边坡的整体安全系数，势必要先求得各网格单元的点安全度。

1.4.2 点安全度

1.4.2.1 层间岩石点安全度

层间岩石应力状态由全局坐标系下的应力状态决定，其点安全度Fsm依据空间应力状态和Mohr-Coulomb剪切包络线的位置关系(图2)采用最小距离法[6]获得，即

Fsm=AC/AB=

[2Cm+(σ1+σ3)tanφm]cosφm/(σ1-σ3) ,

(4)

式中，cm为层间岩石黏聚力,MPa；φm为层间岩石内摩擦角,(°)；σ1和σ3分别为指定网格单元的最大主应力和最小主应力,MPa。

1.4.2.2 层面点安全度

层面应力状态由层面产状决定的局部坐标系上的应力状态所决定，其安全度Fsb基层面应力状态由层面产状决定的局部坐标系上的应力状态所决定。层面的抗剪安全度基于强度储备法定义，即

图2 全局坐标系下的应力状态与Mohr-Coulomb剪切强度包络线关系

Fsb=(cb+σbntanφb)/τb,

(5)

式中，cb为指定网格单元的层面黏聚力,MPa；φb为指定网格单元的层面内摩擦角,(°)；τb和σbn分别为指定网格单元层面上的剪应力和法向应力,MPa。

1.4.2.3 层状岩体点安全度

由于岩体的破坏总是沿着最薄弱路径破坏，因此，层状岩体点安全度Fs取Fsm和Fsb的最小值，即

Fs=min(Fsm,Fsb) .

(6)

2 典型算例对比研究及分析

基于FLAC3D[6]采用2个经典算例开展上述4种层状岩体边坡整体稳定性数值分析方法的对比研究。比较对象为潜在滑动面的形状、位置及整体安全系数。在强度折减法中，失稳判据为计算不收敛。

2.1 算例1

图3为文献[1]提供的一个处于临界状态的顺层岩质边坡几何模型。该边坡层理面间距为10 m，层面顺倾且倾角为40°,边坡角为65°。

图3 算例1几何模型

4种边坡整体稳定性数值分析方法的计算模型采用单元厚度模型，以等效进行平面应变分析。数值计算模型几何尺寸为水平方向125 m，竖直方向95 m,见图4。4种安全系数求解方法所采用的物理力学参数根据文献[1,7]给定并列于表1。

图4 算例1数值计算模型

方法容重/(kN/m3)层间岩石层面弹性模量层间岩石/GPa层面泊松比层间岩石层面黏聚力/MPa层间岩石层面内摩擦角/(°)层间岩石层面一25-10kn=10 000 kN/mks=3 846 kN/m0.2-1.00.123824二2517100.01 GPa0.20.31.00.123824三、四25-10-0.2-1.00.123824

注：kn为层面法向刚度，ks为层面剪切刚度。

图5为4种求解方法给出的整体安全系数和潜在滑动面(方法一以位移云图给出，方法二、三以等效塑性应变云图给出，方法四以点安全度云图给出)。可以看出，4种方法得到的整体安全系数大小存在差异，方法一、二和四与文献[1，7]给出的原解答1.05十分接近，方法三则与该解答有较大差异。此外，这4种方法给出的潜在滑动面非常接近。从滑动面也可以看出，该边坡表现为沿第3组层面滑动，其尾部后缘存在一定的拉破坏并往上延至坡顶，使得边坡的潜在滑动面为非圆弧形滑动。这与岩体各向同性处理后获得的圆弧形深层滑动模式(图6)有显著差别。

图5 算例1中4种求解方法的边坡整体安全系数及滑动面

图6 各向同性处理后的算例1边坡滑动面

2.2 算例2

文献[3]基于Sarma 极限平衡法和有限元强度折减法探讨了多个复杂多层软弱夹层边坡岩体的破坏机制及稳定性特征。选取该文中的层面倾角和坡角都为45°、层面间距为30 m的顺层岩质边坡(图7)为研究对象来进一步验证这4种安全系数数值求解方法给出的边坡整体安全系数和潜在滑动面的差异。仍采用单元厚度的模型，计算网格模型几何尺寸为水平方向500 m，竖直方向200 m，见图8。4种方法所采用的物理力学参数根据文献[3]给定并汇总于表2。

图7 算例2几何模型

图8 算例2数值计算模型

方法容重/(kN/m3)层间岩石层面弹性模量层间岩石/GPa层面泊松比层间岩石层面黏聚力/MPa层间岩石层面内摩擦角/(°)层间岩石层面一26.5-2kn=555 000 kN/mks=15 000 kN/m0.2-0.30.053019.3二26.51020.26 GPa0.20.450.30.053019.3三、四26.5-2-0.2-0.30.053019.3

注：kn为层面法向刚度，ks为层面剪切刚度。

图9为4种求解方法给出的整体安全系数和潜在滑动面(方法一以位移云图给出，方法二、三以等效塑性应变云图给出，方法四以点安全度云图给出)。可以看出，这4种方法得到的整体安全系数大小存在差异，方法一、二和四与原解答1.22十分接近，方法三则与该解答有较大差异。此外，这4种方法给出的潜在滑动面非常接近。方法一、二给出的滑动面表明，该顺层边坡表现为沿第一个层面滑动，而后在下部穿切层间岩石，从坡脚滑出的“顺层-圆弧”复合滑动模式。基于隐式节理等效连续模型的强度折减法也完全反映了这一复合破坏模式。这与各向同性处理后的边坡圆弧形滑动(图10)差异明显 。

图9 算例2中4种求解方法的边坡整体安全系数及滑动面

图10 各向同性处理后的算例2边坡滑动面

综上所述，4种层状岩体边坡安全系数的数值求解方法给出的结果都能满足实际需求，但是每一种方法都存在固有缺陷。例如基于层面接触单元的强度折减法和基于层面薄层单元的强度折减法由于要将层面在网格模型中直接用几何形式表现出来，故只能针对层面间距较大、边坡岩体结构相对比较简单的顺层岩质边坡，否则数值求解比较复杂的以薄层岩体(图11)为主的顺层岩质边坡安全系数时，网格模型建立比较困难。

图11 两河口水电站厂址区域陡倾薄层岩体

现有算例也表明，基于隐式节理等效连续模型的强度折减法求解的层状岩体边坡稳定的整体安全系数与其他方法差异较大，这很可能是由于在层状岩体破坏过程中，层面强度与层间岩石强度的调动与发挥过程存在明显不同，在进行强度折减时同时除以相同的折减系数，忽略了这个事实。方法四相对而言适用性更广，理论也相对严密。因为它是基于未折减强度参数获得的应力场，对于岩体中的层面也不要求在网格模型中予以显性表达，因此，求解层状岩体边坡不存在理论上的明显缺陷。但是该方法使用前需要确知潜在滑动面的位置和形状。算例1和2之所以能采用该方法求得安全系数，是因为这2个边坡都已接近临界状态，基于未折减强度参数给出的位移场或等效塑性应变场已经能较为明确地给出潜在滑动面。因此，如何在基于未折减强度参数的应力场条件下获得非临界状态下的边坡潜在滑动面是该方法求解安全系数的一个关键问题。上述4种安全系数数值求解方法的优劣对比见表3。

3 结 论

(1)层状岩体是一类典型的横观各向同性材料，其力学性质由层面和层间岩石协同控制，并受临空面与层面产状组合关系的制约，因此,其边坡潜在滑动面与均质各向同性边坡的圆弧形滑动面差异较大。在对其进行边坡稳定性评价时，不能笼统地以各向同性材料对待，即无论是采用何种力学模型或者模拟手段，均应通过不同的方式表现这种材料力学性质的层面依赖性特征。

表3 层状岩体边坡4种整体安全性数值分析方法主要特点

(2)将层状岩体中的层面处理为实体单元和将层面处理为接触面单元，再基于各向同性力学模型的强度折减法可以获得相对接近的安全系数计算结果，但这类方法要将层面在网格模型中予以显性表示，因此，只能针对层面间距较大、岩体结构相对比较简单的层状边坡，对于比较复杂的层状边坡(薄层岩体边坡等)，因网格模型建立比较困难则难以应用。基于隐式节理等效连续模型的强度折减法可以采用均一化网格处理复杂层状边坡稳定性分析，但是其获得的整体安全系数计算结果与其他方法差别较大，表明层状岩体在破坏过程中，其层间岩石和层面的强度调动过程和机理存在较大差异，求解安全系数时对这二者折减相同系数的做法有待商榷，尚需进一步研究。

(3)在顺层岩质边坡潜在滑动面上，基于接触面单元的强度折减法、基于层面薄层单元的强度折减法、基于隐式节理等效连续模型的强度折减法均能给出比较一致的形状。层状岩体点安全系数的面积加权平均法则取决于边坡本身的状态，若处于临界状态附近，则破坏模式与基于隐式节理等效连续模型的强度折减法结果相同；若边坡处于非临界状态，则需借助其他方式来确定潜在滑动面。