基于光滑粒子流体力学法分析坡高和坡角对边坡破坏特征的影响

冉光炯* 刘 勇 白 皓

（四川高速公路建设开发集团有限公司，四川 成都610000）

在我国地形陡峭的山区，由于地质环境极为特殊，地质作用非常复杂，极易导致边坡失稳破坏[1]。

在已有的文献中，研究人员常采用FEM（有限元方法）[2]、DEM（离散元法）[3]和FDM（有限差分法）[4]来研究失稳破坏后边坡的变形特征。然而，在模拟边坡失稳破坏时，有限元法和有限差分法都具有一定的缺点，如网格变形过大土体破坏引起的变形。DEM法虽然有利于模拟大变形问题，但DEM在处理工程尺度的模拟问题时，其计算效率还有待提高，且DEM的微观参数难以标定。

近年来，无网格平滑粒子流体力学（SPH）[5-8]已应用于探究岩土材料的大变形问题。已有研究表明，该方法可以有效地预测边坡破坏引起的土体大变形。本文采用SPH 方法对该问题进行了研究，采用弹塑性Drucker-Prager 模型[9]来模拟土体，探究了不同坡高和坡角的边坡的破坏特征。

1 光滑粒子流体力学法

在SPH 方法中，使用一组运动粒子来模拟连续介质的运动；每个都分配有恒定的质量，并在相应位置“携带”场变量。通过加权求和，将连续场从各种粒子中插值出来，其中权重根据假定的核函数随距离而下降。然后，“粒子间相互作用”的范围对应于内核的“支持域”的半径。类似地，空间导数是通过在支撑域上的插值来计算的。SPH 插值理论可参考Liu & Liu[10]和Monaghan[11]的论文。

以下简要介绍土体运动方程的SPH 离散化。连续体的运动可以用以下方程来描述：

在SPH 框架中，此等式通常离散为

其中i 表示考虑中的粒子；ρi和ρj分别为粒子i 和j 的密度；N 是“相邻粒子”的数目，即在粒子i 的支撑域中的粒子的数目；mj是粒子j 的质量；Cαβij是一个稳定项，用于消除应力波动和拉伸不稳定[12]；Wij是核函数，这里选择的是三次样条函数[13]。公式(4)中的稳定项由两个分量组成：人工粘度和人工应力，其计算方法与文献[12]相似。但人工粘性项的声速在此由下式计算

式中，E 是土的杨氏模量。

公式(4)在有效应力张量已知的情况下，可以使用标准的蛙跳算法[14]进行积分。

本构模型描述了给定材料的应力和应变之间的关系。根据经典塑性理论，将弹塑性材料的总应变率张量ε˙分解为两部分：弹性应变率张量ε˙e和塑性应变率张量ε˙p：

式中，S˙'αβ为偏有效剪应力张量；G 为剪切模量，σ˙'m为平均有效应力。

塑性应变率张量ε˙αβp按塑性流动法则计算：

其中λ˙是塑性比例系数的变化率，gp是塑性势函数。本文采用非关联塑性流动规律的Drucker-Prager 模型[9]，假设屈服面固定在应力空间中，只有当应力状态达到屈服面时才会发生塑性变形。因此，只有在满足下列屈服准则的情况下，才会发生塑性变形：

其中I1和J2是应力张量的第一和第二不变量；αø和kc是Drucker-Prager 常数，这是根据库仑材料常数c(内聚力)和φ(内摩擦角)计算得出的。

上述土体本构模型需要6 个土体参数：粘聚力系数、摩擦角、剪胀角、杨氏模量E、泊松比和土体密度ρ。Bui 等人[12]在SPH 框架内对当前土体本构模型进行了验证。其中SPH 关于土体材料重力流动和粘性土上的地基问题的计算结果与实验，和有限元法一致。

2 模型构建与计算流程

本文基于开源的Dual-SPH 软件[15]，采用编写命令流的方式来进行不同坡高和坡脚的边坡模型的构建。其中，坡高和坡脚的定义如图1 所示，采用Dual-SPH 构建出来的边坡模型如图2所示。注意在建立模型的时候，我们将模型沿垂直直面方向拉伸了一个单位长度(1mm)的宽度，使其形成一个假三维模型。

图1 边坡坡高y(m)与坡脚θ(°)的示意图

图2 Dual-SPH 构建的边坡模型

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

图3 Dual-SPH 边坡模型中粒子的位移图

确定了滑动粒子后，在水平方向上用滑动距离S 来量化评价边坡失稳后的破坏特征，滑动距离S 是指从边坡破坏前的坡脚到沉积稳定后的边界点的距离。

为了控制变量，在所有边坡模型中，其土体参数的取值都一致，见表2。

表2 计算中用到的土体参数取值

表3 25 个边坡的坡脚和坡高及其滑动距离的计算结果

3 坡高与坡脚对边坡稳定性的影响分析

本小节在第二小节的基础上，主要针对边坡的坡脚和坡高对其稳定性的影响展开分析。首先，采用第2 小节所述的方法，在中分别构建了25 个坡脚和坡高不同的边坡，具体每个边坡的坡高和坡脚如表3 所示。然后采用Dual-SPH 分别计算每个边坡的滑动距离，其结果表3 所示。

为了便于直观地观察边坡高度与边坡坡度对滑动距离的影响，我们将边坡高度、边坡坡度与滑动距离绘制在同一个三维图上，如图4 所示。可以看到，随着边坡坡度的增大和边坡高度的增高，边坡的滑动距离随之显著增大。为了研究三者之间的定量关系，我们用二次曲面对这些散点进行拟合，得到边坡坡度、边坡高度与滑动距离的关系式如公式（17）所示，其中x 表示边坡坡度，y 表示边坡高度，S 表示边坡的滑动距离。从图中可以看出，该公式的拟合优度R2为99.7%。这表明：在本文的分析框架与参数条件内，边坡高度、边坡坡度与滑动距离的关系可以很好地用二次曲面来进行拟合。图5 分别展示了不同边坡高度与边坡坡度的SPH 粒子位移云图。

图4 边坡高度、边坡坡度与滑动距离的关系

图5 不同边坡高度与边坡坡度的SPH 边坡的位移云图

4 结论

本文基于开源的Dual-SPH 软件来对边坡坡高和坡脚对边坡失稳后滑动距离的影响进行研究，分别构建了25 个不同坡高和坡脚的边坡，并对其失稳后的滑动距离进行了计算。计算结果表明：随着边坡坡度的增大和边坡高度的增高，边坡的滑动距离随之明显增大。进一步分析发现，在本文所建立的模型与参数条件下，边坡坡度、坡高与滑动距离的关系可以用二次曲面很好的拟合，且拟合优度R2为99.7%。本文的研究结果为定量分析边坡坡高、坡脚与边坡失稳后滑动距离的关系奠定了基础。后续将继续针对边坡参数对其失稳后破坏特征参数的影响展开定量研究。