基于无网格方法的边坡支挡结构变形特征研究

徐晓东

（内蒙古包钢钢联股份有限公司巴润矿业分公司）

随着露天矿开采规模的扩大及开采程度的深入，采矿边坡滑坡俨然已成为一种主要的地质灾害［1］，对矿产开采、采矿设备和人员安全造成严重影响。目前的研究对边坡稳定性和失稳原因的探究较为成熟，但对于边坡失稳后所造成的冲击影响研究较少［2］。因此，明确边坡失稳后的滑动位移和冲击效应是采矿工程领域重要的研究课题。

边坡失稳时通常会出现大变形特征，当发生流动大变形时，变形不再符合小变形理论，而是进入了非线性的大变形状态［3-4］。此时，传统的数值模拟方法，如有限元法（FEM）和有限差分法（FDM），由于其网格属性，容易在计算大变形时产生网格扭曲，导致对边坡失稳后的滑移大变形分析存在一定困难［5］。为了克服目前常用的数值模拟方法中的不足，离散元法（DEM）、非连续变形分析法（DDA）、无网格伽辽金法（EFG）以及细胞自动机法（CA）等离散介质力学方法和无网格方法得到了广泛应用。然而，这些方法在研究边坡变形方面还存在许多缺点，如难以确定非物理参数和直接描述土体的应力-应变关系等［6-10］。

针对滑体冲击影响的问题，国内外也进行了一定的研究。在理论方面，通常采用速滑计算方法计算滑体的冲击能，如能量法［11］和条分法［12］等。根据相似性原理，利用物理模型实验模拟滑体的冲击过程，评估冲击的影响也有涉及［13-15］。然而，岩土体是1 种散体，理论方法在经过过多的简化和假设后，不能准确描述滑体滑动过程及冲击破坏过程［16］。基于物理模型的相似性设计相对复杂麻烦，其制作周期长、费用高［17］。

基于此，本文采用1种无网格方法—光滑粒子流体动力学法（SPH）对边坡变形特征及冲击效应进行数值模拟研究。结合HBP（Herschel-Bulkley-Papanastasiou）本构模型，模拟了经典的Poiseuille 流问题，同时对比分析了Bingham 模型的流体运动问题。进一步建立露天矿边坡数值模型，对边坡失稳演化过程及滑体对支护结构的受力变形规律以及支护效果进行模拟分析，分析柔性支护结构在边坡失稳过程中的冲击作用及受力特征，为后期同类工程的施工提供一定的借鉴。

1 SPH基本理论

1.1 基本控制方程

SPH方法是一种具有拉格朗日特征的粒子法，该方法的基本思想是对研究材料进行离散化处理，将其离散为能够代表研究材料本身的粒子（图1）。每个粒子携带自身的材料属性，如质量、密度、速度和能量等，并通过跟踪这些粒子特征反映研究材料的状态，同时遵循物质的运动守恒控制方程。

图1 中h为光滑长度；k值为影响系数，表示影响域大小，一般根据所研究问题来确定k值，随着k值的增大，计算精度会提高，但同时也会降低计算效率；W为核函数；i为计算粒子；j表示相邻粒子；rij表示粒子i到j的距离。

SPH计算方法中，研究具有材料强度的控制方程有粒子密度、速度和位置变化方程，如下：

式中，ρi为第i个粒子密度，g/m3；xα为空间坐标；vα为速度矢量；σαβ总应力张量，Pa，其拉应力为正，压应力为负，一般由各向同性压力P和剪应力τ构成，计算公式为σαβ= -pδαβ+ταβ；α，β为坐标方向。

经过离散化处理，式（1）和式（2）在SPH 方法下的计算公式：

式中，mj为第j个粒子的质量，g；P为正应力，Pa；g为重力，N；Wij为核函数，通常情况下，3次B样条函数在计算精度和计算效率均比较好。

1.2 HBP本构模型

