CN110750860B - 一种土坡滑坡全过程分析方法 - Google Patents

Abstract

本发明提供了一种土坡滑坡全过程分析方法。所述方法在土坡处于小变形状态时采用理想弹塑性固体模型从而可以适用于稳定性分析，而在土坡失稳产生滑坡后自动转化为宾汉姆流体模型从而可以适用于滑坡大变形分析。从而实现土坡从稳定‑临界失稳‑失稳滑动‑最终淤积的全过程分析。

Description

一种土坡滑坡全过程分析方法

技术领域

本发明涉及计算力学领域，尤其涉及一种基于伪弹簧固液转化SPH方法的土坡滑坡全过程分析方法。

背景技术

近年来，滑坡造成的重大事故屡见不鲜，如2015年12月深圳滑坡、2017年6月24号四川阿坝州茂县滑坡等都造成了大量的人员伤亡和财产损失。然而，目前对于土坡稳定性及滑坡分析的研究，往往是将滑坡灾害前的机理研究即边坡稳定性研究，以及滑坡后的运动过程研究独立开来的。然而，土坡从稳定到滑坡结束的整个过程，即：稳定状态——临界稳定状态——失稳——滑坡运动——形成最终堆积形态的全过程，是一个密不可分的整体，因此将滑坡前的稳定性研究与滑坡后的运动过程研究分割开来是与实际灾害情况不符的。目前分开研究的原因，是由于土坡滑坡的运动过程及影响范围研究具有难以实验的特性，非常依赖于计算机辅助技术，而传统的数值模拟方法，诸如极限平衡法(LEM)、有限元法(FEM)、有限体积法(FVM)等，仅能获取边坡的安全系数及潜在滑移面，即进行边坡的稳定性分析，而对于其失稳后的滑坡运动过程分析存在一定的困难。另外，对于目前逐渐发展起来的诸如离散元法(DEM)、无网格伽辽金法(EFG)等无网格方法虽然可以进行大变形分析，但也由于其发展较晚，且存在诸如难以直接描述应力—应变关系、难以确定非物理参数等缺陷，使其应用受到了一定限制。因此，构建一种新型的数值方法，对土坡从稳定到滑坡结束的全过程进行分析，是具有十分重要的理论及实际意义的。

光滑粒子流体动力学方法(Somoothed Paticle Hydrodynamics，简称SPH)是一种纯拉格朗日无网格计算方法，近年来在岩土力学领域中开始得以成功的应用，其优点是既可以分析土坡的稳定性系数及潜在滑移面，又可以对土坡失稳后的滑坡运动过程进行研究，因此特别适合求解土坡稳定性及滑坡问题。然而，目前SPH方法在土坡稳定性及滑坡分析研究领域内仍存在不足，且较为热点的方法主要有两种：一种是将土坡视为弹塑性体，该方法在进行边坡滑坡前的稳定性分析时具有较高的准确性，但由于边坡失稳后的高速运动属于大变形过程，而继续用基于小变形理论的弹塑性本构模型计算该过程，其适用性尚需进一步探讨；另一种是将土坡视为一种非牛顿流体，该法较为适用于边坡失稳后的滑坡运动过程分析，但由于采用的是流体模型，因此无法分析边坡失稳前的稳定性系数及潜在滑移面，而需要人为指定滑移面及滑坡区域。因此，两种分析方法虽然各有其优势，但由于其选取的本构模型的特点使其应用范围也受到了一定的限制，难以真实的表达边坡从稳定状态——临界稳定状态——失稳——滑坡运动——到形成最终堆积形态的全过程。综上所述，由于目前的SPH方法在进行土坡滑坡灾害全过程分析时尚具有局限性，因此寻求一种合适的方法，即可以准确的计算土坡的稳定性安全系数及潜在滑移面，又能够准确预测土坡失稳后的滑坡运动过程是非常有必要的。

发明内容

为解决上述问题，本发明公开了一种土坡滑坡全过程分析方法，包括如下步骤：

根据所需研究土坡的截面信息，生成粒子模型，根据所述粒子模型的物理力学参数对其不同位置的相关信息分别赋值；

确定计算中需要采用的数值处理技术，并设置屈服准则、流动法则、时间步长；

进行粒子搜索；

对密度求解；

