基于FLAC的极限状态黄土高边坡强度参数反演

万黎明，蔺晓燕，余宏明，曾红彪，赵 超

(1.中国地质大学(武汉) 工程学院，武汉 430074； 2.西安石油大学 地球科学与工程学院，西安 710065)

基于FLAC的极限状态黄土高边坡强度参数反演

万黎明1，蔺晓燕2，余宏明1，曾红彪1，赵 超1

(1.中国地质大学(武汉) 工程学院，武汉 430074； 2.西安石油大学 地球科学与工程学院，西安 710065)

评价边坡稳定性关键问题是确定滑动面和可靠的抗剪强度参数。基于数值模拟技术，提出一种新的方法，可在无需事先假定滑动面的情况下得到可靠的强度参数。通过对兰州地区30个极限状态黄土边坡的坡高和坡宽参数进行统计，发现二者呈线性关系，且拟合度较高。根据所得到的回归方程，采用FLAC建立边坡模型，利用强度折减法对兰州地区边坡进行强度参数反演，将4个邻近高度边坡得到的参数反演曲线绘于同一坐标系中，得到该坡段综合强度参数值。反演结果表明：随着坡高的增加，黄土内摩擦角逐渐减小，黏聚力增大。根据计算结果，绘制60～120 m边坡高度下强度参数综合反演曲线，可直接进行黄土高边坡稳定性评价，提高了工作效率，也为高边坡几何形态的设计提供可靠的数据支持。

黄土高边坡；FLAC；强度折减法；参数反演；安全系数；兰州地区

1 研究背景

黄土边坡稳定性一直是工程地质领域研究的重要内容之一，其关键问题是确定边坡的危险滑动面和安全系数。而边坡安全系数是通过抗剪强度参数计算所得，因此合理地确定强度参数值决定着边坡稳定性的正确评价。工程中确定非饱和土体强度指标的方法主要有经验法、土工试验法和反分析法。工程类比法是一种经验估算方法，由于滑坡成因、几何边界、岩土体性质和研究者的工程经验存在着差异性，工程类比法难以得出准确的抗剪强度参数值。采用试验法确定边坡抗剪强度参数是基于Bishop(1959年)提出的非饱和土有效应力原理来进行抗剪强度计算的[1]，但由于在取样、试验中的局限和干扰，土-水特征曲线及基质吸力参数确定的困难和技术的不完善性，导致试验值与实际的抗剪强度参数不相符合，影响边坡稳定性的评定。因此，本文利用反演分析法来计算兰州地区黄土高边坡的抗剪强度参数值。

进行强度参数反演的前提是边坡安全系数值为1.0，据此李萍等[2]提出了极限状态边坡的概念及其4条判定标准。赵尚毅等[3]将有限元静力平衡方程组是否有解、有限元计算是否收敛作为边坡破坏的依据，定量判定极限状态边坡的存在。但是陈力华等[4]发现这种判据与其他方法对陡边坡的判断结果存在差异，并对判定结果产生分歧的原因进行了研究。房智恒等[5]利用FLAC3D建立边坡模型，根据时步最大节点位移曲线计算单元安全度相结合的边坡稳定性分析方法来判断弧状边坡极限状态。由于现场试验强度参数的不可靠性，基于数值模拟定量判定极限状态边坡安全系数是否为1.0也会有较大的误差，得出的结论与现实情况不太符合，因此本文利用自然历史分析法对研究区边坡极限状态进行评定。

刘超[6]对甘肃地区黄土高边坡采用Morgenstern-Price法进行强度参数反演研究，并计算了该区域边坡的失效概率，但是Morgenstern-Price法计算安全系数需要事先假定边坡滑移面的形态和位置，且无法提供边坡节点的位移信息，具有一定局限性。刘文红等[7]利用Bishop法对边坡进行参数反演，但Bishop法是基于圆弧状滑动面进行计算的，对滑动面的形状有较高的要求。因此，本文提出了一种新的模型建立方法和稳定性计算方法对黄土高边坡的强度参数进行反演分析，无需假定潜在滑移面，并且能较好地反映土体受力情况和运动过程，获取节点位移变化和变形信息，更为直观地判断边坡的贯通破坏区。最终绘制边坡强度参数随坡高变化规律图，建立不同边坡高度和宽度下的强度参数综合反演曲线，为边坡设计和加固提供可靠的数据支持。

2 参数反演原理

2.1 边坡变形状态判定

