CN106503421A - 一种边坡渐进破坏全过程计算方法 - Google Patents

Abstract

本发明提出了一种边坡渐进破坏全过程计算方法，在假设地质材料全过程主应力‑主应变满足相关本构关系条件下，计算材料全过程主应力‑主应变的任意应力和应变，利用坐标转换与可逆应变和不可逆应变无关的特点，获得任意应力状态下的应力和应变值；在滑面破坏满足某一峰值应力准则条件下，由于破坏区摩阻力的下降，产生剩余下滑力，增加滑体前面单元的法向应力和沿滑面的驱动剪应力，致使边坡一步一步向前发生破坏；通过计算滑面上滑体单元应变，决定滑体与滑床的相对应变值，利用该应变值计算滑面摩阻力。本发明直接计算滑面剪应力和剪应变，不存在条块之间变形协调假设，也可以计算滑体的应力和应变。

Description

一种边坡渐进破坏全过程计算方法

技术领域

本发明涉及土木工程和地质灾害、基坑、锚杆锚索与抗滑桩设计等的防治、评价和预测预报等技术领域，特别涉及一种边坡渐进破坏全过程计算方法。

背景技术

边坡的稳定性分析和渐进破坏计算大多数采用的是理想弹塑性模型，这种模型在应力达到峰值应力时，应变无数个，另外摩阻力也恒定等，该模型特征决定了难以描述边坡的渐进破坏过程的。另外，边坡沿滑面在破坏后区，其剪应力和剪应变不连续，传统的数值方法难以较好描述其破坏特征，急需建立一种描述边坡渐进破坏特征的数值计算方法。

发明内容

本发明所要解决的技术问题是提供一种边坡渐进破坏全过程计算方法，不仅可以计算边坡的渐进全过程变形破坏特征，同时可以用于决定边坡的潜在滑动面。

为解决上述技术问题，本发明采用的技术方案是：

一种边坡渐进破坏全过程计算方法，包括以下步骤：

当滑体和滑面为同种材料时，

步骤1：对材料实施主应力-应变全过程曲线试验，试验获得峰值主应力、应变和全过程曲线特征，计算滑体材料的物理力学特征方程，包括主应力-主应变方程、峰值应力方程、临界主应变方程、弹性模量变化方程、软化系数方程；

步骤1.1：主应力-主应变方程

主应力-主应变为四参数本构方程：

式中，σ_i、ε_i分别为主应力和主应变，E_ii为初始弹性模量，m_i、S_i、ρ_i为主应力方向的常系数，σ_i、E_ii的单位为MPa或kPa或Pa，m_i、S_i、ρ_i为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ρ_i≤0和1+m_iS_i≠0；临界应变空间满足关系式：

式中，ε_i ^peak为临界主应力对应的主应变；

步骤1.2：假设临界主应力满足摩尔库伦准则，即

式中，C为凝聚力，C和σ_i的单位为MPa或kPa或Pa，φ为材料摩擦角；

步骤1.3：设临界主应变空间仅相关于第三主应力，临界主应变ε_i ^peak采用关系式：

或

式中，a_i,1,a_i,2,a_i,3,ζ_i,N,为常系数；a_i,1,a_i,2单位为MPa或kPa或Pa，a_i,3,ζ_i,N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²；

步骤1.4：初始弹性模量方程：

E_ii＝E_ii,0+b_i,1σ₃+b_i,2σ₃ ²

式中，E_ii,0为第三主应力σ₃为零值的E_ii值，b_i,1,b_i,2为常系数，单位分别为无量纲和1/MPa或1/kPa或1/Pa；

步骤1.5：对无量纲参数ρ_i，软化系数演化方程表示为：

式中，ρ_i,0为第三主应力σ₃为零值的ρ_i值，ρ_i,c为σ_i,3等于σ_i,3 ^c时的ρ_i值，为常系数；

步骤1.6：计算滑体主应变在滑面单元上产生的应变ε_ij'表达式为：

主应变ε₁的方向与x轴间所夹角度按式子和求得，当γ_xy为正直，且ε_xx>ε_yy时，角为负值；计算滑面上的单元的应力，计算公式为：

