CN105335607A - 一种边坡渐进破坏潜在滑动面的计算方法 - Google Patents

Abstract

本发明提出了一种边坡渐进破坏潜在滑动面的计算方法——简称破坏角转动法。针对现有技术的缺陷，在假设地质材料破坏时满足临界应力状态所对应的最大剪切应力面与最小主应力轴夹角为（）条件，基于边坡在不同外加荷载和重力荷载作用时主应力方向在不同位置发生转动的事实，实施边坡潜在滑动面的搜索计算，从而决定潜在滑动面。并定义了破坏率和破坏比概念。破坏角转动法保证了边坡破坏过程中，破坏点的应力状态处于临界应力状态，且破坏过程中，破坏路径随应力的变化而变化，结合破坏率和破坏比概念，破坏路径的本构关系考虑不同法向应力作用下的软化特征，可以在数值计算的基础上实施边坡潜在滑动面的求解。

Description

一种边坡渐进破坏潜在滑动面的计算方法

技术领域

本发明涉及土木工程和地质灾害、基坑等的防治、评价和预测预报等技术领域，特别涉及地质灾害、基坑等的稳定性分析、评价、预测预报、防治措施的制定，本发明实现了地质灾害和基坑等的渐进破坏过程的的潜在滑动面决定和稳定性评价，对边坡和基坑等的防治及预测预报等具有巨大的推动作用。

背景技术

边坡的稳定性评价均是建立在极限平衡状态的假设上，其广泛采用的方法为：瑞典法、简化Bishop法、Janbu法、传递系数法、Sarma法、楔形体法、Fellenius法、有限元强度折减法等十几种边坡稳定性计算方法。其潜在滑动面的决定也是建立在临界应力状态理论之上，然而现场边坡破坏是渐进的，滑面部分处于临界应力状态，部分处于破坏后区或峰值前应力状态区，现行的极限平衡状态法获取的潜在滑动面是难以与现场实际相符的，鉴于此，本发明提出了一种边坡渐进破坏潜在滑动面的计算方法—以下简称破坏角转动法，该方法使潜在滑动面的决定向实际现场推动了一大步。

发明内容

本发明的目的在于提出一种边坡渐进破坏潜在滑动面的计算方法，基于边坡破坏是渐进的，其主应力轴是发生转动的，但最大剪切面上的破坏角相对最小主应力是恒定的，在此基础上获得破坏角的转动规律，实施边坡潜在滑动面的搜索计算，从而决定潜在滑动面(如图1)。并定义了破坏率(滑体作用于滑面的下滑剪应力(或拉应力)除以滑床作用于滑面临界摩阻应力(或临界拉应力)的绝对值，当大于100％时取100％)和破坏比(沿可能滑动的滑面面积与破坏率乘积之和除以总面积)概念。破坏角转动法保证了边坡破坏过程中，破坏点的应力状态处于临界应力状态，且破坏过程中，破坏路径随应力的变化而变化，结合破坏率和破坏比概念，破坏路径的本构关系考虑不同法向应力作用下的软化特征，可以在数值计算的基础上实施边坡潜在滑动面的求解。

本发明一种边坡渐进破坏潜在滑动面的计算方法，包括如下步骤：

(1)对滑体材料实施剪应力--剪应变全过程曲线试验，试验获得峰值应力、应变和全过程曲线；

(2)通过峰值应力决定凝聚力C、滑面摩擦角值，以峰值应变决定常系数a₁,a₂,a₃大小，变化曲线特征决定剪切模量G、临界法向应力σ_n ^crit、常系数ξ、α、k_n；

(3)按现行方法建立数值计算模型，在考虑剪破坏分布的同时，同时考虑拉破坏分布区；

(4)在考虑应变软化本构模型数值计算的基础上，计算现状边坡每点的破坏率、破坏面和破坏比，以组合的方式提出不同的可能破坏路径；

(5)对于单元，利用单元破坏剪应力面与最小主应力的夹角为：计算最大主应力相对竖直方向的转动角δ，从而决定滑面相对于水平面的转角：所述转动角δ，二维计算公式为tan2δ＝-2τ_xy/(σ_xx-σ_yy)，三维计算公式为tan2δ_xx＝-2τ_xy/(σ_xx-σ_yy)，tan2δ_yy＝-2τ_zy/(σ_yy-σ_zz)，tan2δ_zz＝-2τ_zx/(σ_zz-σ_xx)；

