基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法

技术领域

本发明涉及斜坡稳定性分析及预测预警技术领域，特别涉及一种基于牵引式斜坡变形破坏机理的稳定性分析及预测预警方法。 

背景技术

牵引式斜坡的稳定性计算及预测预警一直采用的是推移式斜坡的经验稳定性计算及预测预警方法，迄今为止，还没有将两者加以区分开来。牵引式斜坡的稳定性计算及预测预警方法还不完善；另外，牵引式斜坡的破坏模型随着坡体的不同组成，其模式并不相同。且破坏时，各处的位移值均不相同，与推移式斜坡相比，破坏时推移式斜坡变形值很小，并不能完全适用于牵引式斜坡，导致计算及预测结果不精确。 

发明内容

本发明的目的就是提出一种基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，在牵引式斜坡变形机理、坡体与滑面的力学特性以及滑面不同点的演化特征分析的基础上，提出了牵引式斜坡的几种破坏模式、基于变形的稳定系数计算方法及斜坡坡面位移的决定方法，从而实施预测预报。 

本发明基于牵引式斜坡变形破坏机理的稳定性分析及预测预警方法，包括如下步骤： 

(1)分析牵引式斜坡的变形机理，将斜坡前缘以临界滑面为基准，定义为破坏后区、临界状态面及稳定状态区，牵引式斜坡高度与变形的关系呈现出S曲线特征；分析可能发生破坏的三种情况，情况I为沿最弱的滑面发生破坏，情况II为斜坡发生新的剪切破坏，情况III为坡体或坡体上的节理或裂隙发生拉破坏； 

(2)计算获得滑面各点的临界应变值，并计算出现行滑面临界应力状态点及位移值； 

(3)针对步骤(1)的三种情况分别计算稳定系数，进行稳定性分析； 

(4)利用边坡稳定性计算的滑面边界法获取斜坡坡面位移，进行预测预警。 

所述情况II的剪切破坏由滑面的剪应力与剪应变的全过程曲线决定，情况III的拉破坏由坡体的拉应力与拉应变全过程曲线决定。 

所述步骤(2)的具体计算过程如下： 

步骤(211)分析斜坡的基本形态、特征，试验获得滑面的基本物理力学参数G、S、m、ρ、C、φ、a₁，a₂，a₃，ξ_N和坡体的弹性模量力学参数，计算相对应的位移场和应力场，通过应力场、位移场决定相对应的稳定系数； 

步骤(2.2)将步骤(2.1)得到的参数代入公式τ＝Gγ[1+γ^m/S]^ρ，式中τ、γ分别表示材料的剪应力和剪应变，τ、G的单位为MPa或kPa或Pa，S、m和ρ为无单位参数，且-1<ρ≤0和1+mρ≠0； 

临界应力空间τ^peak采用摩尔库伦准则τ^peak＝C+σ_ntanφ，式中C为凝聚力，σ_n为法向应力，C和σ_n的单位为MPa或kPa或Pa，φ为滑面摩擦角； 

临界应变空间γ_peak采用公式式中σ_n为法向应力，单位为MPa或kPa或Pa； 

临界应力空间与临界应变空间的关系为τ^peak/γ_peak＝G[1-1/(1+mρ)]^ρ，临界应变空间符合公式S+(1+mρ)γ^m_peak＝0； 

参数ρ＝ρ₀/(1+(ρ₀/ρ_c-1)(σ_n/σ_n^c)^ζ)，式中ρ₀为法向应力σ_n为零值的ρ值，ρ_c为σ_n等于σ_n^c时的ρ值，ζ为常系数。 

所述步骤(3)可利用不平衡拉力法计算稳定系数： 

(3.1)先对选定坡体的滑面，结合坡面特征点决定滑面特征点，对滑面特征点的切线作垂线，交滑体坡面形成不同特征点； 

(3.2)连接各特征点形成直线，形成不平衡拉力法的计算条块； 

(3.3)不平衡拉力法的假设同不平衡推力法假设； 

(3.4)计算时，按条块编号从下往上进行； 

(3.5)计算方法为使斜坡体最上面的一个条块推力为零，从而获得不平衡拉力法的各条块下滑力、阻力、剩余拉力和稳定系数； 

上述不平衡拉力法的条块划分是为了决定坡体各特征点对应的最小断面，比较各断面所对应的剩余拉力P_i和滑体相应断面对应的最大抗拉力大小，从而决定各断面是否发生拉破坏，当时，滑体发生拉破坏； 

式中σ_c^t，S_d分别为滑体或节理、裂隙的抗拉强度和对应的断面面积乘积。 

所述步骤(3)还可利用综合下滑力-抗滑力法计算稳定系数： 

对于情况I或情况II选定的滑体的每个条块，对沿滑面的单元的下滑力P_sⁱ和临界摩阻力T_i分别求其矢量和，形成综合下滑力P_s和 综合摩阻力T，且分别与水平轴形成的最小夹角为α_s和α_f，定义稳定性系数为： 