步骤2：在峰值主应力满足摩尔库伦准则条件下，其破坏面上剪应力与最小主应力的夹角为45°+φ/2，φ为材料摩擦角，单元剪应力摩尔库伦准则为：

τ_peak＝C+σ_n tanφ

其中，τ_peak,σ_n,C分别为破坏面最大剪应力、法向应力和凝聚力，单位为MPa或kPa或Pa；

步骤3：滑面上的剪应力和剪应变满足的本构方程，具体如下：

步骤3.1：剪应力-剪应变为四参数本构方程：

τ＝Gγ[1+γ^q/p]^ξ

式中，τ、γ分别为剪应力和剪应变，G为剪切模量，p、q、ξ为在不同法向应力下的常系数，τ、G的单位为MPa或kPa或Pa，p、q、ξ为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ξ≤0和1+qξ≠0；临界应变空间满足关系式：

p+(1+qξ)γ^q _peak＝0

式中，γ_peak为临界应力对应的应变；

步骤3.2：滑体主应变在滑面上产生的法向应变ε_n'和剪应变γ'表达式：

γ'＝(ε₁-ε₃)sin2α

式中，α为滑面与最小主应力方向的夹角；

步骤3.3：初始剪切模量方程：

式中，ν为泊松比；

步骤3.4：对于无量纲参数ξ，软化系数演化方程表示为：

ξ＝ξ₀/(1+(ξ₀/ξ_c-1)(σ_n/σ_n ^c)^ψ)

ξ₀为法向应力σ_n为零值的ξ值，ξ_c为σ_n等于σ_n ^c时的ξ值，ψ为常系数；

对于同一种地质材料，直接利用单元在滑面上求解的驱动剪应力、剪应变和法向应力求解软化系数：求解方法为：

利用破坏面上所得到的法向应力，对应地求解出主应力σ₁,σ₃，方程为：

滑面上的法向应力：

临界剪应力方程：

求解得临界主应力和围压：

利用方程 获得滑面上临界剪应变、临界剪应力和对应的法向应力，计算常系数：

p＝-(γ_peak)^q(1+qξ)

再利用滑面上的实际计算剪应力，以方程τ＝Gγ[1+γ^q/p]^ξ计算软化系数ξ；

若滑面材料与滑体不一致，则滑面计算方法为：

步骤4：对于不同于滑体材料滑面的剪应力-剪应变四参数本构方程、临界应力方程、临界剪应变方程和软化系数演化方程均与方程：

τ＝Gγ[1+γ^q/p]^ξ

τ_peak＝C+σ_ntanφ

或者二者之一

ξ＝ξ₀/(1+(ξ₀/ξ_c-1)(σ_n/σ_n ^c)^ψ)一致，其初始剪切模量方程为：

G＝G₀+b₁σ_n+b₂σ_n ²

式中，G₀为法向应力σ_n为零值的G值，b₁,b₂为常系数，单位为无量纲和1/MPa或1/kPa或1/Pa；

步骤5：沿滑面的变形分布特征

当M条块或单元为临界状态时，沿滑面的力分布特征为：在M+1至N条块间，下滑力等于摩阻力，压应力等于反压应力，在这个阶段，沿滑面的剪应变γ以下滑应力等于抗滑剪应力加以计算；在M条块或单元下滑应力等于抗滑临界摩阻应力，以此计算该条块的临界剪应变γ；

在第1至M-1条块或单元间，下滑应力大于摩阻应力，其任一条块或单元L存在两个位移：第一种位移为第L条块或单元本身产生的位移，即：由第L-1的右边和第L+1的左边位移之差加以计算，第二种位移为第L条块或单元本身产生的位移加上第L条块或单元的底边中心相对该点沿滑床原始位置滑动的位移；第一种位移用于计算滑体的驱动下滑剪应力，第二种位移用于计算滑床提供给滑体的摩阻应力；

L条块或单元的第一种底边位移S_L,j ¹计算：

S_L,j ¹＝ε_j ^Ll^L，j∈(x,y,z)

L条块或单元的第二种底边位移S_L,j ²计算：