(6)对于可能施加的荷载或位移工况，分步施加相应的工况，对在不同工况下可能破坏模式进行搜寻，将潜在滑动面转动角连续化，计算对应的边坡稳定系数，从而决定潜在滑动面；

(7)对于具有软化和硬化特征的滑面剪应力--剪应变满足如下本构方程，本发明的滑面本构方程如下：

(7.1)剪应力-剪应变方程

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

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

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

(7.2)软化特征

对于具有软化特征的材料行为，则有：-1＜ξ≤0和1+qξ≠0。临界应变空间(临界应变定义为：临界应力对应的应变)满足如下关系式：

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

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

假设临界应力空间τ_peak满足摩尔库伦准则(注：也可以满足其它相关准则)：

式中：C为凝聚力，σ_n为法向应力，C和σ_n的单位为MPa或kPa或Pa，为滑面摩擦角；

临界应变空间可以假设仅相关于法向应力，临界应变γ_peak采用如下关系式：

${({\mathrm{\γ}}_{peak}/a_3)}^2+{\left(\right({\mathrm{\σ}}_n-a_2)/a_1)}^{{\mathrm{\ζ}}_N}=1---\left(\mathrm{7.4.1}\right)$

或 ${{\mathrm{\γ}}_{peak}}^2=a_1^0+a_2^0{\mathrm{\σ}}_n+a_3^0{{\mathrm{\σ}}_n}^2---\left(\mathrm{7.4.2}\right)$

式中：a₁,a₂,a₃,ζ_N,为常系数；a₁,a₂单位为MPa或kPa或Pa，a₃,ζ_N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²；

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

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

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

式中，ξ₀为法向应力σ_n为零值的ξ值，ξ_c为σ_n等于σ_n ^c时的ξ值，为常系数；该关系式可由不同的法向压力试验曲线而获得。

(7.3)硬化特征

当地质材料的法向应力大于临界法向应力时，则没有明显的峰值应力，本发明提出两种计算方法:

(7.3.1)方法一

将本构方程(7.1)取ξ＝-1，q＝1，则a’＝1/(Ga”)，b′＝1/(Gp)，其方程形式与邓肯—张模型一致，此时只能描述材料的弹塑性硬化行为特征。

$\mathrm{\τ}=\frac{\mathrm{\γ}}{a^{,}+b^{,}\mathrm{\γ}}---\left(7.7\right)$

式中：a’、b’、a”为常系数。

在峰值应力条件下，方程(7.7)变为：

$a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{{\mathrm{\τ}}_{peak}/{\mathrm{\γ}}_{peak}}---\left(7.8\right)$

定义割线模量 $k_{scant}=\frac{{\mathrm{\τ}}_{peak}}{{\mathrm{\γ}}_{peak}}---\left(7.9\right)$

则 $a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{K_{cant}}---\left(7.10\right)$

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

$G_i=\frac{a^{,}}{{(a^{,}+b^{,}\mathrm{\γ})}^2}---\left(7.11\right)$

利用方程(7.11)，在最大应力时的切线模量G_t则有：

G_i＝a’K_cant ²(7.12)

众所周知，对于没有明显峰值应力的塑性硬化行为，目前试验难以获得峰值应力，峰值应力的选取则必须满足目前的各种应力准则(如：摩尔库伦准则)，对应的剪应变也满足本发明提出的应变空间方程；在峰值应力时，研究试验曲线的切向模量，记为G_t，假设其具有如下特征：

$G_t=\mathrm{\α}({\mathrm{\σ}}_n-{\mathrm{\σ}}_n^{crit}){({\mathrm{\σ}}_n/{\mathrm{\σ}}_n^{crit})}^{k_n}---\left(7.13\right)$

${\mathrm{\σ}}_n^{crit}\≤{\mathrm{\σ}}_n\≤{\mathrm{\σ}}_n^{\max },$ α、k_n为常系数。

方程(7.13)的特点为：

