基于广义Hoek-Brown准则的多参数强度折减法在边坡稳定性分析中的应用

刘昱君， 刘明扬， 杜文杰

(1.中建国际投资(湖北)有限公司， 武汉 430071； 2.中国建筑国际集团有限公司， 香港 999077； 3.中国科学院 武汉岩土力学研究所,岩土力学与工程国家重点实验室， 武汉 430071； 4.中国科学院大学， 北京 100049)

一直以来，道路交通、水利水电等工程在国家发展战略中占有重要地位。中国中西部多山，地质条件复杂多变，工程施工过程中需要对山体进行开挖，这就对原山体造成了一定扰动，岩体稳定性因此降低。很容易在外力作用下(二次开挖、极端恶劣天气等)发生山体滑坡、泥石流等自然灾害。保障这些工程的安全重点在于保障道路、铁路沿线开挖山体以及水利水电工程附近山体的稳定，这些都可以归结为边坡稳定性问题[1-3]。

近年来，随着岩土体的本构模型和有限元理论的发展，有限元强度折减法已被广泛应用于边坡稳定性分析和其他相关工程问题上，成为边坡稳定性分析应用最广的方法之一。强度折减法的发展得益于子技术的发展：失效准则已不再局限于单一的计算不收敛准则，目前已经发展到塑性区贯通、塑性应变能突变、最大水平和竖直位移突变等多种失稳判据可以选用。同时，三维有限元模型强度折减法得到了很大的发展，这对于解决岩土工程中很多复杂边坡问题、提高计算精度是一个全新的思路。不过，有限元强度折减法在非连续介质，比如支挡结构与岩土介质共同作用方面，计算精度的问题依然有待提高。

1999年，Griffith首次采用有限元强度折减法进行边坡稳定性分析，求解得到的安全系数与传统极限平衡分析法得到的安全系数在数值上十分接近，表明用有限元强度折减法进行边坡稳定性分析是可以有效指导实践的[4-5]。当然，随着有限元的发展和广泛的分析应用，现有强度折减法的问题也暴露了出来，仅仅对强度参数黏聚力c、内摩擦角φ进行折减的方案忽略了边坡强度参数和变形参数之间存在的对应关系；另外对抗剪强度参数黏聚力c、内摩擦角φ进行折减模拟得到的结果往往与实际边坡破坏相差较大，如果以塑性区贯通为失效准则，模型底部出现一部分不合理的塑性区，甚至有时会先于破坏面贯通，另外实际边坡顶部的拉应力区也无法在模型中体现出来，这与实际边坡破坏是有差别的。于是，在抗剪强度折减法的基础上有学者相继提出了对其他参数，主要是变形参数进行折减的方案，这些学者的研究成果进一步推动了强度折减法的完善和发展，更好地与工程实际相结合并正确指导工程实际[6-10]。

郑宏等[11]指出如果只考虑对强度参数的折减，塑性区将集中在边坡深部，模型底部在约束的影响下发生低应力破坏而出现塑性区，这与实际工程情况不符；袁维等[12]证明了弹性模量和泊松比对安全系数和塑性区分布有一定影响，对模型中下部的剪切塑性区影响较大，抗拉强度对边坡的安全系数和塑性区分布的影响较大，因此有必要对变形参数进行折减；施建勇等[13]指出，只对抗剪强度参数进行折减并没有从土体的本构关系上全面反映强度折减法的内涵，对此提出了在抗剪强度参数折减的基础上，同时调整初始切线模量以满足土体的变形性质。

现针对有限元强度折减法研究成果及其存在的不足，基于Hoek-Brown准则，考虑岩体弹性模量的沿深度递增关系，提出一种针对弹性模量进行折减的方案；并结合边坡算例进行稳定性计算验证本文方案的合理性。

1 基于Hoek-Brown准则的折减方法探讨

1.1 广义Hoek-Brown准则

20世纪80年代，Hoek和Brown基于岩体性状理论，结合岩体工程分类理论，统计并分析了大量现场试验资料和三轴试验资料，总结出一个岩体破坏的半经验公式，Hoek-Brown准则[14-15]其极限主应力表达式为

(1)

(2)

(3)

(4)