式中，ε_j ⁱ,lⁱ分别为沿滑面第i条块或单元在x、y、z轴方向的应变和底边长度；在滑面方向的位移为三者的矢量和，滑面剪应变γ为上述第一、二种位移除以滑面长度；

步骤6：硬化特征描述

当滑体仅表现出硬化特征时，其计算采用简化的模型，将本构方程取ξ＝-1，q＝1，则a,＝1/(Ga”)，b'＝1/(Gp)，在主应力条件下其方程形式与邓肯—张模型一致，此时描述材料的硬化行为特征：

式中，a_i’、b_i’为常系数；

在单元主应力满足方程条件下，相对应的单元应变定义为峰值应变，满足方程或在峰值应力和应变下，方程变为：

定义单元割线模量则

对方程求导，相应的导数为切线模量，在任一应力状态条件下，切向模量E_i表示为：

得到初始弹性模量为1/a_i'，利用方程在主应力满足方程时的切线模量定义为E_t,i，则有：

E_t,i＝a_i’K_scant,i ²

在峰值应力时，研究试验曲线的切向模量为E_t,i，设其具有如下特征：

式中，α_i,k_i为常系数，p_a为大气压；

方程的常系数利用三个不同法向应力试验加以确定；

步骤7：针对条分法，边坡渐进破坏计算方法包括如下步骤：

步骤7.1：利用现行方法进行条块划分；

步骤7.2：竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力满足应力分布条件；计算相对应的应力，同时计算主应力和旋转角，根据主应力计算对应主应变，通过坐标旋转计算滑动方向和其它方向的应变和位移；

步骤7.3：根据现场实测或假设确定第一条块水平和竖直方向的位移，该位移加上步骤7.2计算的位移，所得位移除以滑面长度即为滑面剪应变或视剪应变，利用该剪应变计算条块底边摩阻应力；

步骤7.4：重复步骤7.2和步骤7.3，直至现场判断的临界状态条块为止；

步骤8：实施边坡渐进破坏的数值计算。

与现有技术相比，本发明的有益效果是：能展现边坡渐进变形破坏过程，解决了破坏后区剪应力和剪应变不连续的难题，能使边坡渐进破坏以手算(excel表格)的形式表现出来，可以表述拉、拉剪和压剪破坏特征，既适宜于推移式、也适宜于牵引式，还适宜于两者的混合式。

附图说明

图1是本发明中推移式边坡渐进破坏位移分布图，S1、SL分别为第1、L条块或单元底边中点(现状b1，bL)相对滑坡初始状态中点(原始状态a1，aL)的位移，M为临界状态条块或单元。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细的说明。本发明提供的一种边坡渐进破坏全过程计算方法，包括如下步骤：

1、当滑体和滑面为同种材料时，对材料实施主应力-应变全过程曲线试验，试验获得峰值主应力、应变和全过程曲线特征。例如，可采用中国专利一种边坡渐进破坏潜在滑动面的计算方法(申请号：201510658880.6)决定滑体材料的物理力学特征方程：1)主应力-主应变方程，2)峰值应力方程，3)临界主应变方程，4)弹性模量变化方程，5)软化系数方程)；

1.1、主应力-主应变方程

主应力--主应变为四参数本构方程：

式中：σ_i、ε_i分别为主应力和应变，E_ii为初始弹性模量，m_i、S_i、ρ_i为主应力方向的常系数，σ_i、E_ii的单位为MPa或kPa或Pa，m_i、S_i、ρ_i为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ρ_i≤0和1+m_iS_i≠0；临界应变空间满足如下关系式：

式中：ε_i ^peak为临界主应力对应的主应变。

1.2、假设临界主应力满足摩尔库伦准则(也可以其它准则)：

式中：C为凝聚力，C和σ_i的单位为MPa或kPa或Pa，φ为材料摩擦角。

1.3、假设临界主应变空间仅相关于第三主应力，临界主应变ε_i ^peak采用如下关系式：

或

式中：a_i,1,a_i,2,a_i,3,ζ_i,N,为常系数；a_i,1,a_i,2单位为MPa或kPa或Pa，a_i,3,ζ_i,N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²。

