基于SPH的土坡滑面分析判定方法、系统、终端及介质与流程

本发明属于建设工程模拟信息化技术领域，尤其涉及一种基于sph的土坡滑面分析判定方法、系统、终端及介质。

背景技术：

在岩土工程领域，边坡稳定性问题是最基本的问题之一，广泛存在于天然边坡、人工边坡、挡土结构、水库大坝等结构的整体及局部稳定性分析当中。正确评价边坡的稳定性对于各类工程的安全运行极具重要意义。然而，以极限平衡法为代表边坡稳定性分析方法，确定滑面的位置和形状是其核心问题，对分析结果具有决定性和控制性作用。

随着计算机性能与计算方法的不断提升与改进，以有限元为代表的数值分析方法逐步应用在岩土工程中。由于有限元方法在模拟计算大变形区域过程中会产生数值不稳定性；从初始状态到后续大变形的整个变形过程难以连续求解。所以，极限平衡状态的判断尤为重要，为了解决此类问题，计算流体动力学(cfd)建模和离散建模(如dem)等方法逐渐发展起来。在cfd中，计算网格固定而无变形，但材料要求特殊类型的流体。故cfd针对流动问题有效，但难以解决静态变形问题。此外，基于固体力学的本构模型不易处理材料在变形过程中的应力-应变历史；而dem离散建模不适合处理基于连续介质近似的土体本构模型。

光滑粒子流体动力学方法(smoothedparticlehydrodynamics，sph)是一种拉格朗日无网格计算方法，可以解决大变形问题而不产生网格畸变。目前sph方法在土坡稳定性及滑坡分析研究领域内仍存在不足，且较为热点的方法主要有两种：一种是将土坡视为弹塑性体，该方法在进行边坡滑坡前的稳定性分析时具有较高的准确性，但由于边坡失稳后的高速运动属于大变形过程，而继续用基于小变形理论的弹塑性本构模型计算该过程其适用性尚需进一步探讨。因此，目前的sph方法在进行土坡滑坡灾害全过程分析时尚具有局限性，需要寻求一种合适的方法来对土坡滑面全过程运动进行分析是有必要的。

技术实现要素：

针对现有技术中的缺陷，本发明实施例提供了一种基于sph的土坡滑面分析判定方法、系统、终端及介质，依据sph方法对以土坡从小变形到大变形连续求解确定土坡滑面的具体位置和形态，能更准确地预测在安全与危险的边界处的岩土材料的变形。

第一方面，本发明实施例提供的一种基于sph的土坡滑面分析判定方法，包括：

根据所需研究土坡的剖面信息，生成粒子模型；

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

对粒子模型进行粒子搜索及配对；

计算粒子密度；

对土坡分别采用弹塑性本构模型进行应力应变的求解得到粒子应力应变结果，在弹塑性本构模型中，岩土体的抗剪强度采用粘聚力和土体的内摩擦角表示；

根据粒子应力应变结果判断岩土体的应力状态，采用修正方法对粒子的应力状态进行调整；

根据sph边坡稳定性分析强度折减法对岩土体抗剪强度参数进行折减，得到边坡的安全系数及潜在滑面；

根据潜在滑面分析出滑面位置和形态。

第二方面，本发明实施例还提供一种基于sph的土坡滑面分析判定系统，包括：粒子模型生成模块、应力应变分析模块、修正调整模块和滑面分析模块，

所述粒子模型生成模块用于根据所需研究土坡的剖面信息，生成粒子模型；

所述应力应变分析模块用于确定计算中需要采用的数值处理技术，并设置屈服准则和时间步长；

对粒子模型进行粒子搜索，计算粒子密度，对土坡分别采用弹塑性本构模型进行应力应变的求解得到粒子应力应变结果，在弹塑性本构模型中，岩土体的抗剪强度采用粘聚力和土体的内摩擦角表示；

所述修正调整模块用于根据粒子应力应变结果判断岩土体的应力状态，采用修正方法对粒子的应力状态进行调整；

所述滑面分析模块用于根据sph边坡稳定性分析强度折减法对岩土体抗剪强度参数进行折减，得到边坡的安全系数及潜在滑面，根据潜在滑面分析出滑面位置和形态。

