基于SPH法的排土场边坡破坏过程数值模拟研究

2022-07-12魏良针王建平应涛涛龚礼岳吴则祥

科技创新与应用 2022年20期

魏良针，王建平，杨磊，应涛涛，龚礼岳*，吴则祥

（1.温州市铁路与轨道交通投资集团有限公司，浙江温州 325000；2.温州大学建筑与土木工程学院，浙江温州 325035）

边坡失稳破坏导致的灾害具有发生频率高、突发性强、经济损失大、社会影响大等特点，排土场边坡的破坏作为其中典型问题频发，使得排土场的滑坡灾害的防治一直以来都受到较多的关注，因而排土场边坡的破坏特征分析是边坡工程中的一个重要研究课题，也一直是岩土工程领域的热点研究方向之一[1]。

目前对排土场边坡稳定性的研究主要集中在边坡失稳的起始阶段[2-4]，然而滑坡的整个过程不仅仅只包括滑坡的起始，还包括滑动土体的迁移和最终的堆积。另外，值得注意的是，滑坡的最终后果往往取决于滑动土体的搬运和最终沉积。这可能是由于以往的常规方法在模拟滑坡发生后土体的大变形方面存在困难所致。然而，随着计算机技术的不断发展与硬件性能的逐渐提高，近年来发展起来的拉格朗日型无网格粒子法却能够很好地解决这类问题，特别适合模拟岩土工程中的大变形问题，例如排土场整体破坏的全过程模拟，而且具有较高的计算效率和较好的数值稳定性。

近年来，基于拉格朗日粒子的无网格方法得到了发展，并越发流行。这些方法在网格法难以建模的问题上更具优势，因为它们或多或少地避免了网格的使用。常见的方法包括：光滑粒子流体力学法（SPH）[5]，物质点法（MPM）[6]，粒子有限元法（PFEM）[7]。在这些方法中，SPH确立已久，并且在天体物理、流体动力学和固体力学中有着悠久的应用历史。近年来，SPH还被用于岩土建模领域中。已知的应用包括大型形变，颗粒流动，土地结构间相互作用，滑坡和泥石流、流体-土壤混合物动力学，开裂和地下爆炸[8-11]等一系列方向。结果表明，SPH能够很好地解决基于网格的方法所不能解决的上述问题，而排土场边坡滑坡问题的模拟也取得了较大的进步。因此，SPH是有意义的，并有潜力应用于更广泛的学术和工程岩土应用[12-14]。

在以往学者的研究中，鲜有针对排土场边坡的破坏特征影响的研究。因此，本文基于SPH（Smoothed particle hydrodynamics）方法，选了25组不同的粘聚力与摩擦角的组合研究其对均质土体边坡的滑动距离的影响，并在此基础上模拟了深圳一处排土场的滑坡现象。在边界条件复杂，破坏面位置和形式未知的条件下，本文基于SPH法模拟了排土场边坡失稳的完整过程，准确清晰地反映了土坡的破坏机理，计算得出土体内部各处的应力与应变，为此类滑坡灾害的防范与治理提供了数据支持。

1 SPH法

SPH是一种基于拉格朗日粒子法，因此，在岩土问题中，SPH控制方程写成如下拉格朗日形式：

其中，ρ表示密度，ν表示速度，σ表示柯西压力张量，g表示由外力引起的加速度，如绝大多数岩土问题中常见的重力加速度，表示材料时间导数，∇ 表示梯度算子。我们没有考虑孔隙压力的情况，所以材料可能是干燥的或者可以被总压力建模。在许多仿真模拟中，方程（1），即连续方程可以删除，因为我们应用了拉格朗日粒子，质量守恒总是满足的。连续方程只有在基于状态方程（EOS）的本构模型中使用。

在SPH中，计算域使用粒子进行划分，粒子携带场变量并随材料移动。控制方程可以通过跟踪粒子的运动和携带变量的变化来求解。要做到这一点，必须使用SPH插值离散化控制方程。场函数可以由以下积分插值近似：

其中，＜＞表示值的近似，Ω是核函数W（x-x′，h）的影响域，在本文之后简写为W。W取决于距离||x-x′||和参数h，即平滑步长。本文的球形支持域W半径为2h。此外，核函数W必须满足一定的条件，如归一化条件，增量函数性质和紧凑支持条件。常用的核函数包括高斯核函数（Gaussian kernel），三次样条核函数（cubic spine kernel）和温德兰核函数（Wendland kernel）。在本文中，由于可以有效阻止粒子聚集，我们只使用C2温德兰核函数。

类似的，场函数的空间导数可以通过以下方程近似：

其中，∇x′表示导数在x′处取值。上述方程仅在支撑域与计算域边界不相交时有效。然而，这个条件在边界附近的区域通常是不满足的。因此，SPH需要一些边界的特殊处理来解决所谓的边界缺陷。

由于可以在支撑域中累加所有粒子的作用，SPH中方程（3）（4）的连续积分插值形式可以重写为以下形式：

其中，方便起见我们省略了＜＞符号；xi是场函数估计的粒子，xj是支撑域中的粒子；∇iWij是 ∇xiW（xi-xj，h）的简写形式。由于核函数的对称性，∇iWij=- ∇jWij；mj表示粒子j的质量，表示粒子的体积。

考虑SPH插值方程（5）（6）的离散形式，控制方程可以重写为以下形式：