1.4、初始弹性模量方程：

E_ii＝E_ii,0+b_i,1σ₃+b_i,2σ₃ ² (1.5)

式中：E_ii，0为第三主应力σ₃为零值的E_ii值，b_i，1,b_i，2为常系数，单位分别为无量纲和1/MPa或1/kPa或1/Pa。

1.5、对于无量纲参数ρ_i，软化系数演化方程表示为：

式中，ρ_i,0为第三主应力σ₃为零值的ρ_i值，ρ_i,c为σ_i,3等于σ_i,3 ^c时的ρ_i值，为常系数。

1.6、计算滑体主应变在滑面单元上产生的应变(ε_ij')表达式为：

主应变ε₁的方向与x轴间所夹角度(滑面倾角)则可按下式求得：

当γ_xy为正直，且ε_xx>ε_yy时，角为负值，故在上式中用负号。

计算滑面上的单元的应力，计算公式如下：

2、在峰值主应力满足摩尔库伦准则条件下，其破坏面上剪应力与最小主应力的夹角为45°+φ/2(φ为材料摩擦角)，单元剪应力摩尔库伦准则为：

τ_peak＝C+σ_n tanφ (2.1)

其中：式中：τ_peak,σ_n,C分别为破坏面最大剪应力、法向应力和凝聚力，单位为MPa或kPa或Pa。

3、滑面上的剪应力和剪应变满足中国专利一种边坡渐进破坏潜在滑动面的计算方法(申请号：201510658880.6)中的本构方程。

3.1、剪应力-剪应变为四参数本构方程：

τ＝Gγ[1+γ^q/p]^ξ (3.1)

式中：τ、γ分别为剪应力和剪应变，G为剪切模量，p、q、ξ为在不同法向应力下的常系数，τ、G的单位为MPa或kPa或Pa，p、q、ξ为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ξ≤0和1+qξ≠0；临界应变空间满足如下关系式：

p+(1+qξ)γ^q _peak＝0 (3.2)

式中：γ_peak为临界应力对应的应变；临界应变以坐标转换方法加以计算。

3.2、滑体主应变在滑面上产生的法向应变(ε_n')和剪应变(γ')表达式为：

γ'＝(ε₁-ε₃)sin2α (3.4)

式中：α：滑面与最小主应力方向的夹角。

3.3、初始剪切模量方程：

式中：ν：泊松比。

3.4、对于无量纲参数ξ，软化系数演化方程表示为：

ξ＝ξ₀/(1+(ξ₀/ξ_c-1)(σ_n/σ_n ^c)^ψ) (3.6)

式中，ξ₀为法向应力σ_n为零值的ξ值，ξ_c为σ_n等于σ_n ^c时的ξ值，ψ为常系数。

对于同一种地质材料，无量刚软化系数不必具体求解方程(3.6)，而直接利用单元在滑面上求解的驱动剪应力、剪应变和法向应力求解软化系数：求解方法为：

利用破坏面上所得到的法向应力，可以对应地求解出主应力(σ₁,σ₃)，方程如下：

滑面上的法向应力为：

临界剪应力方程为：

求解得临界主应力和围压：

利用方程(3.9，3.10，1.1)，可以获得滑面上临界剪应变、临界剪应力和对应的法向应力。利用这些结果计算常系数：

p＝-(γ_peak)^q(1+qξ) (3.12)

再利用滑面上的实际计算剪应力，以方程(3.1)计算软化系数ξ。

上述计算方法适宜于同种地质材料的渐进破坏过程分析，如果滑面材料与滑体不一致，则滑面计算方法如下：

4、滑面上的剪应力和剪应变满足中国专利一种边坡渐进破坏潜在滑动面的计算方法(申请号：201510658880.6)中的本构方程。

对于不同于滑体材料滑面的剪应力-剪应变四参数本构方程、临界应力方程、临界剪应变方程和软化系数演化方程均与方程(3.1)、(2.1)、(1.4.1或1.4.2)和(3.6)一致。其初始剪切模量方程为：

G＝G₀+b₁σ_n+b₂σ_n ² (4.1)

式中：G₀为法向应力σ_n为零值的G值，b₁,b₂为常系数，单位为无量纲和1/MPa或1/kPa或1/Pa。

5、沿滑面的变形分布特征

对于推移式滑坡沿滑面的位移分布见图1，当M条块(或单元)为临界状态时，沿滑面的力分布特征为：在M+1至N条块间，下滑力等于摩阻力，压应力等于反压应力，在这个阶段，沿滑面的剪应变(γ)以下滑应力等于抗滑剪应力加以计算；在M条块(或单元)下滑应力等于抗滑临界摩阻应力，以此计算该条块的临界剪应变(γ)；在第1至M-1条块(或单元)间，下滑应力大于摩阻应力，其任一条块(或单元)L存在两个位移：第一种位移为第L条块(或单元)本身产生的位移(即：由第L-1的右边和第L+1的左边位移之差加以计算，第二种位移为第L条块(或单元)本身产生的位移加上第L条块(或单元)的底边中心相对该点沿滑床原始位置滑动的位移。第一种位移用于计算滑体的驱动下滑剪应力，第二种位移用于计算滑床提供给滑体的摩阻应力。

L条块(或单元)的第一种底边位移(S_L,j ¹)计算：

S_L,j ¹＝ε_j ^Ll^L，j∈(x,y,z) (5.1)

L条块(或单元)的第二种底边位移(S_L,j ²)计算：

式中：ε_j ⁱ,lⁱ分别为沿滑面第i条块(或单元)在x,y,z轴方向的应变和底边长度。在滑面方向的位移为三者的矢量和，滑面剪应变(γ)(或视剪应变)为上述第一、二种位移除以滑面长度。

6、硬化特征描述

当滑体仅表现出硬化特征时，其计算也可以采用简化的模型，按照中国专利一种边坡渐进破坏潜在滑动面的计算方法(申请号：201510658880.6)，将本构方程(1.1)取ξ＝-1，q＝1，则a,＝1/(Ga”)，b'＝1/(Gp)，在主应力条件下其方程形式与邓肯—张模型一致，此时可以描述材料的硬化行为特征；

式中：a_i’、b_i’为常系数。

在单元主应力满足方程(1.3)条件下，相对应的单元应变定义为峰值应变，满足方程(1.4.1或1.4.2)，在峰值应力和应变下方程(6.1)变为：

定义单元割线模量

则

对方程(6.1)求导，相应的导数为切线模量，在任一应力状态条件下，切向模量E_i(亦即：切线弹性模量)表示为：

从上述方程知，初始弹性模量为1/a_i'，用方程(6.4)，在主应力满足方程(1.3)时的切线模量定义为E_t,i则有：

E_t,i＝a_i’K_scant,i ² (6.6)

在峰值应力时，研究试验曲线的切向模量为E_t,i，假设其具有如下特征：

式中：α_i,k_i为常系数，p_a气压。

方程(6.7)的常系数可以利用三个不同法向应力试验加以确定。

7、针对条分法，边坡渐进破坏计算方法包括如下步骤：

7.1、利用现行方法进行条块划分。

7.2、竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力假设满足一定的应力分布条件；计算相对应的应力，同时计算主应力和旋转角，根据主应力计算对应主应变，通过坐标旋转计算滑动方向和其它方向的应变和位移。

7.3、根据现场实测(或假设)确定第一条块水平和竖直方向的位移，该位移加上上述7.2计算的位移，所得位移除以滑面长度即为滑面剪应变(或视剪应变)，利用该剪应变计算条块底边摩阻应力。

7.4、重复上述7.2和7.3，直至现场判断的临界状态条块为止，当然整个计算可以根据现场实测资料加以对比，从而修正模型参数。

8、数值计算

对于数值计算，按照上述理论，可以按照中国专利一种边坡稳定性计算的滑面边界法(申请号：2014100250810)实施边坡渐进破坏计算。

Claims (1)

1.一种边坡渐进破坏全过程计算方法，其特征在于，包括以下步骤：

当滑体和滑面为同种材料时，

步骤1：对材料实施主应力-应变全过程曲线试验，试验获得峰值主应力、应变和全过程曲线特征，计算滑体材料的物理力学特征方程，包括主应力-主应变方程、峰值应力方程、临界主应变方程、弹性模量变化方程、软化系数方程；

步骤1.1：主应力-主应变方程

主应力-主应变为四参数本构方程：

${\mathrm{\&sigma;}}_i-{\mathrm{\&sigma;}}_3=E_{ii}{\mathrm{\&epsiv;}}_i{(1+{{\mathrm{\&epsiv;}}_i}^{m_i}/S_i)}^{{\mathrm{\&rho;}}_i}$

式中，σ_i、ε_i分别为主应力和主应变，E_ii为初始弹性模量，m_i、S_i、ρ_i为主应力方向的常系数，σ_i、E_ii的单位为MPa或kPa或Pa，m_i、S_i、ρ_i为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ρ_i≤0和1+m_iS_i≠0；临界应变空间满足关系式：

$S_i+(1+m_i{\mathrm{\&rho;}}_i){\left({{\mathrm{\&epsiv;}}_i}^{peak}\right)}^{m_i}=0$

式中，ε_i ^peak为临界主应力对应的主应变；

步骤1.2：假设临界主应力满足摩尔库伦准则，即

$\frac{1}{2}({\mathrm{\&sigma;}}_1-{\mathrm{\&sigma;}}_3)=Ccos\mathrm{\&phi;}+\frac{1}{2}({\mathrm{\&sigma;}}_1+{\mathrm{\&sigma;}}_3)sin\mathrm{\&phi;}$

式中，C为凝聚力，C和σ_i的单位为MPa或kPa或Pa，φ为材料摩擦角；

步骤1.3：设临界主应变空间仅相关于第三主应力，临界主应变ε_i ^peak采用关系式：

${({{\mathrm{\&epsiv;}}_i}^{peak}/a_{i,3})}^2+{\left(\right({\mathrm{\&sigma;}}_3-a_{i,2})/a_{i,1})}^{{\mathrm{\&zeta;}}_{i,N}}=1$

或

式中，a_i,1,a_i,2,a_i,3,ζ_i,N,为常系数；a_i,1,a_i,2单位为MPa或kPa或Pa，a_i,3,ζ_i,N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²；

步骤1.4：初始弹性模量方程：

E_ii＝E_ii,0+b_i,1σ₃+b_i,2σ₃ ²

式中，E_ii,0为第三主应力σ₃为零值的E_ii值，b_i,1,b_i,2为常系数，单位分别为无量纲和1/MPa或1/kPa或1/Pa；

步骤1.5：对无量纲参数ρ_i，软化系数演化方程表示为：

式中，ρ_i,0为第三主应力σ₃为零值的ρ_i值，ρ_i,c为σ_i,3等于σ_i,3 ^c时的ρ_i值，为常系数；

步骤1.6：计算滑体主应变在滑面单元上产生的应变ε_ij'表达式为：

${\mathrm{\&epsiv;}}_1=\frac{1}{2}\&lsqb;({{\mathrm{\&epsiv;}}_{xx}}^{\&prime;}+{{\mathrm{\&epsiv;}}_{yy}}^{\&prime;})+\sqrt{{({{\mathrm{\&epsiv;}}_{xx}}^{\&prime;}-{{\mathrm{\&epsiv;}}_{yy}}^{\&prime;})}^2+{\mathrm{\&gamma;}}^{\&prime;}{{}_{xy}}^2}\&rsqb;$

${\mathrm{\&epsiv;}}_2=\frac{1}{2}\&lsqb;({{\mathrm{\&epsiv;}}_{xx}}^{\&prime;}+{{\mathrm{\&epsiv;}}_{yy}}^{\&prime;})-\sqrt{{({{\mathrm{\&epsiv;}}_{xx}}^{\&prime;}-{{\mathrm{\&epsiv;}}_{yy}}^{\&prime;})}^2+{\mathrm{\&gamma;}}^{\&prime;}{{}_{xy}}^2}\&rsqb;$

主应变ε₁的方向与x轴间所夹角度按式子和求得，当γ_xy为正直，且ε_xx>ε_yy时，角为负值；计算滑面上的单元的应力，计算公式为：

步骤2：在峰值主应力满足摩尔库伦准则条件下，其破坏面上剪应力与最小主应力的夹角为45°+φ/2，φ为材料摩擦角，单元剪应力摩尔库伦准则为：

τ_peak＝C+σ_n tanφ

其中，τ_peak,σ_n,C分别为破坏面最大剪应力、法向应力和凝聚力，单位为MPa或kPa或Pa；

步骤3：滑面上的剪应力和剪应变满足的本构方程，具体如下：

步骤3.1：剪应力-剪应变为四参数本构方程：

τ＝Gγ[1+γ^q/p]^ξ

式中，τ、γ分别为剪应力和剪应变，G为剪切模量，p、q、ξ为在不同法向应力下的常系数，τ、G的单位为MPa或kPa或Pa，p、q、ξ为无单位参数；

其行为描述为：对于具有软化特征的材料行为，则有-1<ξ≤0和1+qξ≠0；临界应变空间满足关系式：

p+(1+qξ)γ^q _peak＝0

式中，γ_peak为临界应力对应的应变；

步骤3.2：滑体主应变在滑面上产生的法向应变ε_n'和剪应变γ'表达式：

${{\mathrm{\&epsiv;}}_n}^{\&prime;}=\frac{1}{2}({\mathrm{\&epsiv;}}_1+{\mathrm{\&epsiv;}}_3)+\frac{1}{2}({\mathrm{\&epsiv;}}_1-{\mathrm{\&epsiv;}}_3)cos2\mathrm{\&alpha;}$

γ'＝(ε₁-ε₃)sin2α

式中，α为滑面与最小主应力方向的夹角；

步骤3.3：初始剪切模量方程：

$G=\frac{E_{11}}{2(1+\mathrm{\&nu;})}=\frac{E_{11}}{2(1+E_{11}/E_{33})}$

式中，ν为泊松比；

步骤3.4：对于无量纲参数ξ，软化系数演化方程表示为：

ξ＝ξ₀/(1+(ξ₀/ξ_c-1)(σ_n/σ_n ^c)^ψ)

ξ₀为法向应力σ_n为零值的ξ值，ξ_c为σ_n等于σ_n ^c时的ξ值，ψ为常系数；

对于同一种地质材料，直接利用单元在滑面上求解的驱动剪应力、剪应变和法向应力求解软化系数：求解方法为：

利用破坏面上所得到的法向应力，对应地求解出主应力σ₁,σ₃，方程为：

滑面上的法向应力：

临界剪应力方程：

求解得临界主应力和围压：

利用方程 获得滑面上临界剪应变、临界剪应力和对应的法向应力，计算常系数：

$q=\frac{{({\mathrm{\&tau;}}_{peak}/({\mathrm{\&gamma;}}_{peak}G\left)\right)}^{1/\mathrm{\&xi;}}}{\mathrm{\&xi;}(1-{\left({\mathrm{\&tau;}}_{peak}/\left({\mathrm{\&gamma;}}_{peak}G\right)\right)}^{1/\mathrm{\&xi;}})}$

p＝-(γ_peak)^q(1+qξ)

再利用滑面上的实际计算剪应力，以方程τ＝Gγ[1+γ^q/p]^ξ计算软化系数ξ；

若滑面材料与滑体不一致，则滑面计算方法为：

步骤4：对于不同于滑体材料滑面的剪应力-剪应变四参数本构方程、临界应力方程、临界剪应变方程和软化系数演化方程均与方程：

τ＝Gγ[1+γ^q/p]^ξ

τ_peak＝C+σ_ntanφ

或二者之一

ξ＝ξ₀/(1+(ξ₀/ξ_c-1)(σ_n/σ_n ^c)^ψ)一致，其初始剪切模量方程为：

G＝G₀+b₁σ_n+b₂σ_n ²

式中，G₀为法向应力σ_n为零值的G值，b₁,b₂为常系数，单位为无量纲和1/MPa或1/kPa或1/Pa；

步骤5：沿滑面的变形分布特征

当M条块或单元为临界状态时，沿滑面的力分布特征为：在M+1至N条块间，下滑力等于摩阻力，压应力等于反压应力，在这个阶段，沿滑面的剪应变γ以下滑应力等于抗滑剪应力加以计算；在M条块或单元下滑应力等于抗滑临界摩阻应力，以此计算该条块的临界剪应变γ；

在第1至M-1条块或单元间，下滑应力大于摩阻应力，其任一条块或单元L存在两个位移：第一种位移为第L条块或单元本身产生的位移，即：由第L-1的右边和第L+1的左边位移之差加以计算，第二种位移为第L条块或单元本身产生的位移加上第L条块或单元的底边中心相对该点沿滑床原始位置滑动的位移；第一种位移用于计算滑体的驱动下滑剪应力，第二种位移用于计算滑床提供给滑体的摩阻应力；

L条块或单元的第一种底边位移S_L,j ¹计算：

S_L,j ¹＝ε_j ^Ll^L，j∈(x,y,z)

L条块或单元的第二种底边位移S_L,j ²计算：

${S_{L,j}}^2=\underset{i=M}{\overset{L}{\&Sigma;}}{{\mathrm{\&epsiv;}}_j}^i{l_j}^i,j\&Element;(x,y,z)$

式中，ε_j ⁱ,lⁱ分别为沿滑面第i条块或单元在x、y、z轴方向的应变和底边长度；在滑面方向的位移为三者的矢量和，滑面剪应变γ为上述第一、二种位移除以滑面长度；

步骤6：硬化特征描述

当滑体仅表现出硬化特征时，其计算采用简化的模型，将本构方程取ξ＝-1，q＝1，则a’＝1/(Ga”)，b'＝1/(Gp)，在主应力条件下其方程形式与邓肯—张模型一致，此时描述材料的硬化行为特征：

${\mathrm{\&sigma;}}_i-{\mathrm{\&sigma;}}_3=\frac{{\mathrm{\&epsiv;}}_i}{{a_i}^{,}+{b_i}^{,}{\mathrm{\&epsiv;}}_i}$

式中，a_i’、b_i’为常系数；

在单元主应力满足方程条件下，相对应的单元应变定义为峰值应变，满足方程或在峰值应力和应变下，方程变为：

${a_i}^{,}+{b_i}^{,}{\mathrm{\&epsiv;}}_{i,peak}=\frac{1}{{({\mathrm{\&sigma;}}_i-{\mathrm{\&sigma;}}_3)}_{peak}/{\mathrm{\&epsiv;}}_{i,peak}}$

定义单元割线模量则

对方程求导，相应的导数为切线模量，在任一应力状态条件下，切向模量E_i表示为：

$E_i=\frac{{a_i}^{,}}{{({a_i}^{,}+{b_i}^{,}{\mathrm{\&epsiv;}}_i)}^2}$

得到初始弹性模量为1/a_i'，利用方程在主应力满足方程时的切线模量定义为E_t,i，则有：

E_t,i＝a_i’K_scant,i ²

在峰值应力时，研究试验曲线的切向模量为E_t,i，设其具有如下特征：

σ₃≠0

式中，α_i,k_i为常系数，p_a为大气压；

方程的常系数利用三个不同法向应力试验加以确定；

步骤7：针对条分法，边坡渐进破坏计算方法包括如下步骤：

步骤7.1：利用现行方法进行条块划分；

步骤7.2：竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力满足应力分布条件；计算相对应的应力，同时计算主应力和旋转角，根据主应力计算对应主应变，通过坐标旋转计算滑动方向和其它方向的应变和位移；

步骤7.3：根据现场实测或假设确定第一条块水平和竖直方向的位移，该位移加上步骤7.2计算的位移，所得位移除以滑面长度即为滑面剪应变或视剪应变，利用该剪应变计算条块底边摩阻应力；

步骤7.4：重复步骤7.2和步骤7.3，直至现场判断的临界状态条块为止；

步骤8：实施边坡渐进破坏的数值计算。