第三方面，本发明实施例还提供一种智能终端，包括处理器、输入设备、输出设备和存储器，所述处理器、输入设备、输出设备和存储器相互连接，所述存储器用于存储计算机程序，所述计算机程序包括程序指令，所述处理器被配置用于调用所述程序指令，执行上述实施例所述的方法。

第四方面，本发明实施例还提供一种计算机可读存储介质，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时使所述处理器执行上述实施例所述的方法。

本发明的有益效果：

本发明实施例提供的一种基于sph的土坡滑面分析判定方法，采用sph方法对土坡从小变形到大变形的连续求解确定土坡滑面的具体位置和形态，能更准确地预测在安全与危险的边界处的岩土材料的变形。

本发明实施例提供的一种基于sph的土坡滑面分析判定系统，采用sph方法对土坡从小变形到大变形的连续求解确定土坡滑面的具体位置和形态，能更准确地预测在安全与危险的边界处的岩土材料的变形。

本发明实施例提供的一种智能终端，与上述基于sph的土坡滑面分析判定方法出于相同的发明构思，具有相同的有益效果。

本发明提供的一种计算机可读存储介质，与上述基于sph的土坡滑面分析判定方法出于相同的发明构思，具有相同的有益效果。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍。在所有附图中，类似的元件或部分一般由类似的附图标记标识。附图中，各元件或部分并不一定按照实际的比例绘制。

图1为本发明第一实施例所提供的一种基于sph的土坡滑面分析判定方法的流程图；

图2为本发明第一实施例中的实例的模型；

图3为第一实施例提供的方法得到的某一个时刻滑面的示意图；

图4为第一实施例提供的方法得到的某一时刻的最大剪切应变率sph粒子位置的分布图；

图5为第一实施例中采用两种方法求解得到的滑面对比图；

图6为本发明另一实施例所提供的一种基于sph的土坡滑面分析判定系统的结构框图；

图7为本发明另一实施例所提供的一种智能终端的结构框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。

还应当理解，在此本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。

还应当进一步理解，本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。

如在本说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。

需要注意的是，除非另有说明，本申请使用的技术术语或者科学术语应当为本发明所属领域技术人员所理解的通常意义。

图1示出了本发明第一实施例所提供的一种基于sph的土坡滑面分析判定方法的流程图，该方法包括以下步骤：

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

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

s2:确定计算中需要采用的数值处理技术，并设置屈服准则和时间步长。

s3:对粒子模型进行粒子搜索及配对。

s4:计算粒子密度。

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

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

核近似场函数的公式为：

