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

技术特征：

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

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

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

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

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

${\mathrm{\σ}}_i-{\mathrm{\σ}}_3=E_{ii}{\mathrm{\ϵ}}_i{(1+{{\mathrm{\ϵ}}_i}^{m_i}/S_i)}^{{\mathrm{\ρ}}_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：实施边坡渐进破坏的数值计算。