尽管其他文献中也有控制方程的不同离散形式，上述方程能够保留线性和角动量，并被广泛使用。

2 模型构建与边坡破坏特征量化方法

本文基于C++编写了SPH计算模拟程序，采用编写命令流的方式来构建边坡模型，模型尺寸如图1所示。所对应的SPH粒子模型如图2所示。

图1 边坡尺寸示意图

图2 基于SPH粒子的边坡模型

对于如图2所示的二维SPH边坡，其最靠底面与侧面的粒子为受约束粒子，这些粒子的x、y两个方向的位移都始终为0。

为了探究土体参数对边坡破坏特征的影响规律，采用Run-out Distance（RD）来定量表征边坡的破坏特征。图3是通过SPH模拟得到的边坡土体破坏后SPH粒子的位移云图，图3用不同的颜色表示不同程度的土体的位移大小，可根据位移大小变化识别失稳滑动带。

图3 边坡破坏的SPH粒子位移云图

确定了滑动粒子后，在水平方向上用Run-out Distance定量描述边坡的破坏特征，Run-out Distance是指从边坡破坏前的坡脚到沉积稳定后的边界点的距离。

3 粘聚力与摩擦角对边坡破坏特征的影响分析

本节主要针对粘聚力和摩擦角对边坡破坏特征的影响展开分析。基于第一节介绍的SPH法与第二节所构建的边坡模型，给边坡赋值不同的粘聚力和摩擦角进行分析。共对25组不同粘聚力和摩擦角的组合进行了分析，每组粘聚力和摩擦角的取值见表1。其余土体参数的取值均相同，即重度为γ=20 kN/m3，弹性模量为E=20 MPa，泊松比为μ=0。

表1 25组边坡的粘聚力和摩擦角取值

�

将25组中的9组不同粘聚力和摩擦角边坡的破坏特征分析示例结果列于图4。图4中No.为对应的边坡编号，RD分别为计算出的边坡破坏特征参数，其定义见第二节。从图4中可以看出，粘聚力和摩擦角的变化对边坡破坏特征的影响很大。对于编号N=1（即c=5 kPa，φ=10°）的边坡，边坡的破坏程度十分明显，其破坏特征参数的数值都很大。当粘聚力和摩擦角均增加到最大值（即c=25 kPa，φ=30°）时，边坡的破坏特征参数趋于0，说明具有该土体参数的边坡处于稳定状态，其变形很小。

图4 9组不同粘聚力与摩擦角边坡的SPH模拟结果

接下来本文分别探究了粘聚力、摩擦角与边坡的各破坏特征参数RD之间的定量关系。为了便于观察粘聚力和摩擦角对边坡破坏参数的影响，本文将粘聚力、摩擦角与对应的边坡破坏特征参数绘制在同一个三维图上，如图5所示。可以看到，随着粘聚力和摩擦角的增大，边坡的破坏特征参数RD随之减小。本文采用线性平面对这些散点进行拟合，得到边坡粘聚力c、摩擦角φ与边坡破坏特征参数的关系式，结果如公式（9）所示。

图5 粘聚力c、摩擦角φ与特征参数RD的关系

研究结果表明：非线性平面可以很好地用于拟合粘聚力、摩擦角和边坡破坏特征参数的关系，拟合优度R2都高达99%。研究结果表明，对于均质土体边坡，在本文的计算条件下，单个因素（粘聚力或摩擦角）对边坡破坏特征参数的影响也是非线性的。该研究结果对于工程实际中的计算和预测边坡的破坏特征、指导设计、施工和后期加固具均有重要的意义。

4 工程实例

本文选取了广东省深圳市某排土场作为工程实例进行滑坡模拟。图6展示了该排土场滑坡前后地形对比图，该滑坡现场可分为3个区域，由高到低分别为：倒梯形的滑坡物源区、脆弱的滑口区和平缓的堆积区。滑坡物源区原为采石场，工程中期作为排土场倾倒渣土，下方堆积区原为平坦的工业园区，现受到严重破坏。缺口出现在该排土场北面，此处为一缓坡，坡高约30 m，滑动土体即余泥渣土从该缺口处涌出，下滑至工业园区，以较高的速度破坏地面建筑物，冲击原地表土体，并最终形成大面积的堆积体。

图6 滑坡前后对比图

以该排土场边坡作为研究对象，根据勘探资料和该排土场剖面图，基于SPH法建立动力显式模型。该排土场模型以平行于边坡走向为X轴，指向排土场内部为正；垂直边坡走向为Y轴；竖直方向为Z轴，以向上为正。排土场设计余泥渣土堆积方案为：分10级堆放，每一级高度10 m，坡度1∶2.5，相邻两级间设宽3 m的马道，滑口处于第1级顶面，滑口以下土体未发生滑动，因此可处理为基岩，滑口以上部分全为余泥渣土，高度为90 m，等效坡角为20°。取纵剖面进行分析，其几何模型如图7所示。

图7 边坡尺寸示意图（单位：m）

在有限元分析软件ABAQUS中建立对应的SPH三维边坡粒子模型，如图8所示。边界的大小直接影响应力与应变的分布，本模型的边界条件为：四周和底部为位移限定边界，坡面为自由边界，只允许竖向沉降。网格疏密对计算精度有着明显影响，网格划分太稀，则计算误差大，太密则需要耗费过多的效能，为兼顾数值模拟的效率和计算精度，本文中单元尺寸设为5 m，划分网格后共49 619个节点，62 318个单元。