$F_s=\frac{T\cos ({\alpha }_f-{\alpha }_s)}{P_s}---\left(1\right)$

式(1)的物理意义为单位下滑力所产生的抗力； 

对于情况III，假定计算所得e点为临界状态点，选取可能破坏的滑体为研究对象，对于选定滑体的每个条块，对沿滑面的单元的下滑力P_sⁱ和可能破坏面的剩余拉力与临界摩阻力T_i和可能破坏面的最大拉力分别求其矢量和，形成综合下滑力P_s和综合抗滑力T，且分别与水平轴形成的最小夹角为α_s和α_f，定义稳定性系数为式(1)； 

对于情况III，以最小稳定系数作为可能发生拉破坏的稳定系数；在此基础上分析下一步可能的破坏形式。 

所述步骤(3)还可利用牵引式斜坡的主拉力法计算稳定系数： 

对于情况I或情况II选定的滑体的每个条块，m-1条块为临界状态条块，求1～m-1条块滑面的剩余拉力P_sⁱ矢量和，形成综合剩余拉力P^m，且与X轴的夹角为α_s^m；求m至n每个条块的下滑力P_sⁱ与临界摩阻力T_i差值的矢量和，形成综合剩余抗滑力T^m-n，其与水平轴的夹角为α_f^m-n，定义富余稳定性系数为： 

$F_s=\frac{T^{m-n}\cos ({{\alpha }_f}^{m-n}-{{\alpha }_s}^m)}{P^m}---\left(2\right)$

式(2)的物理意义为沿主滑方向单位下滑力的剩余抗力； 

对于情况III，假定计算所得e点为临界状态点，选取可能破坏的滑体为研究对象，m-1条块为临界状态条块，求1～m-1条块滑面的剩余拉力P_sⁱ矢量和，形成综合剩余拉力P^m，且与X轴的夹角为α_s^m；求m条块至bf断面对应的每个条块的下滑力P_sⁱ与临界摩阻力T_i差值，以及bf断面最大抗拉力与此时具有的拉力差值矢量和，形成综合剩余抗滑力T^m-n，其与水平轴的夹角为α_f^m-n，定义富余稳定性系数为式(2)。 

所述步骤(3)还可利用牵引式斜坡变形稳定系数法计算稳定系数： 

定义牵引式边坡现状稳定系数，X轴方向的现状稳定系数为边坡后缘关键块发生破坏时沿滑面在X轴方向的位移矢量和S_p-t除以现状边坡前缘至滑面现状临界状态单元或条块在X轴方向的位移矢量和S_c-t，采用相同的方法定义Y、Z轴方向的现状稳定系数；综合稳定系数定义为边坡后缘关键块发生破坏时沿滑面在X、Y、Z轴方向的位移 矢量和除以现状边坡前缘至滑面现状临界状态单元或条块在X、Y、Z轴方向的位移矢量和；存在三个不同方向的稳定系数，X、Y、Z轴三方向的稳定系数分别为F_s-x＝S_c-t^x/S_p-t^x、F_s-y＝S_c-t^y/S_p-t^y、F_s-z＝S_c-t^z/S_p-t^z。 

步骤(4)中所述边坡稳定性计算的滑面边界法，对于牵引式斜坡，取滑体作为计算对象，潜在滑动面上的力与位移以边界条件加以处理，滑面以梁单元或其它单元加以处理，基本假设为梁单元沿横向没有变形或变形遵循不同本构关系，坡体采用现行的各种单元；具体计算步骤如下： 

第一步，选取坡体作为有限元计算对象，除滑面的正应力与剪应力外，坡体表面边界条件计算方法与常规一致，以滑面底边梁单元为第1单元，沿着滑面底边直至滑面顶面的最后一个梁单元依次排序定义为第N单元； 

第二步，在滑面垂直方向赋应变为零，并在此基础上进行斜坡的第一次计算，从而获得滑面梁单元底面第一次计算的法向应力、切向应力和应变，如式(3，4)， 

${T_{n,s}}^1=\left(\begin{array}{cc}{{\sigma }_n}^{\mathrm{1,1}}& {{\tau }_s}^{\mathrm{1,1}}\\ ...& ...\\ {{\sigma }_n}^{m,1}& {{\tau }_s}^{m,1}\\ ...& ...\\ {{\sigma }_n}^{n,1}& {{\tau }_s}^{n,1}\end{array}\right)---\left(3\right)$${S_{n,s}}^1=\left(\begin{array}{cc}{{ϵ}_n}^{\mathrm{1,1}}=0& {{ϵ}_s}^{\mathrm{1,1}}\\ {{ϵ}_n}^{\mathrm{2,1}}=0& {{ϵ}_s}^{\mathrm{2,1}}\\ ...& ...\\ {{ϵ}_n}^{n-\mathrm{1,1}}=0& {{ϵ}_s}^{n-\mathrm{1,2}}\\ {{ϵ}_n}^{n,2}=0& {{ϵ}_s}^{n,2}\end{array}\right)---\left(4\right)$

