CN108334719B - 一种基于sph法的土质边坡稳定性及滑坡运动过程分析方法 - Google Patents

Abstract

本发明公开了一种基于SPH法的土质边坡稳定性及滑坡运动过程分析方法。本发明可以计算边坡的安全系数及潜在滑移面，还能够对边坡失稳后的滑坡过程进行分析，获取边坡失稳滑坡后的影响范围、淹没深度及冲击力等指标。该发明可对边坡进行失稳前的稳定性分析，且提出了两种失稳判据可有效的获取边坡的安全系数及潜在滑移面，可以在不事先指定滑移面（或滑移范围）的前提下，自动计算失稳滑移的区域。

Description

一种基于SPH法的土质边坡稳定性及滑坡运动过程分析方法

技术领域

本发明涉及计算力学领域，尤其涉及一种基于SPH法的土质边坡稳定性及滑坡运动过程分析方法。

背景技术

近年来，滑坡造成的重大事故屡见不鲜，如2015年12月深圳滑坡、2017年6月24号四川阿坝州茂县滑坡等都造成了大量的人员伤亡和财产损失，而造成滑坡的重要原因之一在于边坡失稳即稳定性不足，因此稳定性分析一直是岩土工程领域内的研究热点之一。

目前对于边坡的稳定性分析分析方法主要有极限平衡法、有限元法有限差分方法等，这些方法在进行边坡稳定性分析时具有较高的准确性，但边坡失稳后的滑坡过程分析却存在困难。因此，迫切需要一种有效的计算方法，既能对边坡失稳前的稳定性进行分析，又能计算边坡失稳后的滑坡运动过程，从而为边坡的安全评估及综合治理提供可靠的依据。

光滑粒子流体动力学(SPH)方法是一种纯拉格朗日无网格方法，特别适合求解大变形领域的问题。但目前SPH在岩土工程领域内主要还是应用在滑坡后影响范围分析，且主要进行已知滑移面、已产生事故的边坡反演分析。与目前已有的SPH方法主要进行滑坡分析不同，该发明可对边坡进行失稳前的稳定性分析，且提出了两种失稳判据可有效的获取边坡的安全系数及潜在滑移面；另外，现有SPH技术在进行滑坡分析时需要事先指定滑移范围，而该发明可以在不事先指定滑移面(或滑移范围)的前提下，自动计算失稳滑移的区域及影响范围。

发明内容

本发明提供了一种基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，包括如下步骤：

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

步骤二：预先确定计算工程中的数值处理技术，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法；

步骤三：设置屈服准则、流动法则、时间步长；

步骤四：进行粒子搜索；

步骤五：进行密度求解；

步骤六：进行应力应变的求解；

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

步骤八：根据总的速度变化率更新质点信息，并进行应力调整，更新后的信息进入下一循环按以上步骤进行循环计算，循环完成后进行步骤九；

步骤九：对输出文件进行处理和分析，得到边坡安全系数及潜在滑移面。

步骤十：将折减系数设置为国家规范值，返回步骤一，计算在该折减系数下边坡可能的滑坡区域及影响范围。

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