边坡从失稳到流滑的动力过程十分复杂，一般认为其经历了失稳—滑动—流滑3 个过程。为准确地描述边坡在大应变状态下剪应力-剪应变非线性关系，本文将变形体视为具有可变黏度的黏性材料，采用HBP 模型表征变形体的流变特性。HBP 模型由Papanastasiou［18］提出，剪切力是通过与应变张量相关的剪应力张量表示，其本构方程如下：

式中，τij为剪应力，Pa；μeff为有效黏滞系数，描述流体的流变特性；εαβ为应变张量，Pa。

其中，有效黏滞系数计算公式：

式中，γ为剪切应变率；μ为表观动态黏度系数；τy表示屈服应力，Pa；m，n为计算常数，可以模拟剪切变薄或剪切增稠行为，当m=0 和n=1 时，模型变为牛顿模型，而m→∞且n=1时，模型简化为Bingham模型。

因此，最终SPH动量方程可表示如下：

1.3 时间步长确定

在SPH 求解过程中，需要对控制方程和本构方程进行离散求解，求解方法有蛙跳算法、预测矫正法以及龙格-库塔法等，在求解过程中一个关键参数是时间步长的确定。一般情况下，蛙跳算法所需存储量低，在每一次计算中只需要进行一次优化估值等优点，其迭代过程如下：

式中，ci代表声速，m/s；vi代表粒子速度，m/s；CFL是稳定性常数，小于1.0。

2 HPB本构模型数值验证

2.1 Poiseuille流动

选取典型的Poiseuille 流为模拟算例，流体模型为在2个固定不动的平板之间充满水，并且初始状态下水处于静止状态。接着，当受到水平方向的外力作用时，水开始流动。采用SPH 方法模拟瞬态流动过程时，建立如图2 所示的Poiseuille 流模型，计算区域为2 m×1 m 的矩形。在此计算中，主要的计算参数：粒子个数1 640 个，粒子的初始间距r0=0.025 m，光滑长度h=1.2r0，水平外力F=1 N，初始参考密度ρ0=1 000 g/m3，时间步长为1.0 s。

该流体模型在初始状态下（t=0）位于2 块固定的无限大平板之间且处于静止状态。然后，在t＞0 时，一个恒定的水平作用力作用于整个流场，导致最初静止的流体开始逐渐流动。流动速度如图3所示。

从图3中可以看出，流体的速度沿着水平轴中心线具有对称性。即靠近平板附近的流体由于黏附在静止的平板上，因此它们的速度较小且近似为零；而流道中心的速度最大，这与日常物理现象相一致。在经历0.8 s 后，流动达到最大状态，其最大流速约为0.006 m/s，如图3（c）所示。为了进一步验证方法的正确性，在流场中选取了20个空间测量点（图2）采用有限容积法（FVM 方法）与所提方法进行对比，得到HBP 本构模型条件下Poiseuille 流的速度对比图，如图4所示。

由图4 可看出，在时间步长分别取1.0，2.0 和5.0 s 时，平行板两侧粒子速度几乎相同，靠近两侧粒子速度逐渐减小，各时刻流速分布图趋势均是抛物线型，并且左右对称。SPH 方法与FVM 方法的所得的结果具有较好的吻合性，验证了SPH 方法在模拟黏弹性流体瞬态流动问题的准确性和有效性。

2.2 Bingham模型对比分析

基于HBP 和Bingham 2 种本构模型，建立三维水槽流动模型，模型两侧为Bingham 本构模型，中间为HBP 本构模型。粒子初始间距0.008 5 m，光滑长度为粒子初始间距的1.2 倍，初始时间步长为1.0×10-4s，对比不同本构关系下的SPH数值模拟结果，得到了各瞬态下粒子的分布情况以及速度场分布情况，如图5所示。如图5（a）所示，水流运动过程中最大传播速度出现于HBP 本构模型的流体前端部位，并逐渐向后端递减，且相对于Bingham 流体模型更为明显。随着时间的推移，流体前端速度继续加大，并且运动距离明显大于上下两侧流体，如图5（b）所示。如图5（c）所示，t=60.0 s 时，流体速度达到峰值后逐渐减小至零，达到静止平衡状态，从横向速度场分布及运动距离上来看，HBP 本构模型可以更好地模拟变形流变特征，最接近水体自由运动。