进行强度参数反演的第1步即确定边坡所处的变形状态，根据安全系数可以将滑坡分为4个不同的变形阶段，如表1[8]。其中，边坡处于极限平衡状态是进行强度参数反演的前提，即安全系数Fs的取值为1.00。

表1 滑坡不同阶段的安全系数

对于安全系数为1.00的边坡，李萍等[2]提出了极限状态边坡的概念及其4条判定标准，即：①坡顶有拉裂隙；②坡面局部滑塌，地形破碎；③现有新滑坡恢复的原始坡形；④新滑坡两侧与其坡高、坡度相当的边坡。根据以上4条判定标准，李萍等[9]对甘肃地区266个极限状态边坡高度和边坡宽度等数据进行了实测。同时基于FLAC数值模拟技术也可定量判断边坡所处状态。利用FLAC强度折减法对边坡达到极限平衡状态的判据主要有3种：①限定求解迭代次数，当超过某限值仍未收敛则认为破坏发生[10]；②限定节点不平衡力与外荷载的比值大小[11]；③利用可视化技术，当剪切塑性区自坡角到坡顶贯通则定义为坡体破坏。

由于黄土的强度性质最主要受其内部结构影响，根据试验所得的强度参数具有不确定性，因此采用FLAC建立模型来定量计算极限状态边坡的安全系数也会有较大的误差。为了确定研究区边坡是否处于极限状态，本文基于李萍等提出的4条判定标准，结合自然历史分析法，对研究区区域地质背景进行研究，对高边坡变形迹象、规模以及形态进行勘察，对部分边坡进行采样测定其物理参数，并分析了控制其变形破坏的主导因素，确定了兰州地区30个处于极限状态的边坡(如图1所示)。表2为兰州会宁地区坡高为69.8 m、坡度35.4°的极限状态边坡的物理力学参数值。

2.2 FLAC参数反演思路

图1 极限状态边坡测点位置分布

表2 研究边坡的物理力学参数

注：兰州地区马兰黄土重度差异性不大，根据统计分析取平均值为γ=1 800 kN/m3

边坡失稳使得滑动体由静止状态变为运动状态，同时沿着滑移面产生较大的位移。FLAC软件通过强度折减法使得边坡达到极限破坏状态，此时的位移和塑性应变都将会产生突变，不再是一个定值。程序无法从数值方程组中找到一个既能满足静力平衡，又能满足应力-应变和强度准则的解，不管是从位移的收敛条件和力的收敛准则来判断数值计算都不收敛[12]。根据边坡破坏之前计算收敛，破坏之后计算不收敛来表征土体的移动，因此可以将数值计算是否收敛作为边坡破坏的依据。本文按照边坡收敛计算其安全系数，设置力不平衡比率R为1×10-5作为收敛的界限值，以R<1×10-5表征收敛，设定安全系数上限值与下限值之间(即k12-k11)差值为0.001作为退出循环的限定迭代次数条件，以保证安全系数的精度要求。最后计算的安全系数上限值与下限值的均值为最终的安全系数值，即

(1)

FLAC参数反演的基础是边坡计算收敛，极限状态边坡安全系数值为1.0。其基本思路为基于试验值和经验值设定一个强度参数，本次设定内摩擦角φ为0～30°范围内的定值，代入所建边坡模型当中；反复调整另一个强度参数值c，利用强度折减法进行计算，直至计算所得的安全系数值为1.0，此时的强度参数即反演参数值。改变内摩擦角φ值，再进行参数反演，得到另一组黏聚力值，这样每一个边坡高度都可以得到多组强度参数反演值，将所得到的数据进行拟合，发现φ与c符合一定的对数关系。对于邻近边坡高度参数反演曲线绘于同一坐标系中，相邻曲线交于一点，为边坡段的综合强度参数反演值。

3 FLAC强度参数反演

3.1 FLAC强度折减法

在判定边坡状态后，利用极限状态边坡安全系数Fs=1.0进行边坡强度参数反演，利用FLAC数值模拟软件并采用强度折减法进行稳定性计算。强度折减法即通过式(2)、式(3)改变折减系数，从而调整边坡强度参数值反复折减计算，直到边坡达到极限状态，此时的折减系数即是边坡的安全系数Fs[13]。

(2)

(3)

式中：ct，φt分别为折减后的黏聚力和内摩擦角；Ft为折减系数。

