基于局部安全系数的边坡滑面搜索

王东英，杨光华，姜 燕，靳晓兵

(1. 广东省水利水电科学研究院，广东 广州 510635；2. 广东省岩土工程技术研究中心，广东 广州 510635；3. 广东省乐昌峡水利枢纽管理处，广东 韶关 512200)

边坡稳定性分析主要包括两方面工作：一是确定潜在的最危险滑面，二是通过边坡局部或整体安全系数，分析指定位置处的稳定状态。极限平衡法和强度折减法是当前边坡稳定性分析中最常用的两种方法。极限平衡法基于力或力矩平衡求解安全系数，一般需事先假定滑面形状，并通过遍历计算各滑面安全系数，遴选出安全系数最小的滑面。该方法原理简单易懂，但计算耗时大，且在假定滑面形状时需要丰富的工程经验[1-4]。强度折减法是在有限元分析的基础上，按假定的衰减规律，折减边坡体的强度参数直至计算不收敛的方法。当前基于强度折减法确定临界滑面的方法繁多且互不关联：①基于塑性区或剪应变增量云图[5-7]；②基于应力或位移等值线[8-9]；③基于位移变化率或等效应变最大值[10-11]。由于各个方法立论依据不一，很难区分哪种方式最适合强度折减法。不仅如此，强度折减法相应的参数折减规律一直是争议较大的问题，其对应的虚拟应力状态是否合理也被质疑。杨光华等认为应更全面的从应力位移场角度来认识边坡稳定问题，可以获得更全面的认识[12-13]。

事实上，岩土工程的破坏多是由局部破坏逐渐扩展为整体破坏。边坡的失稳也是因个别点超出容许应力首先破坏，而后由点及面逐步演化、扩展，最终导致边坡整体滑动的渐进破坏过程。边坡的整体稳定性受各局部位置稳定性的影响，且当各处安全系数均达到最小值时，边坡的整体安全系数也达到最小值。因此，通过局部点的安全性评估边坡的临界滑面及其安全性理论依据充分，合理可行。

基于上述认识，探究了一种基于局部点的安全性由点及面确定边坡滑面的方法，并基于各局部安全系数推导了整体安全系数表达式。其中，由点及面的搜索过程应与边坡的破坏类型相符合。牵引式破坏应自坡脚至坡顶逐步搜索；推动式破坏则相反。以牵引式破坏为例，通过算例分析验证了滑面搜索和安全系数求解的可靠性，并给定布点间距的合理范围。

1 基础理论

1.1 边坡中的Cauchy问题

以二维边坡为例，若定义曲线为滑面方程，并给定边坡滑面的起始点(x0，y0)以及边坡滑面扩展的控制方程，则可唯一的获取滑面曲线方程[14]。

边坡滑面的起始点应与边坡坡面局部安全系数最小的位置相对应，滑面扩展控制方程的选取应以最优性原理为原则。控制论中的最优性原理表明，最危险滑面应为滑面上各处安全系数均达到最小时对应的滑面。因此，控制方程应以局部安全系数取最小值的方向为准则。

假定边坡滑面起点坐标为(x0，y0)，以局部安全系数取最小值的方向作为控制准则，则边坡滑面扩展的控制方程应为：

(1)

式中ψ为局部安全系数最小值方向与x方向的夹角。

由于ψ的数值与该点的应力状态(σ1，σ3)有关，因此式(1)实质为一阶拟线性偏微分方程[15]。只要准确给出滑面起点坐标及合适的布点间距，在应力场精度满足要求的前提下，可以唯一的确定临界滑面位置。

1.2 局部安全系数及其方向

假设二维边坡中某点的应力状态为(σ1，σ3)(如图1所示)，则该点在任意方向的应力状态可表示为：

(2)

假设岩土体服从摩尔库伦屈服准则，并假定以压为正，边坡局部安全系数应为：

(3)

将式(2)代入式(3)有：

(4)

对(4)式中θ求导，即：

(5)

所求θ值即为使式(4)中安全系数K取最小的方向角。分析图1，自摩尔圆A点到B点，τf与τn均在增大，自B点到C点，τf仍在增大而τn则在减小，因此，使安全系数取最小时对应的方向点应介于A、B之间，2θ的取值为90 °～180°之间，也即安全系数取最小的方向与最大主应力σ1的方向夹角在0～90°间。由于上述表达式求导复杂，采用MATHEMATICA推导其计算结果为：

(6)

进一步化简并计算对应安全系数，结果为：

(7)