式中T_n，sⁱ，S_n，sⁱ分别表示滑面第i次计算所得的法向应力、切向应力和法向应变、切向应变矩阵，σ_n^m，i，τ_s^m，i，ε_n^m，i和ε_s^m，i分别表示第m单元第i次计算所得的法向应力、切向应力、法向应变和切向应变；n表示法向，s表示切向，i表示计算次数； 

比较每一单元第一次计算所得的切向应力与计算所得的法向应力相对应的临界摩阻力的大小及方向，当切向应力与临界摩阻力方向相反时，且切向应力绝对值大于临界摩阻力时，相应单元切向应力作为第二次计算时的边界条件取为，临界摩阻力加上该单元计算所得切向力作为该单元下次计算的初始切向应力边界条件，假定1～k单元切向应力与摩阻力方向相反；对于某单元计算所得切向应力方向与摩阻力方向一致或相反，但绝对值小于临界摩阻力，下次计算时不作为初始条件，假定为k+1～N单元，如式(5)， 

△τ_s^k，1＝τ_crit^k，1+τ_s^k，1，L∈(1，k)   (5) 

第三步，进行第二次计算，相对应1～k及1～N个单元力的初始应力及应变边界条件分别如下式， 

按式(6)和式(7)，进行第二次计算所得的滑面单元的法向应力及切向应力矩阵见式(8)， 

${T_{n,s}}^2=\left(\begin{array}{cc}{{\sigma }_n}^{\mathrm{1,2}}& {{\tau }_s}^{\mathrm{1,1}}+{{\tau }_{\mathrm{crit}}}^{\mathrm{1,1}}\\ ...& ...\\ {{\sigma }_n}^{k,2}& {{\tau }_s}^{k,1}+{{\tau }_{\mathrm{crit}}}^{k,1}\\ {{\sigma }_n}^{k+\mathrm{1,2}}& {{\tau }_s}^{k+\mathrm{1,2}}\\ ...& ...\\ {{\sigma }_n}^{n,2}& {{\tau }_s}^{n,2}\end{array}\right)---\left(8\right)$

首先对1～N单元进行判定，按式(8)的法向应力以下式计算对应的临界摩阻力： 

τ_crit^i，2＝C_i+σ_n^i，2tanφ_i  i∈(1，N)   (9) 

式中，τ_crit^i，j为第i单元第j次计算所得临界摩阻力，σ_n^i，j为第i单元第j次计算所得法向应力； 

对于已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时计算的法向应力之差大于某定值D₁时，重新对该单元切向应力按式(11)赋值， 

|σ_n^L，2-σ_n^L，1|≥D₁，L∈(1，k)   (10) 

△τ_s^L，2＝τ_s^L，1+τ_crit^L，1+(τ_crit^L，2-τ_crit^L，1)/2， 

L∈(1，k)   (11) 

对于未赋临界摩阻力的单元，当切向应力与临界摩阻力方向相反时，且绝对值大于临界摩阻力时，对这样的单元施加临界摩阻力如下式(12)， 

如第i单元第2次计算，施加的切向应力((τ_crit^i，2+τ_s^i，2)，假定i∈(k，m))； 

△τ_s^i，2＝τ_crit^i，2+τ_s^i，2，i∈(k，m)   (12) 

按照上述计算，在进行第三次计算时，1～N个滑面应力及应变边界条件分别为： 

第四步：重复第三步，计算直至多次，对于已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时与上一次计算的法向应力之差最大值小于某定值D₁时，则可以结束计算， 

|σ_n^L，MM+1-σ_n^L，MM|_max<D₁，L∈(1，m)   (15) 

式中，MM为计算次数，m为施加临界摩阻力单元数； 

此时记下法向应力σ_n^i，MM+1和切向应力τ_s^i，MM+1，其中i∈(1，N)，以及切向应变ε_s^i，MM+1，其中i∈(1，N)。 

在考虑渐进破坏，边坡的临界状态稳定性系数计算，则施加的临界应力场应力可以利用申请人提出的“基于边坡变形破坏机理的临界位移预测预警方法”的临界应力。便可以获得考虑渐进破坏边坡临界状态稳定性系数。 

牵引式斜坡有限元的稳定性计算结果，可以利用方法一(不平衡拉力法)、方法二(综合下滑力-抗滑力法)、方法三(主拉力法)和方法四(变形稳定系数法)进行牵引式斜坡的稳定性评价，使斜坡预测预报与斜坡变形紧密相关。 

本发明基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法的有益效果是： 

1)提出的牵引式边坡不平衡拉力法可以计算临界状态和渐进破坏状态的稳定系数； 

