邹 烨，严东方

(中交第二公路勘察设计研究院有限公司， 湖北 武汉 430056)

临界滑动面的确定是边坡稳定性分析过程中的基础问题，目前大部分边坡稳定性分析方法都是以某一滑体或滑面为对象进行的。定义边坡临界滑面的方法有多种，邵国建等[1]研究了利用干扰能量等值面图确定临界滑动面。孙冠华等[2]提出临界滑面上的点为沿深度方向上等效塑性应变的极大值点。而目前在工程中应用最广泛的是基于最小安全系数的临界滑面定义方式。

自然界中的滑坡多呈现三维状态，三维极限平衡分析方法直接从二维情况下发展而来，在这一过程中为求解方程组引入的假设削弱了方法的理论与适用性[3]。陈祖煜等[4-5]建立了三维极限平衡法的上、下限分析体系。郑宏[6]提出了能适应任意形状滑动面并能满足所有平衡条件的严格三维极限平衡法。矢量和法最早由葛修润院士于1983年提出，该方法抓住力是矢量这一基本属性，在极限平衡分析与数值分析中均有应用[7-9]，且能很方便的扩展到三维情形。

目前利用矢量和法进行边坡的稳定性分析时，或结合地质调查情况与研究者经验确定潜滑面，或参照极限平衡分析法的临界滑面，而并没有通过滑面搜索确定矢量和法的临界滑面。对此，本文基于边坡应力场，研究利用矢量和法进行三维边坡临界滑面的搜索问题。

1 三维矢量和法理论

矢量和法用滑面上滑动力与抗滑力的矢量和来求解安全系数，求解图示见图1。

其安全系数定义为：

(1)

图1三维矢量和法计算图示

其中：R为滑面上的抗滑力矢量和在滑动趋势方向上的投影；T为滑面上的下滑力矢量和在滑动趋势方向上的投影。R与T的计算见式(2)、式(3)：

(2)

(3)

式中：dS为滑面上的面元。

σr与σt为滑面上任一点的抗滑力与下滑力，分别按式(4)、式(5)进行计算：

σr=σs1+σn1

(4)

σt=σ·n

(5)

σt同时可以表示为：

σt=σs2+σn2

(6)

式中：σn1、σs1为该点的抗滑正应力与极限抗滑剪应力；σn2、σs2为该点的下滑正应力与下滑切应力；σ为该点的应力张量；n为滑面上该点的外法向量(以指向滑体内部为正)。σn2可由σt表示为：

σn2=(σt·n)n

(7)

σn1与σn2为一对作用力与反作用力，两者大小相等，方向相反，即：

σn1=-σn2

(8)

极限抗滑剪应力σs1的计算公式为：

σs1=(c-σn1·tanφ)·t

(9)

式中：σn1为该点正应力的大小，并采用拉正压负的假定；t为该点极限抗滑剪应力的方向向量。

郭明伟等[8]对三维情形下的矢量和法极限抗滑剪应力方向t进行了研究。t的计算方式表述为：初始滑动趋势方向D0在该点切平面上投影方向的反方向。其中，初始滑动趋势方向D0定义为各dS面上实际剪应力的合力矢方向，即：

(10)

D为滑动趋势方向，其定义为：滑面上极限抗滑剪应力方向的反方向，即：

(11)

2 搜索算法

临界滑面的确定本质上是一个区域极值的求解问题，目前很多研究工作集中在搜索算法的选择上，常用搜索算法有遗传算法[10-11]、蚂蚁算法[12]、神经网络算法[13]、模拟退火法[14]等。本文采用的单纯型调优法是一种局部化的搜索方法，具有简单高效的特点[15]。

以安全系数Fs为目标函数，对于确定的边坡应力场，影响安全系数大小的因素就是滑面的位置。

初始滑动面一般通过一系列综合分析即可得到。事实上，对实际工程来说，利用勘察的方法往往能够确定潜在滑动面的大致位置。

对初始滑面X0：

X0={x1，x2，…xn}T

(12)

式中：x1，x2，xn可以为滑面上点的坐标；滑面形状已知的情况下，也可以是滑面方程中的参数。

利用初始单纯型X0构造初始单纯型S0：

S0={X1，X2，…X1…Xn}

(13)

其中，

(14)

(15)

(16)

式中：t为单纯型的步距，可以根据实际情况进行调整。

在初始单纯型S0的基础上，经过反射、扩张、压缩、收缩等手段，让单纯型进行翻滚、变形并逐步地向目标点靠拢，直到目标函数值(即安全系数Fs)满足停机准则方可停止。

本文的停机准则为：除去令安全系数取最大值的滑面Xh(即单纯型中令目标函数取最大值的点)，单纯型中其余各顶点的目标函数值(安全系数)的均方差小于等于1×10-4(该数可以根据需要进行选取，也可选其他合适的数)时，停止搜索，见式(17)。同时，取此时单纯型中目标函数的最小值为最小安全系数，对应的滑面即为临界滑面。

(17)

(18)

本文的基本搜索流程见图2。

图2搜索流程图

本文利用单纯型法，从二维边坡入手，利用两个算例分别进行二维与三维边坡的滑面搜索研究。

3 算 例

3.1 算例1-二维算例

许多边坡稳定性分析方法都用澳大利亚计算机应用协会(ACADS)的边坡考题做过验证。该例取自ACADS第一个考题EX1(a)。EX1(a)是一个均质边坡，数值计算所用的材料参数见表1。