式中：σ1为岩体破坏时的最大主应力；σ3为作用的最小主应力；σc为完整岩石单轴抗压强度；m为反映岩石的软硬程度的系数；s为反映岩体的破碎程度的系数；m与s均为由岩体工程分类和三轴试验数据总结得到的与岩体破损程度有关的经验参数；mb表示与岩体特征相关的经验参数；mt表示岩石量纲为1的经验参数；a表示与岩体特征相关的经验参数；D是扰动系数；GSI为地质强度指标。

1.2 多参数协调折减方案的理论推导

假定岩土材料充满自重应力场作用下的半无限空间，各方向主应力可以由下式得出

(5)

(6)

σc=ε0E

(7)

式中：ν为泊松比；h为深度；γ为岩体重度。文献[16]提出了完整岩石抗压强度与弹性模量的关系，由式(7)表达，其中ε0为应变参数，将式(5)～式(7)代入式(1)后得

(8)

文献[17]中利用国内外关于地应力的实测资料，结合实验室三轴试验数据，通过数值方法得到弹性模量沿深度的变化规律：10 km埋深以上的土体弹性模量是沿深度线性增加的，因此假定埋深与岩土体弹性模量满足

E=E0+ah

(9)

将式(9)变形为h=(E-E0)/a代入式(8)，得

(10)

化简式(10)得

(11)

求解二元二次方程并化简方程根，舍去负根，得

(12)

化简得到Ei的表达式为

(13)

式(13)中除了νi、Ei外均是系数(参数取值参考文献[16-18])，因此Ei可以认为是νi的函数。基于以上推导，建立了弹性模量与泊松比之间的协调折减关系。

1.3 算法

1.2节中所提出的弹性模量和泊松比关系式计算安全系数的建议算法如下：

1)确定岩体的参数c、φ、m、s、E0、a。

2)根据已有参数求得β=sinφ/(1-2ν)。

3)取折减系数Ki，黏聚力、内摩擦角按ci=c/Ki、φi=arctan(tanφ/Ki)进行折减。

4)对泊松比按νi=(1-sinφ/β)/2进行折减。

5)对弹性模量按式(13)进行折减。

6)得到的ci、φi、νi、Ei代入模型计算，若边坡达到失稳破坏，输出Ki为安全系数；否则调整Ki重复第3)～6)步直至边坡发生失稳破坏。

2 非线性强度折减法算例分析

2.1 边坡算例介绍

边坡算例如图1所示，模拟过程中使用弹性和塑性两种岩土材料，采用非关联流动法则和Mohr-Coulomb屈服准则，边坡尺寸如图1所示。对该边坡算例用本文所提方法进行稳定性分析，为避免深部塑性区率先贯通使计算终止，模拟时以塑性区贯通和计算不收敛作为复合失稳判据。该边坡材料参数见表1，式(13)中各参数取值参考文献[18]。

图1 边坡模型

表1 边坡模型材料参数

按照弹性力学规定，对于边坡这样截面尺寸相比长度很小的实体，可以简化为平面应变问题。假定边坡所受外力不随纵深方向变化，位移和应变都发生在截面内。

对边坡模型沿深度分层，以深度800 m处作为分界线，模型上部岩体为理想弹塑性材料，下部为弹性材料(图1)，弹性模量随深度逐层增大，给出如下理由：

1)模型分层考虑实际边坡深部土体的压实，设置合理的弹性模量变化来模拟弹塑性变形，其塑性区的发展、应力的分布符合实际情况。

2)在式(13)的推导过程中涉及弹性模量随深度变化的关系式，因此对边坡模型沿深度分层以保证推导过程中各弹性模量保持一致。

3)考虑模型分层，塑性区下部的弹性材料单元可以产生一定的垂直变形和水平变形，减小了由于边界效应在边坡下部出现的塑性区。

2.2 折减方案对比设计

为了证明多参数协调折减的合理性，根据两种折减办法设计实验方案进行比对：①抗剪强度折减法，仅考虑折减黏聚力c和内摩擦角φ；②基于本文提出的泊松比和弹性模量折减关系，同时考虑黏聚力c、内摩擦角φ、泊松比ν和弹性模量E的多参数协调折减方法。针对以上折减方法设计如下方案：