上式中为最小安全系数方向与最大主应力的夹角，顺时针转至σ1时取正值，逆时针转至σ1时取负值。

边坡滑面搜索时需要确定最小安全系数方向与x轴正向的夹角，也即式(1)中的ψ。因而需要探究与式(7)中的关系。

根据材料力学应力状态相关理论，最大主应力与x轴或y轴的夹角公式为：

(8)

假定顺时针旋转为正，则根据不同点的应力状态，φ的大小可根据以下几种情况确定。

①当σx>σy且τxy>0，其对应的摩尔应力圆如图2(a)所示，x轴正向与σ1的夹角为α0(顺时针旋转，取正)，则对应极值角表达式为：

(9)

②当σx<σy且τxy>0，其对应的摩尔应力圆如图2(b)所示，y轴正向与σ1的夹角为α0(逆时针旋转，取负)，x轴正向与σ1的夹角为π/2-α0，则对应极值角表达式为：

(10)

③当σx>σy且τxy<0，其对应的摩尔应力圆如图2(c)所示，x轴正向与σ1的夹角为α0(逆时针旋转，取负)则对应极值角表达式为：

(11)

④当σx<σy且τxy<0，其对应的摩尔应力圆如图2(d)所示，y轴正向与σ1的夹角为α0(顺时针旋转，取正)，则x轴正向与σ1的夹角为α0-π/2，其对应极值角表达式为：

(12)

其中β1、β2分别为x轴正向与上、下半应力圆安全系数取最小的方向的夹角。

1.3 基于局部安全系数的边坡滑面搜索方法

众多数值分析算例表明[15-16]，对于均质边坡，多数情况下是坡脚附近区域最先破坏，而后引起邻近区域的相继失稳，也即前述的牵引式破坏。以牵引式破坏为例，基于局部安全系数搜索边坡滑面的步骤阐述如下：

1) 建立数值分析模型，通过有限元或有限差分方法求解自然状态下边坡的应力场；

2) 设置搜索间距：在边坡模型内布置n条等间距的垂直线段，垂直线段要求覆盖整个边坡(如图3所示)。各垂直线段间距为ΔL=L/n，其中L为边坡总长度；

3) 确定初始离散点位置：通过式(7)计算坡脚附近边界各局部安全系数值及其极值角ψ，并将安全系数最小者对应的位置坐标设为起点A(x0，y0)，其对应的极值角正切值记为tanψA；

4) 确定滑面其他离散点位置：根据A点坐标及其最小安全系数方向，通过式(13)求解B点坐标(xB，yB)。并应用插值手段获取B点应力状态，代入(7)求解B点的安全系数及极值角正切值tanψB，继而求得下一搜索点C的坐标，重复应力插值及安全系数求解，获取C点的安全系数及极值角大小。依次类推，直至求解的第m个点不在边坡内为止；

(13)

5) 形成滑面：连接各个搜索点得到边坡的临界滑面。搜索流程如图4所示。

对于推动式破坏，初始点应为坡顶边界局部安全系数取最小值处。相邻点坐标求解应相应调整。

2 边坡整体安全系数求解

如前所述，边坡滑面为由点及面逐渐演化、扩展而成。局部区域的安全性势必影响边坡整体的稳定性。A Kourdey等[15-16]均指出，推导基于局部安全系数的边坡整体安全系数意义重大。

若将整个滑面按布点间距划分为若干微段，根据极限平衡条件，对于每一微段应有下式成立：

(14)

式中ΔLi为第i微段的长度，τfi为微段抗剪强度，τi为微段剪切力，Fi为微段局部安全系数。

则整个滑面上的下滑力应为各微段下滑力之和。若考虑力的方向，该和为各下滑力在投影方向的代数和；若不考虑力的方向，则该和为各微段下滑力的简单代数叠加，此时整个滑面上的下滑力可表示为：

(15)

根据安全系数的定义，边坡整体安全系数应为：

(16)

若考虑力的方向性，按矢量和法求得的整体安全系数公式为：

(17)

相较于文献[17-18]有关局部安全系数与整体安全系数的假定，本文公式推导从安全系数定义出发，物理意义相对明确。

3 算例分析

3.1 方法可靠性验证

为验证所提方法的合理性，分别选用文献[10]和文献[11]提及的均质边坡算例分析边坡滑面，并将结果与经典的的Bishop法及参考文献的结果作对比。其中文献[11]中的滑面位置通过最大位移变化率的方法确定；文献[10]则是基于等效塑性应变方法获得的滑面位置。两种方法均基于强度参数同比例折减后的虚拟应力场，且对滑面点进行拟合时对拟合基函数的要求较高。