2)可以通过计算等决定牵引式边坡的破坏形式、临界状态点或面，临界状态和渐进破坏状态的综合下滑力-摩阻力稳定系数、主拉力富余系数和基于变形的稳定系数；决定牵引式斜坡破坏时滑面、滑体及坡面上不同点的变形值； 

3)可以实施牵引式边坡渐渐破坏变形及力的全过程动态描述； 

4)结合变形监测，可以实施边坡的防护措施的可靠性评价。 

附图说明

图1a为斜坡剪应力(τ)～剪应变(γ)全过程曲线； 

图1b为牵引式斜坡变形机理破坏图； 

图1c为坡体拉应力(σ)～应变(ε)全过程曲线； 

图2(a)为滑面的剪应力与剪应变的全过程曲线； 

图2(b)为牵引式斜坡高度与变形的关系曲线； 

图2(c)为牵引式斜坡变形机理破坏图； 

图2(d)为坡体的拉应力与拉应变全过程曲线； 

图2(e)为牵引式斜坡的剪切破坏变形随时间的关系曲线； 

图3为不平衡拉力法计算框图； 

图4为本发明方法的边坡有限元计算单元划分图； 

图5为本发明方法的力边界处理图； 

图6为本发明方法的位移边界处理图。 

图中，τ为剪应力，γ为剪应变，σ为拉应力，ε为拉应变，I、II、III为破坏三种形式，S为位移，t为时间，H为高度，T_sⁱ为第i单元阻力，T_s，critⁱ为第i单元临界阻力，U_sⁱ为第i单元位移。 

具体实施方式

首先分析牵引式斜坡的变形机理，从而决定斜坡前缘为破坏后区(如图1中a₁′～a₃′段，)、临界状态面(如图1中a₃′点)及稳定状态区(如图1中a₃′～a₆′段)，提出牵引式斜坡高度与变形的关系呈现出“S”曲线特征(如图2中的b图)，分析可能破坏的三种情况(情况I：沿最弱的滑面发生破坏(如图1中的a₁′～a₆′段)，情况II：斜坡产生新的剪切破坏(如图1中的a₁′～b₆′段)，情况III：坡体(或坡体上的节理或裂隙)发生拉破坏(如图1中的a₃′～c₃段)。指出剪切破坏由滑面的剪应力与剪应变的全过程曲线决定(如图2中的a图)，拉破坏由坡体的拉应力与拉应变全过程曲线决定(如图2中的d图)。牵引式斜坡的剪切破坏变形随时间的关系遵循图2中的e图，拉破坏遵循图2中的d图，即拉破坏时，时间很短。由此可见，牵引式斜坡破坏时，变形很小，时间也很短。 

计算分析如下： 

本构方程及参数决定： 

步骤1、分析斜坡的基本形态、特征，试验获得滑面的基本物理力学参数G、S、m、ρ、C、φ、a1，a₂a₃，ξ_N和坡体的弹性模量等力学参数，计算相对应的位移场和应力场，通过应力场、位移场决定相对应的稳定系数； 

步骤2、将步骤1得到的参数代入公式τ＝Gγ[1+γ^m/S]^ρ，式中τ、γ分别表示材料的剪应力和剪应变，τ、G的单位为MPa或kPa或Pa，S、m和ρ为无单位参数，且-1<ρ≤0和1+mρ≠0； 

临界应力空间τ^peak采用摩尔库伦准则τ^peak＝C+σ_ntanφ，式中C为凝聚力，σ_n为法向应力，C和σ_n的单位为MPa或kPa或Pa，φ为滑面摩擦角； 

临界应变空间γ_peak采用公式式中σ_n为法向应力，单位为MPa或kPa或Pa； 

临界应力空间与临界应变空间的关系为τ^peak/γ_peak＝G[1-1/(1+mmρ)]^ρ，临界应变空间符合公式S+(1+mρ)γ^m_peak＝0； 

参数ρ＝ρ₀/(1+(ρ₀/ρ_c-1)(σ_n/σ_n^c)^ζ)，式中ρ₀为法向应力σ_n为零值的ρ值，ρ_c为σ_n等于σ_n^c时的ρ值，ζ为常系数； 

步骤3、经步骤2可以获得滑面各点的临界应变值，并可以计算出现行滑面临界应力状态点及位移值。 

稳定性分析： 

方法一：不平衡拉力法 

a)先对选定坡体的滑面，结合坡面特征点决定滑面特征点，对滑面特征点的切线作垂线，交滑体坡面形成不同特征点，b)连接各特征点形成直线，由(a)和(b)两步形成不平衡拉力法的计算条块，c)不平衡拉力法的假设同不平衡推力法假设，d)计算时，按条块编号从下往上进行，如图(3)，e)计算方法与不平衡推力法相同，即使最后一个条块(斜坡体最上面的一个条块)推力为零，从而获得不平衡拉力法的各条块下滑力、阻力、剩余拉力和稳定系数。 