<f((1)

光滑函数w用来表示临近的β粒子对α粒子的影响。在本实施例中，光滑核函数采用三次样条核函数。光滑核函数的计算公式为：

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

另外，粒子近似将粒子的信息用光滑域内的其他粒子表示，其特征如下：

式(3)中，wij＝w(xi-xj，h),mj和ρj分别表示粒子的质量和密度，n为整数。

s5:对土坡分别采用弹塑性本构模型进行应力应变的求解得到粒子应力应变结果，在弹塑性本构模型中，岩土体的抗剪强度采用粘聚力和土体的内摩擦角表示。

本实施例采用的基于固体力学的sph方法，连续性方程和运动方程可以表示为:

其中ui是速度矢量,ρ是密度,σij是应力张量和fi是外力向量。将sph插值理论应用于状态方程，上式可以写为:

其中为人工黏性项与人工应力项，表示为：

式(8)中，δij是克罗内克δ函数，πi是人工黏性，和是人工应力张量，f^β是排斥力项，n^as是常数，可设置为4。人工黏性项与人工应力项是为了解决sph算法中数值不稳定性、拉伸不稳定性。这样可能导致误差，所以本实施例将岩土体假设为理想弹塑性材料，采用采用drucker-prager屈服准则确定岩土体的塑性流动状态。

s6:根据粒子应力应变结果判断岩土体的应力状态，采用数值处理方法对粒子的应力状态进行调整。

当岩土体的应力状态处于塑性屈服状态时，它随屈服函数的变化而变化。然而，数值误差会导致应力状态超过屈服函数。因此，模拟计算过程中需要及时纠正计算粒子的应力状态。本采用chen和mizuno提出的应力映射调整算法。如果计算粒子的应力状态在计算中出现误差，处在了拉伸状态，采用式(9)将静水压力应力分量进行调整，修正到临界拉伸状态。

式(9)中，为应力拉伸状态处理调整前的粒子应力张量，δ'ij为应力拉伸状态处理调整后的粒子应力张量，kc为粘聚力，为内摩擦角，n为边坡土体的应力状态位于屈服面上时的时间步数，i1为静水应力张量第一分量。如果在计算模拟过程中出现计算误差，粒子应力状态超出了屈服表面，采取应力调整的过程。按比例系数r减小偏切应力分量，静水应力张量第一分量

i1保持不变，如下式所示:

式(10)中，σ′ij是经粒子伸拉状态处理后的粒子应力张量；r是调整比例系数，由下式计算：

式(11)中，其中j2为应力偏量第二不变量。除了上述应力调整之外，本文还采用了jaumann应力速率来考虑刚体自旋运动的影响。αφ和kc是德鲁克-普拉格的常量,两者可以与摩尔库伦材料参数粘聚力c与内摩擦角φ相互转换求得。在平面应变状态下，德鲁克-普拉格的常量定义为：

s7:根据sph边坡稳定性分析强度折减法对岩土体抗剪强度参数进行折减，得到边坡的安全系数及潜在滑面。

根据边坡稳定的极限平衡方法，安全系数的定义为边坡抗滑力与下滑力的比值，即滑面的抗剪强度与滑面实际剪应力的比值。因此，以弹塑性本构岩土体而言，在极限平衡状态下，滑面上任意位置土体的剪应力都达到了其屈服强度，均处于塑性破坏的临界状态。

与基于有限元法或有限差分法的强度折减法相似，基于sph的边坡稳定性分析强度折减法通过对岩土体抗剪强度参数c和进行一系列折减得到ct和其中，折减方法如式(12)、式(13)所示。

其中，c和为岩土体实际抗剪强度参数；c表示为粘聚力，表示内摩擦角，ct为强度折减后的黏聚力值；为强度折减后内摩擦角值；fost为试算的安全系数(折减倍数)。通常，初始时刻fost设置得足够小，以便系统稳定。随后，fost值逐渐增加，直到出现贯穿的滑面并持续呈现，边坡失稳。使边坡出现临界失稳状态的fost的最终值被定义为边坡的安全系数。

基于有限元法或有限差分法的强度折减法不同在确定边坡的极限平衡状态时有多种判别标准，如计算收敛、塑性区贯通、最大位移等等。与此不同的是，基于sph的强度折减法在模拟计算过程中，可以连续求解从小应变到大变形的整个变形过程，即得到潜在滑面，而且还可以求解破坏区、破坏时间、移动距离等。

s8:根据潜在滑面分析出滑面位置和形态。

当土体承受的剪切力达到最大抗剪强度后，随着土体的连续变形，土体抗剪强度便降到残余强度，假设土体的抗剪强度在从峰值强度向残余强度过渡的过程中保持不变，不考虑滑面上岩土体材料的剪切强度的减弱，也就是说，切面残余剪切强度仍等于峰值剪切强度。因此，最大剪切应变率具有明确的物理意义，对于时间δt内发生的剪切位移是最大的，以此作为判定滑面通过的标准是合理可行的。通过连接每个单次搜索区间内最大剪切应变率sph粒子，便得到了滑面的位置和形态。

采用sph方法判断滑面位置和形态的具体方法包括：

根据潜在滑面确定滑面的搜索区域；

设置单次sph粒子搜索区间；

计算单次搜索区域内的最大剪切应变率sph粒子；

判断单次搜索区域内的最大剪切应变率sph粒子是否连续贯通；

若是，则得到滑面的位置和形态；

若否，则在下一特定间隔时间步长时继续计算单次搜索区域内的最大剪切应变率sph粒子。

本发明实施例提供的一种基于sph的土坡滑面分析判定方法，采用sph方法对土坡从小变形到大变形的连续求解确定土坡滑面的具体位置和形态，能更准确地预测在安全与危险的边界处的岩土材料的变形。

以下采用具体的实例对上述方法进行详细说明：

如图2所示，示出了模拟的计算模型，边坡由均质土体构成，为了尽量消除竖直边界(h)条件对模拟结果的影响，h0是滑面的高度，模型加大了底面(b0)的尺寸，模型底面总长度为b。模型左侧是90°垂直边界，右侧是坡角是45°的斜坡，表1列出了本次实例模拟使用的参数。模拟计算中，使用了drucker-prager屈服准则，模型底部的边界使用了无滑移边界处理，模拟中考虑了jaumann应力速率，没有考虑孔隙水压力的影响，在静止土压力的基础上对边坡进行重力加载。

表1数值算例模型参数

根据本申请的方法，通过模型尺寸与参数建模计算，在初始强度参数(c＝12.31kpa，φ＝31.30°)下，计算结果中模型随时间推移并没有发生明显的形变，说明边坡处于稳定状态。按照强度折减的理念，逐渐增加折减系数，依次进行。经过传统的瑞典条分法计算，该算例边坡的安全系数为0.855。土坡处于不稳定状态，即将在重力作用下发生失稳破坏、滑动。图3示出了采用该实例的方法得到的某一个时刻滑面的示意图。图4示出了采用该实例的方法得到的某一时刻的最大剪切应变率sph粒子位置的分布图。

目前，极限平衡法是边坡滑面确定、稳定性分析中最成熟、最方便、应用最广的方法。常用的极限平衡法有瑞典条分法、毕肖普条分法等，由于几种方法所算得的最危险滑面位置与安全系数较为接近，本申请以瑞典条分法为代表求解滑面，并与本方法进行比较分析。根据两种方法求解得到的滑面对比图如图5所示。1号线表示用瑞典条分法获得的最小安全系数的圆弧，2号线代表采用本申请的方法获得最大剪切应变率粒子滑面，可以看出：由于滑面形状符合预先的假设，瑞典条分法滑面是完美的圆弧形状，而本申请的方法求解的滑面为平滑曲线。

根据安全系数来判断剪切应变的实验数据参考表2。

表2实验数据

根据表2的数据得到：当安全系数接近1.1时，虽然没有观察到大的变形，但是剪切应变是可以观察到的。在安全系数最小的情况下，在坡冠和块体滑移处形成裂缝在实例中被证实。对于内聚摩擦材料(算例2)，当安全系数超过1.0时，一个小的剪切带可以看到并且δ/h小于2％，但是不存在大变形。反之，当安全系数小于1.0时，剪切带发育明显，边坡变形明显。

综上所述，本申请依据sph方法对土坡变形过程连续求解的方式，可以确定土坡滑面的具体位置并且评估土坡的稳定性，结果与经典瑞典条分法比较有相似的趋势。此外，常规的圆弧条分法计算只能得到安全系数，而本申请的方法不仅可以评估稳定性，还可以评估变形。在情况下，虽然瑞典条分法最小安全系数值大于1.0，达到1.1，但采用本申请的sph方法中土坡还是出现了少许变形。模拟结果表明，本申请的sph法比极限平衡法在计算中更加趋于安全。造成这种差别的原因在于，传统的极限平衡法预先设定的滑面的形状，只考虑了满足力和力矩的平衡条件，忽略了岩土体应力-应变的关系，特别是变形的发展。也没有考虑土体的强度条件，因此无法给出位移的信息，而实际上判别边坡失稳的条件中，最重要的即是发生较大位移。因此，本申请的基于sph的土坡滑面分析判定方法依据最大剪切应变率原则提出了滑面sph粒子确定方法，得到了理想弹塑性土体下土坡的滑面。为土坡滑面确定与分析提供了一种新的、有效的方式。该方法更能准确地预测在安全与危险的边界处的岩土材料的变形。为边坡稳定性问题中提取出各种有用的信息提供数据支撑。

在上述的第一实施例中，提供了一种基于sph的土坡滑面分析判定方法，与之相对应的，本申请还提供一种基于sph的土坡滑面分析判定系统，请参考图6，其为本发明第二实施例提供的一种基于sph的土坡滑面分析判定系统的结构框图。由于装置实施例基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。下述描述的装置实施例仅仅是示意性的。

如图6所示，示出了本发明另一实施例提供的一种基于sph的土坡滑面分析判定系统的结构框图，包括：粒子模型生成模块、应力应变分析模块、修正调整模块和滑面分析模块，所述粒子模型生成模块用于根据所需研究土坡的剖面信息，生成粒子模型；所述应力应变分析模块用于确定计算中需要采用的数值处理技术，并设置屈服准则和时间步长；对粒子模型进行粒子搜索及配对，计算粒子密度，对土坡分别采用弹塑性本构模型进行应力应变的求解得到粒子应力应变结果，在弹塑性本构模型中，岩土体的抗剪强度采用粘聚力和土体的内摩擦角表示；所述修正调整模块用于根据粒子应力应变结果判断岩土体的应力状态，采用修正方法对粒子的应力状态进行调整；所述滑面分析模块用于根据sph边坡稳定性分析强度折减法对岩土体抗剪强度参数进行折减，得到边坡的安全系数及潜在滑面，根据潜在滑面分析出滑面位置和形态。

滑面分析模块根据潜在滑面分析出滑面位置和形态的具体方法包括：

据潜在滑面确定滑面的搜索区域；

设置单次sph粒子搜索区间；

计算单次搜索区域内的最大剪切应变率sph粒子；

判断单次搜索区域内的最大剪切应变率sph粒子是否连续贯通；

若是，则得到滑面的位置和形态；

若否，则在下一特定间隔时间步长时继续计算单次搜索区域内的最大剪切应变率sph粒子。

以上，为本发明第二实施例提供的一种基于sph的土坡滑面分析判定系统的实施例说明。

本发明提供的一种基于sph的土坡滑面分析判定系统与上述基于sph的土坡滑面分析判定方法出于相同的发明构思，具有相同的有益效果，此处不再赘述。

如图7所示，示出了在本发明另一实施例提供的一种智能终端的结构框图，该终端包括处理器、输入设备、输出设备和存储器，所述处理器、输入设备、输出设备和存储器相互连接，所述存储器用于存储计算机程序，所述计算机程序包括程序指令，所述处理器被配置用于调用所述程序指令，执行上述实施例描述的方法。

应当理解，在本发明实施例中，所称处理器可以是中央处理单元(centralprocessingunit，cpu)，该处理器还可以是其他通用处理器、数字信号处理器(digitalsignalprocessor，dsp)、专用集成电路(applicationspecificintegratedcircuit，asic)、现成可编程门阵列(field-programmablegatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

输入设备可以包括触控板、指纹采传感器(用于采集用户的指纹信息和指纹的方向信息)、麦克风等，输出设备可以包括显示器(lcd等)、扬声器等。

该存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据。存储器的一部分还可以包括非易失性随机存取存储器。例如，存储器还可以存储设备类型的信息。

具体实现中，本发明实施例中所描述的处理器、输入设备、输出设备可执行本发明实施例提供的方法实施例所描述的实现方式，也可执行本发明实施例所描述的系统实施例的实现方式，在此不再赘述。

在本发明还提供一种计算机可读存储介质的实施例，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时使所述处理器执行上述实施例描述的方法。

所述计算机可读存储介质可以是前述实施例所述的终端的内部存储单元，例如终端的硬盘或内存。所述计算机可读存储介质也可以是所述终端的外部存储设备，例如所述终端上配备的插接式硬盘，智能存储卡(smartmediacard,smc)，安全数字(securedigital,sd)卡，闪存卡(flashcard)等。进一步地，所述计算机可读存储介质还可以既包括所述终端的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述终端所需的其他程序和数据。所述计算机可读存储介质还可以用于暂时地存储已经输出或者将要输出的数据。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的终端和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

在本申请所提供的几个实施例中，应该理解到，所揭露终端和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口、装置或单元的间接耦合或通信连接，也可以是电的，机械的或其它的形式连接。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。