表1 考题EX1(a)材料参数

首先建立模型：在有限元软件ABAQUS中建立边坡计算模型，采用ABAQUS中的平面应变单元(CPE4)，单元数目为5 459个，边坡材料服从摩尔-库仑强度准则。考题EX1(a)的有限元计算模型见图3。

图3算例1有限元模型

利用SLIDE软件，可得利用极限平衡法求解得到的临界滑动面位置，并以此滑动面为初始滑动面进行滑面的搜索。

本文使用多项式拟合二维滑面曲线，为保证精度，可采用了尽可能多的点。本文取初始滑面上的12个点来构造初始单纯型。

初始单纯型即为：

X0={x1，y1，x2，y2，…x12，y12}T

(19)

t取0.5，根据式(13)—式(16)即可求得初始单纯型S0。为节省篇幅，本文不再给出S0的具体计算过程。

滑面搜索过程见图4与图5。

图4 滑面搜索过程

图5二维滑面搜索过程收敛曲线

矢量和法与其他三种常用极限平衡条分法计算所得的临界滑面位置对比见图6。

图6临界滑面位置对比

安全系数的计算结果可与澳大利亚岩土工程协会公布的裁判答案进行对比，见表2。

表2 二维边坡搜索结果对比分析

需指出的是，表2中滑动趋势方向θ为与水平方向的夹角，而极限平衡法是边坡稳定性分析软件SLIDE进行的；笔者认为，该软件在计算时没有考虑材料的弹性，利用极限平衡法求安全系数时，可能会使结果偏小。

3.2 算例2-三维算例

在二维分析的基础上，采用相同的思路进行三维边坡滑动面的搜索。采用Zhang[16]的三维边坡算例，对自编的三维矢量和法程序进行验证。有许多研究人员曾利用该算例来验证各自的边坡计算方案，因而这里使用该算例也具有更高的可信度。

该边坡的坡比为1∶5，坡高为12.2 m，滑动面的形状为一个椭球面，其方程为：

(20)

边坡的尺寸以及滑动面在边坡中的位置见图7；边坡模型在Y轴方向的尺寸为200 m，滑面位于边坡的中间位置。

图7三维边坡外形轮廓(单位：m)

在ABAQUS中建立有限元模型，见图8。三维有限元模型采用的是八节点六面体单元(C3D8)，共划分了71 190个单元。

图8三维边坡有限元模型

该三维边坡为均质边坡，数值计算所用材料参数见表3。利用ABAQUS软件进行有限元分析得到边坡的应力场。

表3 三维边坡算例材料参数

采用单纯型优化法搜索滑面并进行安全系数的计算，滑面搜索仍然以椭球面进行。

单纯型法的优化变量为椭球面的球心坐标以及各轴的轴长。考虑到对称性，球心Y坐标不变，初始滑面的确定共使用5个变量，即：

X0={36.6，27.4，24.4，24.4，66.9}T

(21)

t取1.0，根据式(13)—式(16)即可求得初始单纯型S0。

滑面的搜索过程见图9与图10。

图9 三维滑面的搜索过程

图10三维滑面搜索过程收敛曲线

将矢量和法的临界滑面与初始滑面以及文献[17]中搜索得到的滑面进行对比，见图11。

图11三种方法临界滑面的对比

其中，林永生等[17]利用遗传算法进行滑面搜索，在搜索过程中，只改变椭球体球心的X坐标与Z坐标，而轴长均不发生改变，也就是说，椭球体的大小及形状均不发生改变，而只改变椭球体的位置，最终所得到的滑面方程为：

(22)

本文的搜索只对滑面的形状做了限定，即要求滑面为一个椭球体，对轴长与椭球的球心位置没有限制。通常限定条件越少，其结果的与实际也会越接近。通过搜索最终得到的滑面方程为：

(23)

比较式(22)、式(23)，可以看出，在限定滑面为椭球面的条件下，三维矢量和法的临界滑动面与Zhang[16]的计算结果有一定的区别，其中矢量和法与Zhang的滑面在坡顶部分比较接近，而坡底区域，矢量和法的滑面更加的靠近坡脚。

将矢量和法计算得到的安全系数与各文献所得结果进行对比，所得结果见表4。从该表可以看出：当矢量和法采用与式(20)的滑面进行计算时，其安全系数偏差为3.86%。而经过滑面搜索，其安全系数略有降低。值得指出的是，表4中关于误差的计算是以Zhang的计算结果为基准进行的，但并不代表其为标准答案。

表4 三维边坡搜索结果对比分析

事实上，郑宏[6]指出，Zhang的方法忽略了三维条块上四棱柱沿滑动方向上两个侧面上的剪应力可能导致其计算结果偏小。可以认为，利用矢量和法通过滑面搜索得到的计算结果与实际情况更为接近。

4 结 论

本文研究了三维矢量和法临界滑动面的搜索问题，给出了滑动面搜索的全过程，得到以下几点结论：

(1) 矢量和法能很方便的进行二维与三维边坡的稳定性分析，在边坡应力场求解准确的情况下，该方法所得结果具有较好的可靠性。

(2) 二维情形下，利用矢量和法求解极限平衡分析临界滑面的结果与矢量和法临界滑面的结果相近；三维情形下，进行滑面搜索后能更加准确的进行边坡的稳定性分析。