滑体拉破坏的决定： 

上述不平衡拉力法的条块划分就是为了决定坡体各特征点对应的最小断面，比较各断面所对应的剩余拉力(P_i)和滑体相应断面对应的最大抗拉力(S_d：分别为滑体(或节理、裂隙等)抗拉强度和对应的断面面积乘积)大小，从而决定各断面是否发生拉破坏，当时，滑体发生拉破坏。 

方法二：综合下滑力-抗滑力法 

对于选定的滑体(情况I：沿最弱的滑面发生破坏，情况II：斜坡产生新的剪切破坏，)的每个条块，对沿滑面的单元的下滑力(P_sⁱ)和临界摩阻力(T_i)分别求其矢量和，形成综合下滑力P_s和综合摩阻力T，且分别与水平轴形成的最小夹角为α_s和α_f，定义稳定性系数为： 

$F_s=\frac{T\cos ({\alpha }_f-{\alpha }_s)}{P_s}---\left(1\right)$

式(1)的物理意义为：单位下滑力所产生的抗力。 

对于情况III(滑体发生拉破坏)，以图3为例，假如计算所得e点为临界状态点，选取可能破坏的滑体为研究对象(如图3中的bcdef滑体)，对于选定滑体的每个条块，对沿滑面的单元的下滑力(P_sⁱ)和可能破坏面(如图3中的fb面)的剩余拉力与临界摩阻力(T_i)和可能破坏面(如图3中的fb面)的最大拉力(P_i^t)分别求其矢量和，形成综合下滑力P_s和综合抗滑力T，且分别与水平轴形成的最小夹角为α_s和α_f，定义稳定性系数为式(1)。对于情况III，可以求其多种可能的破坏形式，以最小稳定系数作为可能发生拉破坏的稳定系数。在此基础上分析下一步可能的破坏形式，亦即斜坡破坏存在初始发生拉破坏，随后发生剪破坏。 

方法三：牵引式斜坡的主拉力法 

对于选定的滑体(情况I：沿最弱的滑面发生破坏，情况II：斜坡产生新的剪切破坏，)的每个条块，求1～m-1条块滑面(m-1条块为临界状态条块)的剩余拉力(P_sⁱ)矢量和，形成综合剩余拉力P^m，且与X轴的夹角为α_s^m；求m至n每个条块的下滑力(P_sⁱ)与临界摩阻力(T_i)差值的矢量和，形成综合剩余抗滑力T^m-n，其与水平轴的夹角为α_f^m-n，定义富余稳定性系数为： 

$F_s=\frac{T^{m-n}\cos ({{\alpha }_f}^{m-n}-{{\alpha }_s}^m)}{P^m}---\left(2\right)$

式(2)的物理意义为：沿主滑方向单位下滑力的剩余抗力。 

对于情况III(滑体发生拉破坏)，以图3为例，假如计算所得e点为临界状态点，选取可能破坏的滑体为研究对象(如图3中的bcdef滑体)，求1～m-1条块滑面(m-1条块为临界状态条块)的剩余拉力(P_sⁱ)矢量和，形成综合剩余拉力P^m，且与X轴的夹角为α_s^m；求m条块至bf断面对应的每个条块的下滑力(P_sⁱ)与临界摩阻力(T_i)差值，以及bf断面最大抗拉力与此时具有的拉力差值矢量和，形成综合剩余抗滑力T^m-n，其与水平轴的夹角为α_f^m-n，定义富余稳定性系数为(2)式。 

方法四：牵引式斜坡变形稳定系数法 

定义牵引式边坡现状稳定系数：X轴方向的现状稳定系数为：边坡后缘关键块发生破坏时沿滑面在X轴方向的位移矢量和(S_p-t)除以现状边坡前缘至滑面现状临界状态单元(或条块)在X轴方向的位移矢量和(S_c-t)，同理采用相同的方法定义Y、Z轴方向的现状稳定系数；综合稳定系数定义为：边坡后缘关键块发生破坏时沿滑面在X、Y、Z轴方向的位移矢量和除以现状边坡前缘至滑面现状临界状态单元(或条块)在X、Y、Z轴方向的位移矢量和。存在三个不同方向的稳定系数，X、Y、Z轴三方向的稳定系数分别为F_s-x＝S_c-t^x/S_p-t^x、F_s-y＝S_c-t^y/S_p-t^y、F_s-z＝S_c-t^z/S_p-t^z。 

牵引式斜坡有限元滑面边界计算法： 

对于牵引式斜坡，取滑体作为计算对象，潜在滑动面上的力与位移以边界条件加以处理，滑面以梁单元(或其它单元)加以处理，基本假设为：梁单元沿横向没有变形或变形遵循不同本构关系，坡体可以采用现行的各种单元，见图4。 

牵引式斜坡稳定的详尽计算方法，具体计算步骤如下： 