算例1为一均质土坡，取自澳大利亚(ACADS)边坡考题，其材料特性为：容重20 kN/m3，粘聚力3.0 kPa，内摩擦角19.6°边坡模型尺寸如图5所示，其中H=5 m，a:b=1:2。

采用有限元软件ABAQUS求解应力场时，边界条件为底边固定约束，两侧边法向约束；而本构条件及屈服准则分别为弹塑性本构、Mohr-Coulomb屈服准则。

求得自然状态下的应力场后，按照1.3节所述步骤开展滑面搜索工作，滑面搜索结果如图6所示。

图6展示了推荐方法、基于最大位移变化率方法及Bishop法得到的滑面结果。对比发现：3种方法所得滑面位置相近，但形状有一定差别。本文方法在剪出口较为平缓，且剪入口处相对靠后，与Bishop法的圆弧滑面有交叉。

在边坡的安全性方面，文献[11]的计算结果基于强度折减法，对应的安全系数为0.989；极限平衡的Bishop法相应的结果为0.986；本文推荐方法所得边坡整体安全系数为0.963(基于式(16))，小于上述两方法的结果，结果相对合理。

算例2为一均质黏土边坡。黏土容重为20 kN/m3，弹性模量为200 MPa，泊松比为0.4，粘聚力为10 kPa，内摩擦角为20°。边坡模型尺寸如图7(a)所示，推荐方法、基于等效塑性应变最大值法及Bishop法所得滑面位置如图7(b)所示。

图7(b)结果表明，推荐方法、极限平衡法及基于等效塑性应变确定的滑面位置基本一致，且推荐方法相应滑面介于其他两方法之间，在坡脚位置推荐方法滑面相对较深，符合黏土滑面的基本规律。各方法对应的边坡稳定性系数分别为：强度折减法1.380，极限平衡法(Bishop)1.361，推荐方法1.348，其数值也小于前两方法，结果相对合理。

推荐方法所得滑面及边坡安全性与当前已有方法结果的一致性表明，应用推荐方法分析边坡稳定性是合理的。因其原理的合理性和操作的简便性，有推广应用的价值。

3.2 方法应用探讨

本文推荐的滑面搜索方法基于边坡真实应力场，只要能够准确计算出边坡的应力场，就可获取较为理想的滑动面位置。但潜在滑动面的精确定位显然依赖于搜索点间距的取值，因此，有必要对搜索点间距的合理取值进行探究。现以算例1为例，分别取用0.2 m、0.8 m、2 m 3个间距对搜索点间距的取值范围开展研究。

不同搜索点间距对应滑面位置如图8所示。由图8可知，间距为0.2 m、0.8 m时搜索得到的滑面位置较为接近，间距为2 m时得到的滑面位置稍深，但相应滑动面位置仍很靠近小间距滑面。笔者也尝试3 m、4 m两个间距，但相应的滑面位置偏差较大，在此未列出。

由于间距2 m约等于边坡长度(即坡面长度的水平投影)的0.1倍，且所得滑面位置与小间距相比偏差不大，推荐搜索点间距的合理取值不大于0.1倍坡长。

4 结论及展望

本文依据控制论中的最优性原理和边坡中的Cauchy问题，探究了一种基于局部安全系数由点及面确定边坡滑面的方法。并从安全系数的定义出发推导了基于局部安全系数的边坡整体安全系数表达式。通过对比分析推荐方法、极限平衡法、基于最大位移变化率或基于等效塑性应变的强度折减法获得的滑面位置及安全系数结果，验证了推荐方法的合理性，所得结论主要有：

1) 对于规则的均质边坡，推荐方法能够高效获取临界滑面位置，并求得可靠的整体稳定性系数；

2) 推荐方法的布点间距控制在0.04倍边坡总长度内所得滑面结果较为可靠；

3) 与极限平衡法相比，滑面的确定无需遍历搜索；与强度折减法相比，滑面搜索及安全系数求解建立在真实应力场基础上，且滑面搜索方法依据充分，规则统一；

4) 本文所提方法是基于边坡应力场进行计算的，适用于均质岩土边坡。对于多层介质边坡，尤其是含软弱夹层边坡，由于应力场在介质交界处存在不连续现象，安全系数及极值角在交界处也存在不连续现象，需要通过加入连续性条件控制滑面搜索结果。在这方面，尚需进一步研究。