当切线模量等于0，此时曲线呈现近似理想弹塑性模型特征，当σ_n达到一定值时，曲线呈现线性特征，此时的法向应力理论上试验可以决定，记为相对应的切向模量应该等于则具有如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^{max}-{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^{max}/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^{max}---\left(7.14\right)$

在法向应力的范围内，取某一法向应力试验确定对应的切向模量：G^a，可以获得如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^a-{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^a/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^a---\left(7.15\right)$

则由方程(7.14、7.15)，可以决定常系数：

$k_n=\frac{ln\left(G^{max}\right({{\mathrm{\σ}}_n}^{\mathrm{\α}}-{{\mathrm{\σ}}_n}^{crit})/(G^a\left({{\mathrm{\σ}}_n}^{\max }-{{\mathrm{\σ}}_n}^{crit}\right))}{ln({{\mathrm{\σ}}_n}^{\max }/{{\mathrm{\σ}}_n}^{\mathrm{\α}})}$

和 $\mathrm{\α}=G^{max}/\left(\right({{\mathrm{\σ}}_n}^{\max }-{\mathrm{\σ}}_n^{crit}\left){\left({{\mathrm{\σ}}_n}^{max}/{\mathrm{\σ}}_n^{crit}\right)}^{k_n}\right)---\left(7.16\right)$

在确定某一法向应力σ_n条件下的峰值应力切向模量G_t之后，按方程(7.12)可以决定a’，按方程(7.10)可以决定b’，至此推广新的邓肯--张模型各参数得以确定。

(7.3.2)方法二

取方程(7.1)的ξ＝-1，则方程变为：

$\mathrm{\τ}=\frac{G\mathrm{\γ}}{1+{\mathrm{\γ}}^q/p}---\left(7.17\right)$

在峰值应力下，则有：

${\mathrm{\τ}}^{peak}/{\mathrm{\γ}}_{peak}=\frac{G}{1+{\left({\mathrm{\γ}}_{peak}\right)}^q/p}---\left(7.18\right)$

${{\mathrm{\γ}}_{peak}}^q/p=\frac{G}{K_{scant}}-1---\left(7.19\right)$

同理对方程(7.17)求导，相应的导数为切线模量：

$\frac{\∂\mathrm{\τ}}{\∂\mathrm{\γ}}=\frac{G}{(1+{\mathrm{\γ}}^q/p)}-\frac{{\mathrm{Gq\γ}}^q/p}{{(1+{\mathrm{\γ}}^q/p)}^2}---\left(7.20\right)$

当峰值应力满足目前摩尔库伦准则时，则峰值应变也满足方程(7.4)，在峰值应力时，其切线模量记为：G_t。

利用方程(7.18、7.19)，在峰值应力时，切线模量则有方程：

$G_t=K_{scant}\[1+\frac{{\mathrm{qK}}_{scant}}{G}(1-\frac{G}{K_{scant}})\]---\left(7.21\right)$

对于峰值应力对应的切线模量，按照方程(7.13)求解。利用方程(7.21)求解参数q，利用方程(7.19)求解p。

(8)针对广泛使用的条分法，破坏角转动法潜在滑动面决定如下：

(8.1)利用现行方法进行条块划分；

(8.2)竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力假设满足一定的应力分布条件(如：直线分布或抛物线分布)；

(8.3)条块底边摩阻应力计算，按照步骤(7)进行计算；

(9)潜在滑面决定条件，分两种情况：

(9.1)数值计算

数值计算决定潜在滑动面，按传统的强度折减和可能的荷载(或位移)工况分步施加两种方法进行。

(9.1.1)传统强度折减

在本发明提出的破坏角转动法的基础上，对临界抗剪强度进行折减，直至位于临空面的破坏单元处于临界应力状态。

(9.1.2)荷载(或位移)工况施加法

在本发明提出的破坏角转动法基础上，对可能破坏发生施加相应的荷载(或位移)工况，直至位于临空面的破坏单元处于临界应力状态。

(9.2)条分法

条分法决定潜在滑动面，按传统的强度折减和施加荷载(或位移)工况两种方法进行。

(9.2.1)传统强度折减

在本发明提出的破坏角转动法的基础上，对条块底边临界抗剪强度进行折减，直至位于临空面的破坏条块处于极限平衡状态。

(9.2.2)荷载(或位移)工况施加法

在本发明提出的破坏角转动法基础上，对可能破坏施加相应的荷载(或位移)工况，直至位于临空面的破坏条块处于极限平衡状态。

上述两种方法之中，强度折减法没有物理意义，计算所得应力和位移逻辑上不能与实际现场相比较。

本发明边坡渐进破坏潜在滑动面的计算方法的优点在于：

现行边坡渐进破坏潜在滑动面的决定方法主要是利用极限平衡状态下的力学参数(如：凝聚力C,摩擦角)，采用现行的极限状态稳定性(如:瑞典圆弧法等)系数搜索法决定潜在滑动面。这种决定潜在滑动面的方法缺点为:1)整个滑动面处于临界应力状态，然而边坡滑动面破坏是渐进的，2)另外边坡破坏时，破坏点处于临界应力状态，其余是处于破坏后区或峰值前应力状态，这种破坏方式当今的方法是难以描述的等。

针对这些缺点，在假设地质材料破坏时满足临界应力状态所对应的最大剪切应力面与最小主应力轴夹角为条件，基于边坡在不同外加荷载和重力荷载作用时主应力方向在不同位置发生转动(转动角δ)的事实，实施边坡潜在滑动面的搜索计算，从而决定潜在滑动面。并定义了破坏率和破坏比的概念，提出了荷载或位移工况施加法。破坏角转动法保证了边坡破坏过程中，破坏点的应力状态处于临界应力状态，且破坏过程中，破坏路径随应力的变化而变化，结合破坏率和破坏比概念，破坏路径的本构关系考虑不同法向应力作用下的软化特征，可以在数值计算的基础上实施边坡潜在滑动面的求解。

附图说明

图1为边坡渐进破坏潜在滑动面破坏角转动决定法示意图。

图中，σ_xx、σ_yy、τ_xy、σ₁₁、σ₂₂、δ分别为X轴方向应力、Y轴方向应力、剪应力、摩擦角、最大主应力、最小主应力、转动角。

具体实施方式

下面结合附图，对本发明进行进一步说明。如图1所示，本发明一种边坡渐进破坏潜在滑动面的计算方法，包括如下步骤：

(1)对滑体材料实施剪应力--剪应变全过程曲线试验，试验获得峰值应力、应变和全过程曲线；

(2)通过峰值应力决定凝聚力C、滑面摩擦角值，以峰值应变决定常系数a₁,a₂,a₃大小，变化曲线特征决定剪切模量G、临界法向应力σ_n ^crit、常系数ξ、α、k_n等；

(3)按现行方法建立数值计算模型，在考虑剪破坏分布的同时，同时也可以考虑拉破坏分布区；

(4)在考虑应变软化本构模型数值计算的基础上，计算现状边坡每点的破坏率、破坏面和破坏比，以组合的方式提出不同的可能破坏路径；

(5)对于单元，利用单元破坏剪应力面与最小主应力的夹角为：计算最大主应力相对竖直方向的转动角δ(如：二维问题：tan2δ＝-2τ_xy/(σ_xx-σ_yy)，三维问题：tan2δ_xx＝-2τ_xy/(σ_xx-σ_yy)，tan2δ_yy＝-2τ_zy/(σ_yy-σ_zz)，tan2δ_zz＝-2τ_zx/(σ_zz-σ_xx)，从而决定滑面相对于水平面的转角：

(6)对于可能施加的荷载或位移工况，分步施加相应的工况，对在不同工况下可能破坏模式进行搜寻，将潜在滑动面转动角连续化，计算对应的边坡稳定系数，从而决定潜在滑动面；

(7)对于具有软化和硬化特征的滑面剪应力--剪应变满足如下本构方程，本发明的滑面本构方程如下：

(7.1)剪应力-剪应变方程

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

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

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

(7.2)软化特征

对于具有软化特征的材料行为，则有：-1＜ξ≤0和1+qξ≠0。临界应变空间(临界应变定义为：临界应力对应的应变)满足如下关系式：

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

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

假设临界应力空间τ_peak满足摩尔库伦准则(注：也可以满足其它相关准则)：

式中：C为凝聚力，σ_n为法向应力，C和σ_n的单位为MPa或kPa或Pa，为滑面摩擦角；

临界应变空间可以假设仅相关于法向应力，临界应变γ_peak采用如下关系式：

${({\mathrm{\γ}}_{peak}/a_3)}^2+{\left(\right({\mathrm{\σ}}_n-a_2)/a_1)}^{{\mathrm{\ζ}}_N}=1---\left(\mathrm{7.4.1}\right)$

或 ${{\mathrm{\γ}}_{peak}}^2=a_1^0+a_2^0{\mathrm{\σ}}_n+a_3^0{{\mathrm{\σ}}_n}^2---\left(\mathrm{7.4.2}\right)$

式中：a₁,a₂,a₃,ζ_N,为常系数；a₁,a₂单位为MPa或kPa或Pa，a₃,ζ_N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²；

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

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

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

式中，ξ₀为法向应力σ_n为零值的ξ值，ξ_c为σ_n等于σ_n ^c时的ξ值，为常系数；该关系式可由不同的法向压力试验曲线而获得。

(7.3)硬化特征

当地质材料的法向应力大于临界法向应力时，则没有明显的峰值应力，本发明提出两种计算方法:

(7.3.1)方法一

将本构方程(7.1)取ξ＝-1，q＝1，则a’＝1/(Ga”)，b′＝1/(Gp)，其方程形式与邓肯—张模型一致，此时只能描述材料的弹塑性硬化行为特征。

$\mathrm{\τ}=\frac{\mathrm{\γ}}{a^{,}+b^{,}\mathrm{\γ}}---\left(7.7\right)$

式中：a’、b’、a”为常系数。

在峰值应力条件下，方程(7.7)变为：

$a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{{\mathrm{\τ}}_{peak}/{\mathrm{\γ}}_{peak}}---\left(7.8\right)$

定义割线模量 $k_{scant}=\frac{{\mathrm{\τ}}_{peak}}{{\mathrm{\γ}}_{peak}}---\left(7.9\right)$

则 $a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{K_{cant}}---\left(7.10\right)$

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

$G_i=\frac{a^{,}}{{(a^{,}+b^{,}\mathrm{\γ})}^2}---\left(7.11\right)$

利用方程(7.11)，在最大应力时的切线模量G_t则有：

G_i＝a’K_cant ²(7.12)

众所周知，对于没有明显峰值应力的塑性硬化行为，目前试验难以获得峰值应力，峰值应力的选取则必须满足目前的各种应力准则(如：摩尔库伦准则)，对应的剪应变也满足本发明提出的应变空间方程；在峰值应力时，研究试验曲线的切向模量，记为G_t，假设其具有如下特征：

$G_t=\mathrm{\α}({\mathrm{\σ}}_n-{\mathrm{\σ}}_n^{crit}){({\mathrm{\σ}}_n/{\mathrm{\σ}}_n^{crit})}^{k_n}---\left(7.13\right)$

${\mathrm{\σ}}_n^{crit}\≤{\mathrm{\σ}}_n\≤{\mathrm{\σ}}_n^{\max },$ α、k_n为常系数。

方程(7.13)的特点为：

当切线模量等于0，此时曲线呈现近似理想弹塑性模型特征，当σ_n达到一定值时，曲线呈现线性特征，此时的法向应力理论上试验可以决定，记为相对应的切向模量应该等于则具有如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^{\max }-{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^{\max }/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^{max}---\left(7.14\right)$

在法向应力的范围内，取某一法向应力试验确定对应的切向模量：G^a，可以获得如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^a-{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^a/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^a---\left(7.15\right)$

则由方程(7.14、7.15)，可以决定常系数：

$k_n=\frac{ln\left(G^{max}\right({{\mathrm{\σ}}_n}^{\mathrm{\α}}-{{\mathrm{\σ}}_n}^{crit})/(G^a\left({{\mathrm{\σ}}_n}^{\max }-{{\mathrm{\σ}}_n}^{crit}\right))}{ln({{\mathrm{\σ}}_n}^{max}/{{\mathrm{\σ}}_n}^{\mathrm{\α}})}$

和 $\mathrm{\α}=G^{max}/\left(\right({{\mathrm{\σ}}_n}^{\max }-{\mathrm{\σ}}_n^{crit}\left){\left({{\mathrm{\σ}}_n}^{\max }/{\mathrm{\σ}}_n^{crit}\right)}^{k_n}\right)---\left(7.16\right)$

在确定某一法向应力σ_n条件下的峰值应力切向模量G_t之后，按方程(7.12)可以决定a’，按方程(7.10)可以决定b’，至此推广新的邓肯--张模型各参数得以确定。

(7.3.2)方法二

取方程(7.1)的ξ＝-1，则方程变为：

$\mathrm{\τ}=\frac{G\mathrm{\γ}}{1+{\mathrm{\γ}}^q/p}---\left(7.17\right)$

在峰值应力下，则有：

${\mathrm{\τ}}^{peak}/{\mathrm{\γ}}_{peak}=\frac{G}{1+{\left({\mathrm{\γ}}_{peak}\right)}^q/p}---\left(7.18\right)$

${{\mathrm{\γ}}_{peak}}^q/p=\frac{G}{K_{scant}}-1---\left(7.19\right)$

同理对方程(7.17)求导，相应的导数为切线模量：

$\frac{\∂\mathrm{\τ}}{\∂\mathrm{\γ}}=\frac{G}{(1+{\mathrm{\γ}}^q/p)}-\frac{{\mathrm{Gq\γ}}^q/p}{{(1+{\mathrm{\γ}}^q/p)}^2}---\left(7.20\right)$

当峰值应力满足目前摩尔库伦准则时，则峰值应变也满足方程(7.4)，在峰值应力时，其切线模量记为：G_t。

利用方程(7.18、7.19)，在峰值应力时，切线模量则有方程：

$G_t=K_{scant}\[1+\frac{{\mathrm{qK}}_{scant}}{G}(1-\frac{G}{K_{scant}})\]---\left(7.21\right)$

对于峰值应力对应的切线模量，按照方程(7.13)求解。利用方程(7.21)求解参数q，利用方程(7.19)求解p。

(8)针对广泛使用的条分法，破坏角转动法潜在滑动面决定如下：

(8.1)利用现行方法进行条块划分；

(8.2)竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力假设满足一定的应力分布条件(如：直线分布或抛物线分布)；

(8.3)条块底边摩阻应力计算，按照步骤(7)进行计算；

(9)潜在滑面决定条件，分两种情况：

(9.1)数值计算

数值计算决定潜在滑动面，按传统的强度折减和可能的荷载(或位移)工况分步施加两种方法进行。

(9.1.1)传统强度折减

在本发明提出的破坏角转动法的基础上，对临界抗剪强度进行折减，直至位于临空面的破坏单元处于临界应力状态。

(9.1.2)荷载(或位移)工况施加法

在本发明提出的破坏角转动法基础上，对可能破坏发生施加相应的荷载(或位移)工况，直至位于临空面的破坏单元处于临界应力状态。

(9.2)条分法

条分法决定潜在滑动面，按传统的强度折减和施加荷载(或位移)工况两种方法进行。

(9.2.1)传统强度折减

在本发明提出的破坏角转动法的基础上，对条块底边临界抗剪强度进行折减，直至位于临空面的破坏条块处于极限平衡状态。

(9.2.2)荷载(或位移)工况施加法

在本发明提出的破坏角转动法基础上，对可能破坏施加相应的荷载(或位移)工况，直至位于临空面的破坏条块处于极限平衡状态。

上述两种方法之中，强度折减法没有物理意义，计算所得应力和位移逻辑上不能与实际现场相比较。

Claims (5)

1.一种边坡渐进破坏潜在滑动面的计算方法，其特征在于包括如下步骤：

(1)对滑体材料实施剪应力--剪应变全过程曲线试验，试验获得峰值应力、应变和全过程曲线；

(2)通过峰值应力决定凝聚力C、滑面摩擦角值，通过峰值应变决定常系数a₁,a₂,a₃大小，通过变化曲线特征决定剪切模量G、临界法向应力σ_n ^crit、常系数ξ、α、k_n；

(3)建立数值计算模型，在考虑剪破坏分布的同时，考虑拉破坏分布区；

(4)在考虑应变软化本构模型数值计算的基础上，计算现状边坡每点的破坏率、破坏面和破坏比，以组合的方式提出不同的可能破坏路径；

(5)对于单元，利用单元破坏剪应力面与最小主应力的夹角为计算最大主应力相对竖直方向的转动角δ，从而决定滑面相对于水平面的转角所述转动角δ，二维计算公式为tan2δ＝-2τ_xy/(σ_xx-σ_yy)，三维计算公式为tan2δ_xx＝-2τ_xy/(σ_xx-σ_yy)，tan2δ_yy＝-2τ_zy/(σ_yy-σ_zz)，tan2δ_zz＝-2τ_zx/(σ_zz-σ_xx)；

(6)对于可能施加的荷载或位移工况，分步施加相应的工况，对在不同工况下可能破坏模式进行搜寻，将潜在滑动面转动角连续化，计算对应的边坡稳定系数，从而决定潜在滑动面。

2.如权利要求1所述的计算方法，其特征在于：对于具有软化和硬化特征的滑面剪应力--剪应变满足如下本构方程，

(7.1)剪应力-剪应变方程

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

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

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

(7.2)软化特征

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

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

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

假设临界应力空间τ_peak满足摩尔库伦准则：

式中：C为凝聚力，σ_n为法向应力，C和σ_n的单位为MPa或kPa或Pa，为滑面摩擦角；

假设临界应变空间仅相关于法向应力，临界应变γ_peak采用如下关系式：

${({\mathrm{\γ}}_{peak}/a_3)}^2+{\left(\right({\mathrm{\σ}}_n-a_2)/a_1)}^{{\mathrm{\ζ}}_N}=1---\left(\mathrm{7.4.1}\right)$

或 ${{\mathrm{\γ}}_{peak}}^2=a_1^0+a_2^0{\mathrm{\σ}}_n+a_3^0{{\mathrm{\σ}}_n}^2---\left(\mathrm{7.4.2}\right)$

式中：为常系数；a₁,a₂单位为MPa或kPa或Pa，a₃,ζ_N为无量纲系数，或的量纲为1/MPa,1/MPa²或1/kPa,1/kPa²或1/Pa,1/Pa²；

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

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

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

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

(7.3)硬化特征

当地质材料的法向应力大于临界法向应力时，则没有明显的峰值应力，采用以下两种计算方法：

(7.3.1)方法一

将本构方程(7.1)取ξ＝-1，q＝1，则a’＝1/(Ga”)，b′＝1/(Gp)，其方程形式与邓肯—张模型一致，此时只能描述材料的弹塑性硬化行为特征；

$\mathrm{\τ}=\frac{\mathrm{\γ}}{a^{,}+b^{,}\mathrm{\γ}}---\left(7.7\right)$

式中：a’、b’、a”为常系数；

在峰值应力条件下，方程(7.7)变为：

$a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{{\mathrm{\τ}}_{peak}/{\mathrm{\γ}}_{peak}}---\left(7.8\right)$

定义割线模量 $k_{scant}=\frac{{\mathrm{\τ}}_{peak}}{{\mathrm{\γ}}_{peak}}---\left(7.9\right)$

则 $a^{,}+b^{,}{\mathrm{\γ}}_{peak}=\frac{1}{K_{cant}}---\left(7.10\right)$

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

$G_i=\frac{a^{,}}{{(a^{,}+b^{,}\mathrm{\γ})}^2}---\left(7.11\right)$

利用方程(7.11)，在最大应力时的切线模量G_t则有：

G_i＝a’K_cant ²(7.12)

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

$G_t=\mathrm{\α}({\mathrm{\σ}}_n-{\mathrm{\σ}}_n^{crit}){({\mathrm{\σ}}_n/{\mathrm{\σ}}_n^{crit})}^{k_n}---\left(7.13\right)$

${\mathrm{\σ}}_n^{crit}\≤{\mathrm{\σ}}_n\≤{\mathrm{\σ}}_n^{\max },$ α、k_n为常系数；

方程(7.13)的特点为：

当切线模量等于0，此时曲线呈现近似理想弹塑性模型特征，当σ_n达到一定值时，曲线呈现线性特征，此时的法向应力理论上试验可以决定，记为相对应的切向模量应该等于则具有如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^{max}/{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^{max}/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^{max}---\left(7.14\right)$

在法向应力的范围内，取某一法向应力试验确定对应的切向模量：G^a，获得如下方程：

$\mathrm{\α}({{\mathrm{\σ}}_n}^a-{\mathrm{\σ}}_n^{crit}){({{\mathrm{\σ}}_n}^a/{\mathrm{\σ}}_n^{crit})}^{k_n}=G^a---\left(7.15\right)$

则由方程(7.14、7.15)，决定常系数：

$k_n=\frac{ln\left(G^{max}\right({{\mathrm{\σ}}_n}^{\mathrm{\α}}-{{\mathrm{\σ}}_n}^{crit})/(G^a\left({{\mathrm{\σ}}_n}^{\max }-{{\mathrm{\σ}}_n}^{crit}\right))}{ln({{\mathrm{\σ}}_n}^{\max }/{{\mathrm{\σ}}_n}^{\mathrm{\α}})}$

和 $\mathrm{\α}=G^{max}/\left(\right({{\mathrm{\σ}}_n}^{\max }-{\mathrm{\σ}}_n^{crit}\left){\left({{\mathrm{\σ}}_n}^{\max }-{\mathrm{\σ}}_n^{crit}\right)}^{k_n}\right)---\left(7.16\right)$

在确定某一法向应力σ_n条件下的峰值应力切向模量G_t之后，按方程(7.12)决定a’，按方程(7.10)决定b’，至此新的邓肯--张模型各参数得以确定；

(7.3.2)方法二

取方程(7.1)的ξ＝-1，则方程变为：

$\mathrm{\τ}=\frac{G\mathrm{\γ}}{1+{\mathrm{\γ}}^q/p}---\left(7.17\right)$

在峰值应力下，则有：

${\mathrm{\τ}}^{peak}/{\mathrm{\γ}}_{peak}=\frac{G}{1+{\left({\mathrm{\γ}}_{peak}\right)}^q/p}---\left(7.18\right)$

${{\mathrm{\γ}}_{peak}}^q/p=\frac{G}{K_{scant}}-1---\left(7.19\right)$

同理对方程(7.17)求导，相应的导数为切线模量：

$\frac{\∂\mathrm{\τ}}{\∂\mathrm{\γ}}=\frac{G}{(1+{\mathrm{\γ}}^q/p)}-\frac{{\mathrm{Gq\γ}}^q/p}{{(1+{\mathrm{\γ}}^q/p)}^2}---\left(7.20\right)$

当峰值应力满足目前摩尔库伦准则时，则峰值应变也满足方程(7.4)，在峰值应力时，其切线模量记为G_t；

利用方程(7.18、7.19)，在峰值应力时，切线模量满足方程：

$G_t=K_{scant}\[1+\frac{{\mathrm{qK}}_{scant}}{G}(1-\frac{G}{K_{scant}})\]---\left(7.21\right)$

对于峰值应力对应的切线模量，按照方程(7.13)求解；利用方程(7.21)求解参数q，利用方程(7.19)求解p。

3.如权利要求2所述的计算方法，其特征在于：针对条分法，潜在滑动面计算方法包括如下子步骤，

(1)利用现行方法进行条块划分；

(2)竖向应力以比重与高度之积加以计算，水平应力和剪应力以剩余推力为水平方向和垂直于水平方向力的矢量和，水平方向力和垂直于水平方向力假设满足一定的应力分布条件；

(3)计算条块底边摩阻应力。

4.如权利要求3所述的计算方法，其特征在于：条分法决定潜在滑动面，按强度折减法进行，具体为对条块底边临界抗剪强度进行折减，直至位于临空面的破坏条块处于极限平衡状态。

5.如权利要求3所述的计算方法，其特征在于：条分法决定潜在滑动面，按施加荷载或位移工况方法进行，具体为对可能破坏施加相应的荷载或位移工况，直至位于临空面的破坏条块处于极限平衡状态。