选取坡体作为有限元计算对象，坡体表面边界条件(除滑面的正应力与剪应力外)计算方法与常规一致，以滑面底边梁单元为第1单 元，沿着滑面底边直至滑面顶面的最后一个梁单元依次排序定义为第N单元(见图4)。 

第一步：鉴于滑面垂直方向位移较小，因此在滑面垂直方向赋应变为零，并以此加以计算； 

第二步：在上述边界条件的基础上，进行斜坡的第一次计算，从而获得滑面梁单元底面第一次计算的法向应力、切向应力和应变(如下式(3，4))。 

${T_{n,s}}^1=\left(\begin{array}{cc}{{\sigma }_n}^{\mathrm{1,1}}& {{\tau }_s}^{\mathrm{1,1}}\\ ...& ...\\ {{\sigma }_n}^{m,1}& {{\tau }_s}^{m,1}\\ ...& ...\\ {{\sigma }_n}^{n,1}& {{\tau }_s}^{n,1}\end{array}\right)---\left(3\right)$${S_{n,s}}^1=\left(\begin{array}{cc}{{ϵ}_n}^{\mathrm{1,1}}=0& {{ϵ}_s}^{\mathrm{1,1}}\\ {{ϵ}_n}^{\mathrm{2,1}}=0& {{ϵ}_s}^{\mathrm{2,1}}\\ ...& ...\\ {{ϵ}_n}^{n-\mathrm{1,1}}=0& {{ϵ}_s}^{n-\mathrm{1,2}}\\ {{ϵ}_n}^{n,2}=0& {{ϵ}_s}^{n,2}\end{array}\right)---\left(4\right)$

式中：T_n，sⁱ，S_n，sⁱ，：分别表示滑面第i次计算所得的法向应力、切向应力和法向应变切向应变矩阵，σ_n^m，i，τ_s^m，i，ε_n^m，i和ε_s^m，i分别表示第m单元第i次计算所得的法向应力、切向应力、法向应变和切向应变；n：表示法向，s：表示切向，i：表示计算次数。 

为了获得滑面的边界条件，比较每一单元第1次计算所得的切向应力(如：第k单元切向应力τ_s^k，1)与计算所得的法向应力相对应的临界摩阻力(如：第k单元临界摩阻力τ_crit^k，1，计算式为：τ_crit^k，1＝C_k+σ_n^k，1tanφ_k(也可以采用其它的应力准则))的大小及方向， 

当切向应力与临界摩阻力方向相反时，且切向应力绝对值大于临界摩阻力时，相应单元切向应力作为第二次计算时的边界条件取为：临界摩阻力加上该单元计算所得切向力(如方程(5))作为该单元下次计算的初始切向应力边界条件(假如，1～k单元切向应力与摩阻力方向相反)，对于某单元计算所得切向应力方向与摩阻力方向一致或相反，但绝对值小于临界摩阻力，下次计算时不作为初始条件(假如为：k+1～N单元)。 

△τ_s^k，1＝τ_crit^k，1+τ_s^k，1，L∈(1，k)   (5) 

第三步：进行第二次计算。相对应1～k及1～N个单元力的初始应力及应变边界条件分别如下式： 

按上述(6)和(7)式，进行第二次计算所得的滑面单元的法向 应力及切向应力矩阵见式(8)。 

${T_{n,s}}^2=\left(\begin{array}{cc}{{\sigma }_n}^{\mathrm{1,2}}& {{\tau }_s}^{\mathrm{1,1}}+{{\tau }_{\mathrm{crit}}}^{\mathrm{1,1}}\\ ...& ...\\ {{\sigma }_n}^{k,2}& {{\tau }_s}^{k,1}+{{\tau }_{\mathrm{crit}}}^{k,1}\\ {{\sigma }_n}^{k+\mathrm{1,2}}& {{\tau }_s}^{k+\mathrm{1,2}}\\ ...& ...\\ {{\sigma }_n}^{n,2}& {{\tau }_s}^{n,2}\end{array}\right)---\left(8\right)$

对于式(8)的计算结果，分析如下： 

首先对1～N单元进行判定，按(8)式的法向应力以下式计算对应的临界摩阻力： 

τ_crit^i，2＝C_i+σ_n^i，2tanφ_i  i∈(1，N)   (9) 

式中：τ_crit^i，j：第i单元，第j次计算所得临界摩阻力，σ_n^i，j：第i单元，第j次计算所得法向应力。 

分析切向应力，对于那些已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时计算的法向应力之差大于某定值时(如：式(10)的D₁值)，重新对该单元切向应力按式(11)赋值。 

|σ_n^L，2-σ_n^L，1|≥D₁，L∈(1，k)   (10) 

△τ_s^L，2＝τ_s^L，1+τ_crit^L，1+(τ_crit^L，2-τ_crit^L，1)/2， 

L∈(1，k)    (11) 