通过折减系数法可以直接计算边坡的安全系数，不用指定滑移面的形状和位置，并可以监测边坡节点的位移情况。该方法弥补了Morgenstern-Price法人为假定滑移面，无法真实反映土体的受力情况的不足[14]，同时还能显示边坡的位移和变形信息，能直观地观察到塑性贯通区和破坏情况，为边坡极限状态判断提供依据。

利用FLAC强度折减正演计算顺序为：根据已知边坡形态建立边坡模型，进行网格剖分，设定计算边界条件和本构模型，利用Fish语言对参数进行赋值，最后计算边坡安全系数。因此参数的取值正确直接关乎到最终安全系数值的计算正确，影响整个边坡稳定性评判。但是由于非饱和土的特殊性(具有负孔隙水压力的基质吸力)，在试验和取样过程中的干扰，参数确定的困难，导致经过试验所得的边坡强度参数计算出来的安全系数往往与实际不符合。为了得到可靠的强度参数值，利用极限状态边坡的安全系数值为1.0，加上黄土边坡各段c值变异性较大，φ值变异性很小，因此通过已知安全系数值，设定φ值来反算c值，是解决黄土强度参数选取的有效方法。

3.2 极限边坡高度和宽度拟合

在进行参数反演时，为便于极限状态边坡模型的建立，本文对现场所得30个极限状态边坡高度和边坡宽度数据进行统计分析，发现其符合直线分布，即y=1.3x-7.3，相关性系数R2=0.879，这样即可绘制极限状态边坡的坡高、坡宽拟合曲线(图2)。

图2 极限状态边坡高度和宽度拟合曲线Fig.2 Fitted curve of slopeheight and slope width at limit state

从图2中可以看出，极限状态边坡高度和宽度拟合度较高，相关性系数接近于1。通过拟合曲线，对任意极限状态边坡高度或者宽度都可以计算其边坡宽度或者边坡高度，从而有利于FLAC边坡模型的参数选取与建立。为便于计算，后文中边坡高度取整数，每间隔一定高度利用回归方程计算所对应的边坡宽度，从而建立不同断面的边坡模型，进行参数反演。

3.3 参数反演过程

均质邻近高度边坡有相似的断面，在相似条件下2个或者多个断面的方程可以联立求解[15]。当边坡高度相差较大时，物理力学参数值相差也较大，含水量、滑移状态和滑移过程等也有较大区别，滑移面条件不再相似，进行断面联立求解会产生较大的误差，反演所得参数曲线可能有多个或者没有交点。

图3 坡高60 m的边坡模型Fig.3 Model of slope of 60 m height

本文对60～120 m极限状态边坡利用FLAC建立边坡模型，为了便于计算，边坡高度统一采用整数值，边坡宽度通过拟合方程计算可得。图3为建立的高度为60 m的黄土高边坡模型，并对所建单元进行网格划分。

研究区黄土高边坡邻近断面具有相似的物理力学性质和滑移面，文中将边坡高度按照12 m进行分段，每间隔4 m进行一次参数反演，对于40 m以上的高边坡，涉及多层黄土，各层土重度取不同的值，反演时的c，φ值按均质土层考虑[16]。每一段由4个极限状态边坡得到一组反演参数值，为该区域段高度边坡的综合强度参数。因滑坡实际抗剪强度实际值应介于天然强度和饱和强度之间，而黄土边坡c值变异性很高，φ值变异性很小，不同试验方法得到的c值变化大，φ值变化较小[17-19]，所以本文进行参数反演时首先设定φ值为定值，根据试验测得力学参数作为参考，设定φ=8°，12°，16°，20°，24°。在所建模型中反复调整c值直至计算结果满足边坡处于极限状态，将所得到的c，φ值统计并绘制于坐标系中，即得到c-φ拟合曲线。这样每4个极限状态边坡的所得到的反演c-φ曲线绘制于统一坐标系中即为该区域段的综合反演参数[16]，如图4所示。

图4 不同高度边坡参数反演结果Fig.4 Parametric inversion results of slope of different heights

由图4综合反演曲线的斜率变化可知：对于同一断面，随着内摩擦角φ值的增大，对应的黏聚力c值减小，但其安全系数保持为1.0，表明同一断面不同的强度参数组合可得到相同的安全系数；对于不同断面，c值和φ值对边坡高度变化的敏感性不同。

4 反演结果分析

根据图4中16个边坡模型反演计算，将所得60～120 m边坡强度参数综合反演值统计于表3中，将所得数据拟合成曲线即得到强度参数随边坡高度变化的规律，如图5所示。