对土坡分别采用弹塑性本构模型和宾汉姆本构模型进行应力应变的求解:首先将土体视为弹塑性体，采用弹塑性本构得到第一应力应变结果，其次将土体视为非牛顿流体，采用宾汉姆模型得到第二应力应变结果；

根据得到的第一应力应变结果和第二应力应变结果,判定边坡土体所处状态，并通过弹塑性固体和宾汉姆流体求解应力应变关系及内力引起的第一速度变化率；

进行人工粘度、人工应力、速度修正以及在外力作用下引起的第二速度变化率的求解；

根据第一速度变化率及第二速度变化率之和更新质点信息，并进行应力调整，根据更新后的质点信息重新执行粒子搜索，并重新求解密度、第一应力应变结果、第二应力应变结果、应力应变关系、第一速度变化率、第二速度变化率；

根据求解所得的数据得到边坡安全系数及潜在滑移面，如果计算所得的边坡安全系数小于国家规范所规定的最低标准则将折减系数设置为国家规范值后按照前述步骤重新进土坡滑坡全过程分析；若安全系数高于国家规范规定值，则停止计算。

进一步的，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法。

进一步的，光滑核函数的计算式子为：

其中W(R,h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值,α_d在二维空间的计算公式为：α_d＝15/(7πh²)。

进一步的，应力调整方法分为两种情况：

第一种情况：当边坡土体的应力状态超出屈服面的顶点时，在D-P屈服准则下应满足以下条件:

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数；

为应力张量第一不变量，n为边坡土体的应力状态位于屈服面上时的时间步数；

此时对应力分量进行以下转换：

其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量；

第二种情况：当边坡土体应力状态位于屈服面上时，在D-P屈服准则下应满足以下条件：

其中J₂为应力偏量第二不变量；

调整的方法是在应力张量第一不变量不变的条件下，使应力偏量第二不变量减小到屈服面位置，其调整方法如下：

其中，其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量,

分别代表对应应力调整点调整前在XY、YZ、XZ方向上的剪应力分量；

rⁿ为比例系数，其表达式如下：

进一步的，D-P屈服准则为内切圆屈服准则或非正交匹配圆屈服准则。

进一步的，内切圆屈服准则屈服准则中,α_φ及k_c取值为：

非正交匹配圆屈服准则中α_φ及k_c取值为：

其中

代表边坡土体的摩擦角，c代表粘聚力。

进一步的，所述对边坡土体采用弹塑性本构模型进行应力应变的求解，得到第一应力应变结果具体是采用Drucker–Prager模型进行第一应力应变结果的求解,具体为：

其中，α和β表示笛卡尔坐标系下的坐标方向，v表示速度，v^α,v^β表示在α、β坐标方向上的速度分量，ε^αβ为应变张量，δ^αβ为狄克拉函数；

第一应力应变结果的公式为：

对于关联性流动法则有：

其中t代表时间，σ^αβ为应力张量，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量，ε^γγ为三向正应变之和，ω^αβ为自旋速率张量，且有：

而对于非关联性流动法则有：

其中

进一步的，当土体处非牛顿流体状态时，采用宾汉姆流体模型进行第二应力应变结果的求解，具体为：

其中

其中τ为第二应力应变结果，η′为等效粘度系数，

代表边坡土体的摩擦角，η₀为粘度系数；τ_y为屈服强度；γ是等效剪应变率；

为相邻粒子i,j间的偏剪切应变率。

进一步的，根据得到的第一应力应变结果和第二应力应变结果,判定边坡土体所处状态，并通过弹塑性固体和宾汉姆流体求解应力应变关系及第一速度变化率是采用伪弹簧理论，具体分为3种情况：

第一种情况：当伪弹簧中的损伤参数D_ij＝0，即接触参数f_ij＝1时，土体处于纯弹塑性状态，采用弹塑性本构计算应力应变关系及第一速度变化率；

第二种情况：当伪弹簧中的损伤参数0＜D_ij＜1，接触参数0＜f_ij＜1时，土体处于弹塑性固体和宾汉姆流体同时作用的状态，第一速度变化率由弹塑性固体和宾汉姆流体两部分计算组成：

其中

和

分别代表弹塑性本构和宾汉姆本构计算得到的速度变化率，且根据动量守恒方程在SPH下的离散形式进行计算；

其中，ρ代表密度，P代表各项同性压力，f^α代表由外力即重力引起的速度变化率,N代表粒子总数，m代表质量,τ^αβ为粘性应力张量，

为核函数的一阶导数；

第三种情况：当伪弹簧中的损伤参数D_ij＝1，接触参数f_ij＝0时，土体处于纯宾汉姆流体状态，采用宾汉姆本构计算应力应变关系及第一速度变化率；

前述三种情况中，损伤参数D_ij和f_ij的计算公式为：

其中

定义为塑性区贯通时，在滑移面上的最大等效塑性应变，Di和Dj分别为相邻两粒子的损伤参数，Dmax定义为塑性区贯通时，滑移面上的最大损伤参数。

进一步的，时间步长的选取方法根据下公式计算：

其中，c为待分析土坡材料的数值声速，h_i为粒子间距，C_cour为柯朗系数，C_cour的取值范围为0.3-0.4。

进一步的，所述根据求解所得的数据得到边坡安全系数及潜在滑移面是采用塑性区贯通或特征点位移不稳定失稳判据得到边坡安全系数及潜在滑移面。

本发明的有益效果为：

1、与目前常用的数值模拟方法，如极限平衡发、有限元法、有限差分等方法相比，SPH方法做为一种纯拉格朗日无网格方法，由于不存在网格畸变的问题特别适合处理大变形问题，本发明除了可以计算边坡的安全系数及潜在滑移面，还能够对边坡失稳后的滑坡过程进行分析，获取边坡失稳滑坡后的影响范围、淹没深度及冲击力等指标。

2、与目前已有的SPH方法主要进行滑坡分析不同，该发明在边坡处于小变形阶段(稳定阶段)采用弹塑性本构模型进行计算，而在边坡处于大变形阶段(失稳滑坡阶段)采用宾汉姆本构模型进行计算，且两种模型是通过一种“伪弹簧”方法进行连接的。

3、在边坡处于稳定状态时，可通过该法对边坡进行失稳前的稳定性分析，有效的获取边坡的安全系数及潜在滑移面，并将安全系数与国家规范值进行对比。若安全系数小于国家规范值，则将折减系数设置为国家规范值，计算在该折减系数下边坡可能的滑坡区域及影响范围，从而进一步为该边坡的风险评估及综合治理提供依据。

附图说明

图1为伪弹簧示意图。

图2为本发明所述方法的流程图。

具体实施方式

下面先对本发明所述的伪弹簧理论进行说明。

“伪弹簧”是一种简单而直观的方法，最初用于模拟SPH框架中任意演化的裂纹(或材料分离)，如不连续富集或粒子分裂。在“伪弹簧”SPH中，相互作用仅限于近邻粒子之间的相互作用，这与标准的SPH局部交互作用发生在整个影响域(由平滑长度h定义)是不同的，其目的是根据粒子的变形状态来研究粒子间的相互作用。在“伪弹簧”算法中，这种相互作用是通过一个阶数参数D如图1所示。(此处的D相当于后述的损伤参数D_ij表征的是粒子i、j之间相互作用的一种参数，用于确定土坡是弹塑型本构还是宾汉姆本构)。图1中间的图中，右边图的断裂键在处于0<D<1时，处于断裂与非断裂之间。

然而，在土坡稳定性及滑坡分析中，为了实现基于弹塑性—宾汉姆本构的固液转化，有必要对传统“伪弹簧”算法进行改进，使其“弹簧断裂”变为“弹簧变异”，并通过阶数参数D来控制弹塑性本构及非牛顿流体(即宾汉姆)本构对粒子的影响，使得粒子从弹塑性状态逐渐转变为非牛顿流体状态，从而实现在小变形分析即稳定性分析时采用弹塑性模型，而大变形分析时采用非牛顿流体模型，达到对土坡滑坡灾害进行全过程分析的目的。

传统的“伪弹簧”算法中以应变为关键变量，并且在土体产生应变的初始阶段就触发“伪弹簧”模型，而在土坡稳定性及滑坡分析中这是不适用的。一方面，并非边坡部分产生塑性变形就会导致边坡整体失稳，而需要从坡底到坡顶产生塑性变形贯通面；另一方面，土体并非产生塑性应变后迅速由固态向液态转变，而是当其塑性应变大于一定数值后(定义该值为ε_g)，才会产生失稳及固液变化。针对上述情况，本发明拟首先编制相应的塑性区贯通判定程序，判定边坡是否具备失稳的条件；其次，在上述基础上通过实验方法分析边坡失稳时滑移面上累积塑性应变的失稳“阈值”ε_g，当边坡产生塑性变形贯通面，且贯通面上的累积塑性均达到该的“阈值”(即ε_g)后，才会产生失稳并触发“伪弹簧”模型；再次，对“伪弹簧”理论进行改进，采用累计等效塑性应变值作为“伪弹簧”模型是否触发的关键变量，并以塑性区产生贯通滑移面时的最大累积塑性应变值作为“伪弹簧”模型的触发阈值(定义该值为ε_max)，即由固态向液态转化的起始值。最后，拟进行多次不同物理参数(如粘聚力、摩擦角、模型尺寸等)的室内滑坡实验，并针对每次物理实验进行多次不同数值参数(如不同的粒子大小、不同的光滑长度等)计算，得到判定塑性区贯通并失稳时的等效塑性应变“阈值”ε_g的经验公式，并最终通过将数值计算结果与室内物理实验的结果进行对比的方式，验证算法及参数选取的有效性与准确性。

用伪弹簧理论对图1进行类比解释。以左下45度的弹簧为例，弹簧两端的小圆形分别代表粒子i、j，土坡处于弹塑性状态时两者之间相互作用完好，则损伤参数D_ij＝0，接触参数f_ij＝1，当土坡处于弹塑性状态向非牛顿流体状态转变时，粒子粒子i、j之间相互作用出现变化时，则损伤参数0＜D_ij＜1，接触参数0＜f_ij＜1，当土坡处于非牛顿流体状态时D_ij＝1，接触参数f_ij＝0。

基于前述的伪弹簧理论，下面对本申请的技术方案进行说明。

本发明的设计构思为：本发明基于光滑粒子流体动力学(SPH)提出一种适用于土质边坡稳定性及滑坡分析的新方法，主要用于解决的技术问题有：

1、目前常用数值模拟方法仅能进行稳定性分析，而对滑坡运动过程分析存在困难的问题。

2、针对目前现有的SPH方法采用单一的本构模型进行稳定性计算或滑坡计算的问题，采用基于伪弹簧的SPH方法，在小变形范围内采用弹塑性本构模型，大变形范围采用宾汉姆本构模型，从而解决目前SPH方法在进行滑坡全过程分析时的不足。

具体的，如图2所示，本发明所述的土坡滑坡全过程分析方法，包括如下步骤：

步骤一：根据所需研究土坡的截面信息，生成粒子模型，根据所述粒子模型的物理力学参数对其不同位置的相关信息分别赋值。

所述相关信息包括密度、弹性模量、泊松比、外力、土体屈服前的粘聚力和内摩擦角，以及土体屈服后的粘度系数、粘聚力、和内摩擦角。

步骤二：确定计算中需要采用的数值处理技术，并设置屈服准则、流动法则、时间步长。

所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法。

目前常用的密度计算方法主要有连续密度法、密度求和法等；光滑核函数主要有三次样条核函数、四次样条核函数、五次样条核函数等；粒子搜索方法主要有全配对搜索方法、链表式搜索方法、树形搜索方法等；固壁边界处理方法主要有边界排斥力法、镜像粒子法、虚拟粒子法等；时间积分方法主要有verlet算法、leap-frog法、predictor-corrector法等。各个技术可采用现有方法。下面着重介绍以下技术。

1.光滑核函数

目前常用的有三次样条核函数、四次样条核函数、五次样条核函数等，本实施例采用三次样条核函数进行计算：

其中W(R,h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值，α_d在二维三维空间的值为α_d＝15/(7πh²)。

2.应力调整方法

在计算过程中当产生塑性变形时可能会出现应力状态位于屈服面之外的情况，这种情况明显与实际情况不符，应力调整技术就是用来处理这种特殊情况的。

第一种情况是土体的应力状态在屈服面的顶点之外，这种错误可称为拉伸裂缝，这种情况与有限元中的情况类似，而在SPH中可能会产生非物理裂缝或粒子集聚，处理的方法是通过调整应力分量将应力状态回归到屈服面的顶点处，当边坡土体的应力状态超出屈服面的顶点时，在D-P屈服准则下应满足以下条件:

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数(当选择不同的匹配准则的时候其取值不同)；

为应力张量第一不变量，n为边坡土体的应力状态超出屈服面的顶点时的时间步数。在第一种错误情况下，调整的办法是使上式等于0而非小于0，则可对应力分量进行以下转换：

其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量，等式左边为调整后的对应应力分量。

第二种情况是当土体产生塑性变形时，真实的土体应力状态应该是始终位于屈服面上，而非屈服面的上方。然而在数值计算过程中，可能会出现应力状态远离屈服面的现象，在这种情况下，应采取相应的数值方法使应力状态回归到屈服的面上。当材料的粒子在第n个时间步的应力状态超过屈服面时，在D-P屈服准则下应满足以下条件：

其中J₂为应力偏量第二不变量,n为边坡土体应力状态位于屈服面上时的时间步数。

调整的方法是在应力张量第一不变量不变的条件下，使应力偏量第二不变量减小到屈服面位置，其调整方法如下：

其中，其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量,

分别代表对应应力调整点调整前在XY、YZ、XZ方向上的剪应力分量，等式左边为调整后的对应应力分量。

rⁿ为比例系数，其表达式如下：

3.屈服准则

本实施例采用D-P屈服准则，该准则共有5种形式：外接外角圆(D-P1)、内角外接圆(D-P2)、内切圆(D-P3)、等面积圆(D-P4)、非正交匹配圆(D-P5)。更优选的，本实施例采用内切圆屈服准则或非正交匹配圆屈服准则，此时进行稳定性及滑坡分析时得到的效果最好。

内切圆屈服准则屈服准则中,α_φ及k_c取值为：

非正交匹配圆屈服准则中α_φ及k_c取值为：

其中

代表边坡土体的摩擦角，c代表粘聚力。

3.时间步长

时间步长根据模型的粒子间距等进行设置

时间步长的选取方法根据下公式计算：

其中，c为待分析土坡材料的数值声速，h_i为粒子间距，C_cour为柯朗系数，一般取0.2，本实施例C_cour优选取0.3-0.4，当为该取值范围时，仍可以保持计算稳定，并且可以减少计算时间。

步骤三：进行粒子搜索。

本步骤还需要根据预先设置的光滑核函数的光滑半径确定邻近粒子形成粒子对。

步骤四：对密度求解。

根据步骤二选择的密度求解方法进行密度求解，随着时间步长更新粒子密度。

步骤五：对土坡分别采用弹塑性本构模型和宾汉姆本构模型进行应力应变的求解:首先将土体视为弹塑性体，采用弹塑性本构得到第一应力应变结果，其次将土体视为非牛顿流体，采用宾汉姆模型得到第二应力应变结果。

用Drucker–Prager模型进行第一应力应变结果的求解：

其中，α和β表示坐标方向，v表示速度，v^α,v^β表示在α、β坐标方向上的速度分量，ε^αβ为应变张量，δ^αβ为狄克拉函数；

第一应力应变结果的公式为：

对于关联性流动法则有：

等式左边为内力的应力张量随时间的变化率，t代表时间，σ^αβ为应力张量，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量，ε^γγ为三向正应变之和，ω^αβ为自旋速率张量，且有：

而对于非关联性流动法则有：

其中

参数解释如前所述。

当土体视为非牛顿流体时，采用宾汉姆流体模型进行第二应力应变结果的求解具体为：

其中η₀为粘度系数；τ_y屈服强度；γ是等效剪应变率；

为相邻粒子i,j间的偏剪切应变率；采用等效粘度公式，将宾汉姆流体的计算方法与摩尔库伦准则结合，得到最终的第二应力应变结果的计算公式为：

其中等式左边为第二应力应变结果，即粘性应力；

其中η′为等效粘度系数，。

步骤六：根据得到的第一和第二应力应变结果,采用伪弹簧理论判定边坡土体所处状态，并通过弹塑性固体及宾汉姆流体求解内力引起的第一速度变化率。

对于内力及速度变化率的计算分为三种情况：

第一种情况：当伪弹簧中的损伤参数D_ij＝0，即接触参数f_ij＝1时，土体处于纯弹塑性状态。

第二种情况：当伪弹簧中的损伤参数0＜D_ij＜1，接触参数0＜f_ij＜1时，土体处于弹塑性固体和宾汉姆流体同时作用的状态，土体的速度变化率由弹塑性固体和宾汉姆流体两部分计算组成：

其中

和

分别代表弹塑性本构和宾汉姆本构计算得到的速度变化率，且根据动量守恒方程在SPH下的离散形式进行计算。

其中，ρ代表密度，P代表静压，f^α代表由外力即重力引起的速度变化率。

第三种情况：当伪弹簧中的损伤参数D_ij＝1，接触参数f_ij＝0时，土体处于纯宾汉姆流体状态，采用宾汉姆本构计算应力应变关系及速度变化率。

D_ij和f_ij的计算公式为：

其中ε_p为等效塑性应变。

步骤七：进行人工粘度、人工应力、速度修正以及在外力作用下引起的第二速度变化率的求解。

步骤八：根据第一速度变化率及第二速度变化率之和更新质点信息(如质点的密度、质量、速度等)，并进行应力调整，根据更新后的质点信息重新执行粒子搜索，并重新执行前述步骤，包括求解密度、第一应力应变结果、第二应力应变结果、应力应变关系、第一速度变化率、第二速度变化率。

步骤九：根据求解所得的数据进行处理好分析，采用“塑性区贯通”或“特征点位移不稳定”失稳判据，得到边坡安全系数及潜在滑移面。该步骤为边坡的安全评估及综合治理提供可靠保障。当采用“塑性区贯通为失稳判据时”，分别采用不同的折减系数进行单独计算，观察其塑性区范围随折减系数的变化规律，如果在某一折减系数下塑性区从坡底到坡顶达到贯通，则该折减系数便为该边坡的安全系数，对应的塑性区为边坡的潜在滑移面；当采用“特征点位移不稳定为失稳判据时”，首先使得折减系数随时间步的增加而增加，通过观察特征点位移的变化情况获取安全系数的一个大致范围，其次再固定折减系数，观察特征点位移随着时间步的变化情况，如果在某一折减系数条件下，特征点位移随着时间步的增加而不断增大，则认为在该折减系数下边坡处于不稳定的状态，而边坡从稳定到不稳定状态过读的临界折减系数，便是该边坡的安全系数。

如果计算所得的边坡安全系数小于国家规范所规定的最低标准则返回步骤一继续对该边坡进行滑坡影响范围分析，做法是将折减系数设置为国家规范值，计算在该折减系数下边坡可能的滑坡区域及影响范围，从而进一步为该边坡的风险评估及综合治理提供依据。若安全系数高于国家规范规定值，则停止计算。

以上所述的，仅为本发明的具体实施方式，但本发明的保护范围并不局限于上述实施例的限制，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述的权利要求的保护范围为准。

Claims (9)

1.一种土坡滑坡全过程分析方法，其特征在于,包括如下步骤：

根据所需研究土坡的截面信息，生成粒子模型，根据所述粒子模型的物理力学参数对其不同位置的相关信息分别赋值；

确定计算中需要采用的数值处理技术，并设置屈服准则、流动法则、时间步长；

进行粒子搜索；

对密度求解；

对土坡分别采用弹塑性本构模型和宾汉姆本构模型进行应力应变的求解:首先将土体视为弹塑性体，采用弹塑性本构得到第一应力应变结果，其次将土体视为非牛顿流体，采用宾汉姆模型得到第二应力应变结果；

根据得到的第一应力应变结果和第二应力应变结果,判定边坡土体所处状态，并通过弹塑性固体和宾汉姆流体求解应力应变关系及内力引起的第一速度变化率；

进行人工粘度、人工应力、速度修正以及在外力作用下引起的第二速度变化率的求解；

根据第一速度变化率及第二速度变化率之和更新质点信息，并进行应力调整，根据更新后的质点信息重新执行粒子搜索，并重新求解密度、第一应力应变结果、第二应力应变结果、应力应变关系、第一速度变化率、第二速度变化率；

根据求解所得的数据得到边坡安全系数及潜在滑移面，如果计算所得的边坡的安全系数小于国家规范所规定的最低标准则将折减系数设置为国家规范值后按照前述步骤重新进土坡滑坡全过程分析；若安全系数高于国家规范规定值，则停止计算；

根据得到的第一应力应变结果和第二应力应变结果,判定边坡土体所处状态，并通过弹塑性固体和宾汉姆流体求解应力应变关系及第一速度变化率是采用伪弹簧理论，具体分为3种情况：

第一种情况：当伪弹簧中的损伤参数D_ij＝0，即接触参数f_ij＝1时，土体处于纯弹塑性状态，采用弹塑性本构计算应力应变关系及第一速度变化率；

第二种情况：当伪弹簧中的损伤参数0＜D _ij＜1，接触参数0＜f _ij＜1时，土体处于弹塑性固体和宾汉姆流体同时作用的状态，第一速度变化率由弹塑性固体和宾汉姆流体两部分计算组成：

其中

和

分别代表弹塑性本构和宾汉姆本构计算得到的速度变化率，且根据动量守恒方程在SPH下的离散形式进行计算；

其中，ρ代表密度，P代表各项同性压力，f^α代表由外力即重力引起的速度变化率,N代表粒子总数，m代表质量,τ^αβ为粘性应力张量，

为核函数的一阶导数；

第三种情况：当伪弹簧中的损伤参数D_ij＝1，接触参数f_ij＝0时，土体处于纯宾汉姆流体状态，采用宾汉姆本构计算应力应变关系及第一速度变化率；

前述三种情况中，损伤参数D_ij和f_ij的计算公式为：

其中

定义为塑性区贯通时，在滑移面上的最大等效塑性应变，Di和Dj分别为相邻两粒子的损伤参数，Dmax定义为塑性区贯通时，滑移面上的最大损伤参数。

2.根据权利要求1所述的一种土坡滑坡全过程分析方法，其特征在于:所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法。

3.根据权利要求2所述的一种土坡滑坡全过程分析方法，其特征在于:光滑核函数的计算式子为：

其中W(R,h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值,α_d在二维空间的计算公式为：α_d＝15/(7πh²)。

4.根据权利要求2所述的土坡滑坡全过程分析方法，其特征在于:应力调整方法分为两种情况：

第一种情况：当边坡土体的应力状态超出屈服面的顶点时，在D-P屈服准则下应满足以下条件:

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数；

为应力张量第一不变量，n为边坡土体的应力状态位于屈服面上时的时间步数；

此时对应力分量进行以下转换：

其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量；

第二种情况：当边坡土体应力状态位于屈服面上时，在D-P屈服准则下应满足以下条件：

其中J₂为应力偏量第二不变量；

调整的方法是在应力张量第一不变量不变的条件下，使应力偏量第二不变量减小到屈服面位置，其调整方法如下：

其中，其中

代表对应应力调整点调整前在X、Y、Z方向上的应力分量,

分别代表对应应力调整点调整前在XY、YZ、XZ方向上的剪应力分量；

rⁿ为比例系数，其表达式如下：

5.根据权利要求4所述的土坡滑坡全过程分析方法，其特征在于，D-P屈服准则为内切圆屈服准则或非正交匹配圆屈服准则。

6.根据权利要求5所述的土坡滑坡全过程分析方法，其特征在于,

内切圆屈服准则屈服准则中,α_φ及k_c取值为：

非正交匹配圆屈服准则中α_φ及k_c取值为：

其中

代表边坡土体的摩擦角，c代表粘聚力。

7.根据权利要求1～6任一项所述的土坡滑坡全过程分析方法，其特征在于:对边坡土体采用弹塑性本构模型进行应力应变的求解，得到第一应力应变结果具体是采用Drucker–Prager模型进行第一应力应变结果的求解,具体为：

其中，α和β表示笛卡尔坐标系下的坐标方向，v表示速度，v^α,v^β表示在α、β坐标方向上的速度分量，ε^αβ为应变张量，δ^αβ为狄克拉函数；

第一应力应变结果的公式为：

对于关联性流动法则有：

其中t代表时间，σ^αβ为应力张量，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量，ε^γγ为三向正应变之和，ω^αβ为自旋速率张量，且有：

而对于非关联性流动法则有：

其中

8.如权利要求1所述的土坡滑坡全过程分析方法，其特征在于，当土体处非牛顿流体状态时，采用宾汉姆流体模型进行第二应力应变结果的求解，具体为：

其中

其中τ为第二应力应变结果，η′为等效粘度系数，

代表边坡土体的摩擦角，η₀为粘度系数；τ_y为屈服强度；γ是等效剪应变率；

为相邻粒子i,j间的偏剪切应变率。

9.如权利要求1所述的土坡滑坡全过程分析方法，其特征在于，时间步长的选取方法根据下公式计算：

其中，c为待分析土坡材料的数值声速，h_i为粒子间距，C_cour为柯朗系数，C_cour的取值范围为0.3-0.4。

Images