对于那些未赋临界摩阻力的单元，当切向应力与临界摩阻力方向相反时，且绝对值大于临界摩阻力时，对这样的单元施加临界摩阻力如下式(12)(如：第i单元第2次计算，施加的切向应力((τ_crit^i，2+τ_s^L，2)，假如：i∈(k，m))； 

△τ_s^i，2＝τ_crit^i，2+τ_s^i，2，i∈(k，m)   (12) 

按照上述计算，在进行第三次计算时，1～N个滑面应力及应变边界条件分别为： 

第四步：重复第三步计算，直至多次，对于那些已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时与上一次计算的法向应力之差最大值小于某定值时(如：式(15))，则可以结束计算。 

|σ_n^L，MM+1-σ_n^L，MM|_max<D₁，L∈(1，m)   (15) 

MM：计算次数，m：施加临界摩阻力单元数。 

此时记下法向应力和切向应力(σ_n^i，MM+1，τ_s^i，MM+1，其中i∈(1，N))，以及切向应变(ε_s^i，MM+1)，其中i∈(1，N))。 

临界状态单元的决定：上述计算过程中，存在一点(或一个单元)其下滑力正好等于临界摩阻力，这一点(或一个单元)称为临界状态单元，然而，在计算过程中，由于是以面代点(或以体代点)，该点的计算往往存在一定的剩余拉力，为了更准确地决定临界状态单元，可以将该单元进行分割为更小的单元，当切向应力绝对值小于D₂分之的临界摩阻力，则认为可以满足要求；当大于D₂分之的临界摩阻力，则将此单元进行分割成二个及以上单元，进行上述计算，直至该单元切向应力绝对值小于D₂分之的临界摩阻力，则可近似认为此单元为临界状态单元。 

上述计算所得的滑体应力场及滑面应力场为“视真实应力场”(因为滑面真实的应力：部分应为残余应力、部分应为破坏后区应力和部分应力位于峰值应力前状态，而本文赋值应力场均为临界状态应力，所以只能称之为“视真实应力场”)，整个计算边界条件赋值见图5。 

分解1～m滑面单元的临界摩阻应力、下滑应力和剩余拉应力分别为： 

临界摩阻应力：τ_s，crit^i，MM+1：τ_s，crit^i，MM+1＝C_i+σ_n^i，MM+1tanφ_i   (16) 

下滑应力：τ_s^i，MM+1    (17) 

剩余拉应力：τ_s，p^i，MM+1＝τ_s^i，MM+1-τ_s，crit^i，MM+1    (18) 

分解m+1～n滑面单元的临界摩阻应力、下滑应力和剩余拉应力分别为： 

临界摩阻应力：τ_s，crit^i，MM+1：τ_s，crit^i，MM+1＝C_i+σ_n^i，MM+1tanφ_i   (19) 

下滑应力：τ_s^i，MM+1，  (20) 

剩余拉应力：0。   (21) 

对于渐进变化过程的真实应力场模拟，其施加的滑面应力场和应变场可以利用申请人提出的“基于边坡变形破坏机理的临界位移预测预警方法”的应力场和位移场，即可以利用满足该本构关系的滑面单元决定滑面边界条件，从而获得真实的应力场和位移场。滑面边界的位移处理如图6。 

提出的牵引式滑面边界法计算稳定系数具有可比性。 

利用此方法，计算所得稳定系数可以和传统条分法加以比较，其稳定系数计算步骤如下： 

以有限单元法为基础，采用传统强度折减法计算出的稳定系数应能与本发明提出的不平衡拉力法求得的结果进行比较，比较两者结果的基础是：两种计算方法的计算实质必须基本一致。不平衡拉力条分法的计算实质是：假定条间力的合力与上一土条底面相平行，根据力 的平衡条件，逐渐变化稳定系数(f)，并逐条向上推求，直至最后一土条块的拉力零。基于这种计算思路，提出有限单元滑面边界法的计算方法如下：上述滑面边界法获得相对应于滑面的滑体单元应力与应变矩阵分别为： 

${T_{n,s}}^{B,\mathrm{MM}+1}=\left(\begin{array}{cc}{{\sigma }_n}^{1,\mathrm{MM}+1}& {{\tau }_s}^{1,\mathrm{MM}+1}\\ ...& ...\\ ...& {{\tau }_s}^{m,\mathrm{MM}+1}\\ ...& {{\tau }_s}^{m+1,\mathrm{MM}+1}\\ ...& ...\\ {{\sigma }_n}^{n,\mathrm{MM}+1}& {{\tau }_s}^{n,\mathrm{MM}+1}\end{array}\right)---\left(22\right)$${S_{n,s}}^{B,\mathrm{MM}+1}=\left(\begin{array}{cc}{{ϵ}_n}^{1,\mathrm{MM}+1}=0& {{ϵ}_s}^{1,\mathrm{MM}+1}\\ {{ϵ}_n}^{2,\mathrm{MM}+1}=0& {{ϵ}_s}^{2,\mathrm{MM}+1}\\ ...& ...\\ {{ϵ}_n}^{n-1,\mathrm{MM}+1}=0& {{ϵ}_s}^{n-1,\mathrm{MM}+1}\\ {{ϵ}_n}^{n,\mathrm{MM}+1}=0& {{ϵ}_s}^{n,\mathrm{MM}+1}\end{array}\right)---\left(23\right)$