表3 60～120 m边坡强度参数综合反演值

图5 边坡强度参数综合反演曲线Fig.5 Comprehensive inversion curves of strength parameters of slope

从图5中可以得出，对于60～80 m边坡，随着边坡高度的增加，内摩擦角不断减小，黏聚力逐渐增大；当边坡高度>80 m时，内摩擦角减小缓慢。随着边坡高度的增加，内摩擦角φ不断减小，黏聚力c不断增大,这是由黄土的内部组成和外部结构所决定的，与黄土地区边坡强度参数对含水率敏感性变化规律一致，也反映了研究区的地层时代和物理力学性质的变化特点。

图6为60 m边坡内摩擦角为20°和24°下的速度矢量图和剪切位移增量等值线图，其对应安全系数值均为1.0，对比发现同一断面强度参数不同，塑性贯通区也不相同，滑移面的位置也发生了变化。

图6 60 m边坡φ=20°，24°时的速度矢量和剪切位移增量等值线Fig.6 Contours of velocity vector and shear displacement increment of slope with height of 60m and internal friction angles of 20° and 24°

图7 最大节点位移-时步曲线Fig.7 Curve of maximum node displacement vs. time step

根据反演参数计算，利用FLAC可以输出其边坡最大节点位移-时步曲线(如图7所示为69.8 m极限状态边坡的最大节点位移-时步曲线)，根据图7的收敛状况可以对边坡的极限状态进行核定。从图6中可明显观察到塑性贯通区验证了边坡的极限状态，这与利用自然历史分析法判定的结果是一致的。因此在物理力学参数准确的情况下，利用FLAC定量进行边坡参数反演是可行的。

综合反演曲线可较好反映边坡临界状态的特点，反演结果考虑了边坡滑体的边界约束条件和侧向约束作用，而不是将滑体作为一个空间无限体进行计算所得到的，因此更接近真实值，更可靠。根据60～120 m强度参数随边坡高度变化的综合反演曲线，可以直接得到任意形态极限状态黄土边坡可靠的强度参数，利用FLAC可在不用假定滑移面的情况下获得边坡的塑性贯通区和应力分布。所得参数可直接用于黄土高边坡稳定性的评判和边坡几何形态的设计，大大节约了研究者的时间，为黄土边坡防治提供了更加可靠的数据支持。

5 结 论

通过FLAC建立边坡模型，进行强度参数反演，输出速度矢量图和剪切位移增量等值线图，可以明显观察到边坡的塑性贯通区和速度矢量，避免了事先假定滑移面所造成的误差。通过本文研究,得到了可靠的边坡强度参数值，分析其变化规律及稳定性，可大大节约研究者的时间，为边坡的加固和防治提供了依据。经过算例计算和研究，可得到如下结论：

(1) 兰州地区极限状态边坡形态拟合度较高，坡高和坡宽拟合呈直线分布，且相关系数达到0.8以上。因此得到的坡高和坡宽拟合方程对大部分边坡形态均适用，为边坡模型建立创造了条件。

(2) 对兰州地区黄土边坡进行参数反演得出，随着坡高的增大，黄土边坡的内摩擦角不断减小，黏聚力逐渐增大；黏聚力和内摩擦角对不同边坡形态的敏感性不同，同一断面不同强度参数可以得到相同的安全系数，但是滑移面的位置会发生变化。

(3) 利用FLAC建立边坡模型，将边坡进行网格剖分，模拟边坡天然应力状态，具有分析渐进的失稳和破坏，尤其适合处理非线性的大变形破坏问题的优势，且不用事先设定可能滑移面，能有效分析土体的受力情况。

(4) 极限状态边坡在兰州地区是存在的并且广泛发育。在可靠参数下，利用FLAC建立边坡模型，并计算随时步变化的最大节点位移，利用所得到的速度矢量图和剪切位移增量云图，可进行边坡失稳破坏判定。根据计算收敛判据定量确定极限状态边坡的存在是可行的。

[1] FREDLUND D G， RAHARDJO H. Soil Mechanics for Unsaturated Soils[M].New York：John Wiley & Sons Inc，1993.

[2] 李 萍，王秉纲，李同录. 自然类比法在黄土路堑边坡设计中的应用研究[J].公路交通科技，2009，26(2)： 1-5.

[3] 赵尚毅，郑颖人，张玉芳.极限分析有限元法讲座——Ⅱ：有限元强度折减法中边坡失稳的判据探讨[J]. 岩土力学，2005，26(2)：332-336.