方案1：根据抗剪强度折减法，首先取Ki=1.2，对黏聚力c、内摩擦角φ两个参数按传统抗剪强度折减法进行折减。分别取度折减系数分别为Ki=1.4、Ki=1.6、Ki=1.8、Ki=2.0、Ki=2.2、Ki=2.4、Ki=2.6、Ki=2.8、Ki=3.0的边坡材料本构模型。折减系数与对应调整的输入参数见表2。

表2 传统抗剪强度折减法折减后的参数

方案2：根据本文提出的多参数折减方案，对照流程图，首先取Ki=1.2，对黏聚力c、内摩擦角φ、泊松比ν和弹性模量E折减。同上，定义强度折减系数分别为Ki=1.4、Ki=1.6、Ki=1.8、Ki=2.0、Ki=2.2、Ki=2.4、Ki=2.6、Ki=2.8、Ki=3.0的边坡材料本构模型，调整后的参数c、φ、ν和E见表3。

表3 本文折减方案折减后的参数

2.3 计算结果与分析

2.3.1 基于折减方案1的计算结果与分析

输入表2中各折减系数对应的参数进行求解，读入不同折减参数下的结果文件，在后处理中显示边坡塑性应变云图，查看不同折减系数对应的模型塑性区分布，如图2所示。

图2 传统抗剪强度折减方案折减系数对应的累积塑性应变分布

折减系数Ki=3.0时，迭代计算不收敛，模拟终止。为找到准确的折减系数，对折减系数在2.8～3.0区间内细分，细化折减参数折减后的黏聚力和摩擦角，见表4。

表4 抗剪强度折减法细化折减后的参数

各折减系数对应边坡模型的塑性区分布如图3所示，可以看到，当折减系数取Ki=2.92时，折减的参数使模型塑性区完全贯通，此时边坡将发生失稳破坏，塑性区上部的土体将沿贯通的塑性区滑动，因此方案1安全系数取为2.92。

由图2和图3可以看出，除了滑动面的塑性区贯通，模型深部塑性应变明显，甚至先于滑动面贯通，这明显不符合实际工程情况。郑宏等在文献[11]中指出，坡底的塑性区是由于低应力水平作用下，由于边界约束效应引起的破坏，与弹塑性材料的选取有关。

2.3.2 基于折减方案2的计算结果与分析

按照方案2的折减方法，对模型输入表3参数进行稳定性分析。在迭代计算收敛的情况下，在后处理中查看不同折减系数下模型塑性区分布，如图4所示。

图4 本文方案折减系数对应的累积塑性应变分布

折减系数Ki=3.0时，模型计算不收敛，为找到使边坡发生失稳时的折减系数，在2.8～3.0区间内对折减系数细分，对应细化折减后的参数见表5。

表5 本文方案细化折减后的参数

细分后各折减系数对应模型的塑性区分布如图5所示，可以看到，当折减系数Ki=2.86时，输入的计算参数可以使边坡模型塑性区完全贯通，沿贯通的塑性区边坡发生失稳破坏，因此方案2所得折减系数取为2.86。

图5 本文方案细化折减系数对应的累积塑性应变分布

2.3.3 小结

对比两种方案的计算结果，方案2求解得到的安全系数偏小；方案1模型破坏过程中底部塑性区明显，方案2有效消除了底部不合理塑性区；对于实际边坡破坏时坡顶出现的拉应力区，两种方案都没有体现出来，对比见表6。

表6 两种方案的对比

3 结论

基于当前有限元强度折减法研究成果及其存在的不足，基于Hoek-Brown准则，考虑岩体弹性模量的沿深度递增关系，提出了一种针对弹性模量进行折减的强度折减法；结合边坡算例进行稳定性计算，通过ANSYS建立有限元模型，以塑性区贯通为边坡失稳判据，模拟岩石边坡的渐进劣化失稳破坏，最后将结果与抗剪强度折减法所得到的安全系数进行比较，得出以下结论：

1)在强度折减计算过程中，黏聚力c、内摩擦角φ、泊松比ν和弹性模量E的折减不是任意的，而是应该按照一定的协调关系进行同步折减。

2)本文提出的多参数协调折减方案得到的安全系数略小于传统抗剪强度折减法计算的安全系数，稳定性评价较为严格，消除了传统方法折减后模型底部出现的不合理塑性区，表明本文提出的多参数协调折减方案相比传统抗剪强度折减法贴近工程实际，可以有效指导边坡防治。