3 支挡变形数值模拟

3.1 一维流体与支挡

考虑自由流体与支护结构相互作用模型，其上部粒子数为12 160 个，表征自由流体，受自重作用，高度4 m，宽度2 m，如图6（a）所示。下部粒子数352模拟固体支挡结构，对其进行封闭。如图6（b）所示，底部支挡结构收到压力主要集中在两侧，支挡结构未产生变形。随着时间增加，如图6（c）所示，底部支挡结构应力集中转移至支挡结构中心位置，且在压应力的作用下，支挡结构发生弯曲，但是未达到损伤、断裂现象。为了更清楚地分析流体与支护结构之间的相互作用，选取验证模型底部中点为监测点，其位移、压力与时间关系如图7所示。

如图7（a）所示，验证模型测点A 位置位移随时间增加而增加，当计算时间达到0.4 s 时，监测点位移大约为0.16 m。如图7（b）所示，应力也具有同样的规律，即随着时间的增加，支挡结构所承受的压应力也随之增大。通过分析该模型，可知上述方法在变形和支挡结构方面的可行性和正确性。

3.2 柔性支挡

在模型初始高度54 m，初始长度180 m，支护结构长度14 m 的条件下，边坡变形失稳运动演化过程及支护冲击相互作用如图8 所示。数值模拟得到了各瞬态下粒子的分布情况以及冲击效应分布情况，如图9所示。

由图8（a）可以看出，当边坡失稳后，滑体与支护结构接触时，支护结构受到的最大冲击力主要作用在土体与接触部位，即支护结构的底部，但是由于滑动速度较小，冲击不足以使支护结构发生变形。随着边坡滑动时间的增加，边坡加速滑动向前，如图8（b）所示，支护结构阻挡了绝大数失稳土体，且受到的冲击力逐渐上移，支护结构发生弯曲倾倒变形。当边坡失稳速度达到最大时（图9（a）），此时滑体体积也是最大，部分失稳土体越过支护结构，此时支护结构变形达到最大值，如图8（c）所示。随着边坡失稳速度逐渐减小直至边坡停止，如图9（b）所示，此时达到新的平衡状态，支护结构所受到的冲击力突然骤降，柔性支护结构恢复自身变形，但是其底部出现较小的残余应力作用（图8（d）），主要是由于失稳土体的土压力造成的。

为了能够进一步证明该方法的可行性，基于上述模型，对支护结构不同锚固深度时，边坡失稳冲击效应进行了研究，得到边坡在滑动过程中，支护结构在不同锚固深度的压力和位移变化数据，如图10 所示。

通过图10 可以看出，随着支护结构锚固深度的增加，边坡在滑动过程中支护结构受到冲击作用力呈现增长的趋势，并且最终由于边坡停止失稳，冲击力均趋于非零稳定值；当锚固深度3.0 m时，支护结构竖向位移变化平稳，表明边坡滑动过程中，支护结构锚固深度越大，支护结构在受到冲击作用时结构越稳定。

4 结 论

本文基于光滑流体粒子动力方法，结合HBP 本构模型，构建了模拟边坡大变形与支挡结构相互作用过程的数值模型，得出以下结论：

（1）选取HBP 模型作为边坡变形本构模型，有效地克服了土体变形在塑性应变过渡段及大变形情况下可能呈现出非线性流变特性的问题，其对于分析变形较大和两相物质相互作用的问题具有一定的可行性。

（2）由于SPH 方法不依赖于网格及其纯拉格朗日特性，自然而然在处理具有变形或者不连续问题上有比较好的优势，为研究该类问题提供了一种新的方法和研究思路。

（3）露天矿边坡在地震荷载作用下的稳定性也是工程设计关注的重点问题，后续工作将基于SPH深入开展边坡动力失稳机理和滑动过程研究。