其中W(R,h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值，α_d在一维、二维、三维空间下分别为

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

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

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数，当选择不同的匹配准则的时候其取值不同；I₁ ⁿ为应力第一不变量；

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

其中

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

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

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

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

其中，其中

代表某一点在X、Y、Z方向上的应力分量,

分别代表在XY、YZ、XZ方向上的剪应力分量；

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

进一步的，屈服准则为内切圆屈服准则及非正交匹配圆屈服准则，

D-P3屈服准则中α_φ及k_c取值为：

D-P5屈服准则中α_φ及k_c取值为：

其中

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

进一步的，步骤六中，应变求解的公式为：

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

应力求解的公式为：

经过一系列推导和变换，对于关联性流动法则有：

其中

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

其中

t代表时间，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量。

进一步的，步骤一中所述相关信息包括密度、弹性模量、泊松比、粘聚力、内摩擦角、外力。

本发明的有益效果为：

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

2、与目前已有的SPH方法主要进行滑坡分析不同，该发明可对边坡进行失稳前的稳定性分析，且提出了两种失稳判据可有效的获取边坡的安全系数及潜在滑移面；另外，现有SPH技术在进行滑坡分析时需要事先指定滑移范围，而该发明可以在不事先指定滑移面(或滑移范围)的前提下，自动计算失稳滑移的区域。

3、在进行稳定性分析时，由于SPH方法与传统数值模拟方法的差异，其安全系数的判定方式也不同。本发明针对SPH方法在边坡稳定性分析中的应用，提出2种安全系数的判定方式，可准确可靠的得到边坡的安全系数及潜在的滑移面。

4、在进行稳定性分析的基础上，可对安全系数过小的边坡进一步进行滑坡运动分析，做法是将折减系数设置为国家规范值，计算在该折减系数下边坡可能的滑坡区域及影响范围，从而进一步为该边坡的风险评估及综合治理提供依据。

具体实施方式

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

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

2、针对目前现有的SPH方法主要进行事故后滑坡分析，而在稳定性分析中存在困难、滑坡分析需要事先指定滑移面的问题，提出稳定性分析的方法及适用于SPH的边坡失稳判据。

本发明包括如下步骤：

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

该步骤可以分为两种情况，一种是对于结构相对简单的边坡可通过简单的编程实现；另一种是对于截面复杂，特别是内部信息较为复杂的边坡可采用商业软件进行处理(如ABAQUS、HYPERMESH等首先进行网格划分，再提取网格节点信息代入SPH程序生成SPH粒子模型)。在生成粒子模型后，根据模型的物理力学参数对其不同位置的相关信息分别赋值。相关信息包括密度、弹性模量、泊松比、粘聚力、内摩擦角及外力等信息。

步骤二：预先确定计算工程中的数值处理技术，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法。

目前常用的密度计算方法主要有连续密度法、密度求和法等；光滑核函数主要有三次样条核函数、四次样条核函数、五次样条核函数等；粒子搜索方法主要有全配对搜索方法、链表式搜索方法、树形搜索方法等；固壁边界处理方法主要有边界排斥力法、镜像粒子法、虚拟粒子法等；时间积分方法主要有verlet算法、leap-frog法、predictor-corrector法等。各个技术可采用现有方法，对于描述上述有关方法的文献较多，这里不再一一赘述，以下仅对本文采用的一种新型的光滑四次核函数及应力调整方法进行介绍：

1.光滑四次核函数

目前常用的有三次样条核函数、四次样条核函数、五次样条核函数等，本发明在考虑了传统核函数的优势和不足的情况了，新建了一种新型的分段四次核函数，该函数可有效的解决压缩不稳定性问题，且其二次导数更加光滑，并具有较高的峰值效

其中W(R,h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值。α_d在一维、二维、三维空间下分别为

2.应力调整方法

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

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

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数(当选择不同的匹配准则的时候其取值不同)。在第一种错误情况下，调整的办法是使上式等于0而非小于0，则可对应力分量进行以下转换：

其中

代表某一点在各个方向上的应力分量。

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

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

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

其中

代表应力调整点在各个方向上的应力分量。

步骤三：设置屈服准则、流动法则、时间步长。

本发明采用的D-P屈服准则，该准则共有5种形式：外接外角圆(D-P1)、内角外接圆(D-P2)、内切圆(D-P3)、等面积圆(D-P4)、非正交匹配圆(D-P5)，经过分析，内切圆(D-P3)或非正交匹配圆(D-P5)屈服准则在进行稳定性及滑坡分析时得到的效果最好，其α和k的取值分别为：

D-P3屈服准则：

D-P5屈服准则：

其中

代表土体的摩擦角。

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

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

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

步骤四：进行粒子搜索，并根据预先设置的光滑半径确定邻近粒子形成粒子对。

根据步骤二选择的粒子搜索方法对各个待计算粒子进行搜索，并根据设置的光滑半径确定邻近粒子形成粒子对，并通过粒子对在后续步骤确定计算粒子在每个时间步后的相关物理量。

步骤五：进行密度求解。

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

步骤六：进行应力应变的求解。

首先通过粒子间的相对速度求得应变张量，其次再结合弹塑性本构模型理论、D-P屈服准则理论、关联或非关联流动法则理论求得应力及内力引起的应力变化率。

应变求解：

其中，α和β表示坐标方向，v表示速度，ε^αβ为应张量，δ^αβ为狄克拉函数。

应力求解：

经过一系列推导和变换，对于关联性流动法则有：

其中

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

其中

在公式(17)～(20)中，t代表时间，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量。

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

步骤八：根据总的速度变化率更新质点信息，并进行应力调整，更新后的信息进入下一循环按以上步骤进行循环计算，循环完成后进行步骤九；

步骤九：对输出文件进行处理和分析，采用“塑性区贯通”或“特征点位移不稳定”失稳判据，得到边坡安全系数及潜在滑移面。

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

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

Claims (5)

1.一种基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，其特征在于，包括如下步骤：

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

步骤二：预先确定计算工程中的数值处理技术，所述数值处理技术包括密度计算方法、光滑核函数、粒子搜索方法、固壁边界计算方法、时间积分方法、应力调整方法；

光滑核函数的计算式子为：

其中W(R，h)为光滑核函数，h为光滑半径，R为粒子间距与光滑半径的比值，α_d在一维、二维、三维空间下分别为

应力调整方法分为两种情况：

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

其中α_φ、k_c分别为与内摩擦角、粘聚力有关的常数；I₁ ⁿ为应力第一不变量；

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

其中

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

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

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

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

其中，其中

代表某一点在X、Y、Z方向上的应力分量,

分别代表在XY、YZ、XZ方向上的剪应力分量；

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

步骤三：设置屈服准则、流动法则、时间步长；

步骤四：进行粒子搜索；

步骤五：进行密度求解；

步骤六：进行应力应变的求解；

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

步骤八：根据总的速度变化率更新质点信息，并进行应力调整，更新后的信息进入下一循环按以上步骤进行循环计算，循环完成后进行步骤九；

步骤九：对输出文件进行处理和分析，得到边坡安全系数及潜在滑移面，如果计算所得的安全系数小于国家规范所规定的最低标准则进行下一步；

步骤十：将折减系数设置为国家规范值，返回步骤一，计算在该折减系数下边坡可能的滑坡区域及影响范围。

2.如权利要求1所述的基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，其特征在于，步骤三中，屈服准则为内切圆屈服准则或非正交匹配圆屈服准则，

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

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

其中

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

3.如权利要求1所述的基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，其特征在于，步骤六中，应变求解的公式为：

其中，α和β表示坐标方向，v表示速度，e^αβ为应张量，δ^αβ为狄克拉函数；

应力求解的公式为：

经过一系列推导和变换，对于关联性流动法则有：

其中

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

其中

t代表时间，G表示剪切模量，K为体积模量，e^αβ为偏剪切应变率张量，S^αβ为偏剪切应力率张量。

4.如权利要求1所述的基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，其特征在于，步骤一中所述相关信息包括密度、弹性模量、泊松比、粘聚力、内摩擦角、外力。

5.如权利要求1所述的基于SPH法的土质边坡稳定性及滑坡运动过程分析方法，其特征在于，步骤九中，采用塑性区贯通或特征点位移不稳定失稳判据，对输出文件进行处理和分析，得到边坡安全系数及潜在滑移面。