[4] 陈力华，靳晓光. 有限元强度折减法中边坡三种失效判据的适用性研究[J]. 土木工程学报，2012，45(9)：136-146.

[5] 房智恒，王李管，彭南良，等. 高陡边坡整体与局部失稳的强度折减及安全度判别分析[J]. 重庆大学学报(自然科学版)，2015，38(6)：8-14.

[6] 刘 超. 黄土高边坡抗剪强度参数反演研究[D]. 西安：长安大学，2014.

[7] 刘文红，张亚卿，李同录，等. 根据参数反演结果分析非饱和黄土边坡破坏机理[J]. 工程地质学报，2014，22(5)：915-920.

[8] 徐汉斌，王 军.反算法中滑坡稳定系数的取值问题[J]. 四川地质学报，1999，19(1)：86-89.

[9] 李 萍，黄丽娟，李振江，等.甘肃黄土高边坡可靠度研究[J]. 岩土力学，2013，34(3)：811-817.

[10]郑颖人，赵尚毅. 有限元强度折减法在土坡与岩坡中的应用[J]. 岩石力学与工程学报，2004，23(19)：3381-3388.

[11]GRIFFTHS D V，LANE P A. Slope Stability Analysis by Finite Elements[J]. Geotechnique，1999，49(3)：387-403.

[12]彭文斌. FLAC3D实用教程[M].北京：机械工业出版社，2007.

[13]郑志勇，余海兵，徐海清. 软硬岩互层边坡的破坏模式及稳定性研究[J]. 长江科学院院报，2016，33(9)：102-106.

[14]王 科，王常明，王 彬，等. 基于Morgenstern-Price法和强度折减法的边坡稳定性对比分析[J]. 吉林大学学报(地球科学版)，2013，43(3)：902-907.

[15]郑明新. 论滑带土强度特征及强度参数的反算法[J].岩土力学，2003，24(4)：528-532.

[16]李同录，刘 超，李 萍. 甘肃天水地区黄土极限状态边坡的统计分析[J]. 地球科学与环境学报，2013，35(2)：107-112.

[17]雷晓峰.黄土边坡强度参数的选取及应用[D]. 西安：长安大学，2005.

[18]张少宏.黄土边坡稳定计算中参数的敏感性分析[J]. 水利与建筑工程学报，2003，1(3)：40-42.

[19]赵永虎，刘 高，毛 举，等.基于灰色关联度的黄土边坡稳定性因素敏感性分析[J].长江科学院院报，2015，32(7)：94-98.

(编辑：黄 玲)

Strength Parameters Inversion of High Loess Slopeat Limit State Using FLAC

WAN Li-ming1，LIN Xiao-yan2，YU Hong-ming1，ZENG Hong-biao1，ZHAO Chao1

(1.Faculty of Engineering，China University of Geosciences，Wuhan 430074，China; 2.School of Earth Sciences and Engineering，Xi’an Shiyou University，Xi’an 710065，China)

Determining sliding surface and reliable shear strength parameters are crucial issues in evaluating slope stability. A new method based on numerical simulation technology is proposed to obtain reliable strength parameters without presuming the sliding surface. The relationship between heights and widths of thirty loess slopes at limit state in Lanzhou region were found to be linear，and the fitting degree was high. According to the obtained regression equation，FLAC was applied to establish the slope model, and strength reduction method was used to inverse the strength parameters of slope at limit state. Furthermore, the parameter inversion curves of four adjacent high slopes were drawn in the unified coordinate system, and the comprehensive strength parameters of the segment were acquired. The inversion results show that with the increase of slope height，the cohesive force of loess increases gradually but internal friction angle decreases. According to the calculation results，the comprehensive inversion curves of strength parameters of slope with heights from 60 m to 120 m are drawn. The method could directly evaluate the stability of high loess slope，improve the working efficiency，and also provides reliable data to support the geometry design of high loess slope.

high loess slope; FLAC; strength reduction method; parameters inversion；safety factor；Lanzhou region

2016-04-20；

2016-05-17

国家自然科学基金项目(41402274)

万黎明(1992-)，男，湖北天门人，硕士研究生，主要研究方向为岩土体工程性质及其稳定性，(电话)13026162550(电子信箱)wan_liming@126.com。

10.11988/ckyyb.20160380

2017,34(7):65-69,76

P642

A

1001-5485(2017)07-0065-05