滑面边界法强度折减计算步骤如下： 

第一步：初次给定计算初始稳定系数(f₁)，在整个强度折减计算中，第一个滑面单元法向应力是不被改变的，即第一个滑面单元法向应力σ_n^1，MM+1永远不变，当然临界摩阻力也是不变的，以下式加以表不： 

τ_s，crit^1，MM+1＝C₁+σ_n^1，MM+1tanφ₁  (24) 

判定滑面第一单元拉力是否与下滑力方向相反，如相反，则滑面下滑应力恒定(τ_s^1，MM+1)，滑面第一单元剩余拉应力表示为： 

τ_s，p^1，MM+1＝τ_s^1，MM+1-τ_s，crit^1，MM+1/f₁，(25) 

而施加的切向应力场改变为： 

${{\tau }_{s,p}}^{1,{\mathrm{NN}}_0+1}={{\tau }_s}^{1,\mathrm{MM}+1}+\frac{f_1-1}{f_1}{{\tau }_{s,\mathrm{crit}}}^{1,\mathrm{MM}+1}---\left(26\right)$

即进行强度折减第一步计算的初始滑面边界条件为： 

按上述(27)(28)计算所得滑面边界值为： 

${T_{n,s}}^{S-B,1-j}=\left(\begin{array}{cc}{{\sigma }_n}^{1,1}={{\sigma }_n}^{1,\mathrm{MM}+1}& {{\tau }_s}^{\mathrm{1,1}}={{\tau }_s}^{1,{\mathrm{NM}}_0+1}\\ {{\sigma }_n}^{\mathrm{2,1}}={{\sigma }_n}^{2,N_0}& {{\tau }_s}^{\mathrm{2,1}}={{\tau }_s}^{2,{\mathrm{NM}}_0+1}\\ ...& ...\\ ...& ...\\ ...& ...\\ {{\sigma }_n}^{n,1}={{\sigma }_n}^{n,N_0}& {{\tau }_s}^{n,1}={{\tau }_s}^{n,{\mathrm{NM}}_0+1}\end{array}\right)---\left(29\right)$${S_{n,s}}^{S-B,1-j}=\left(\begin{array}{cc}{{ϵ}_n}^{\mathrm{1,1}}=0& {{ϵ}_s}^{\mathrm{1,1}}={{ϵ}_s}^{1,{\mathrm{NM}}_0+1}\\ {{ϵ}_n}^{\mathrm{2,1}}=0& {{ϵ}_s}^{\mathrm{2,1}}={{ϵ}_s}^{2,{\mathrm{NM}}_0+1}\\ ...& ...\\ {{ϵ}_n}^{n-\mathrm{1,1}}=0& {{ϵ}_s}^{n-\mathrm{1,1}}={{ϵ}_s}^{n-1,{\mathrm{NM}}_0+1}\\ {{ϵ}_n}^{n,1}=0& {{ϵ}_s}^{n,1}={{ϵ}_s}^{n,{\mathrm{NM}}_0+1}\end{array}\right)---\left(30\right)$

第二步计算初始边界条件只对滑向第二个边界单元进行改变，判定滑面第二单元拉力是否与摩阻力方向相反，如相反，则保持滑面第二单元法向应力场不变，修改切向应力场为： 

临界摩阻应力： 

${{\tau }_{s,\mathrm{crit}}}^{2,{\mathrm{NN}}_0+1}=C_2+{{\sigma }_n}^{2,{\mathrm{NN}}_0+1}{\tan \phi }_2---\left(31\right)$

而施加的切向应力场为： 

${{\tau }_{s,p}}^{2,{\mathrm{NN}}_0+1}={{\tau }_s}^{2,{\mathrm{NN}}_0+1}-\frac{1}{f_1}{{\tau }_{s,\mathrm{crit}}}^{2,{\mathrm{NN}}_0+1}---\left(32\right)$

第二步计算的初始滑面边界应力条件为： 

依次按上述步骤进行，直到滑面第n单元的切向应力正好等于该单元f₁分之一的临界摩阻应力(或之差绝对值小于一个定值)，则停止计算，此时的f₁即为稳定系数；否则重新增大或减小f₁值，直至滑面第n单元的切向应力正好等于该单元f_NM分之一的临界摩阻应力(或之差绝对值小于一个定值)(如：34式D₃值)。 

|τ_s^n，NM-τ_s，crit^n，NM/f_NM|<D₃   (34) 

式中：NM、f_NM：分别表示迭代次数和强度折减稳定系数。