滑坡-碎屑流冲击导引结构的离散元模拟

张睿骁，樊晓一,2，姜元俊，杨海龙

(1.西南科技大学土木工程与建筑学院，四川 绵阳 621010；2.工程材料与结构冲击振动四川省重点实验室，四川 绵阳 621010；3.中国科学院成都山地灾害与环境研究所，四川 成都 610041；4.成都理工大学地质灾害防治与地质环境保护国家重点实验室， 四川 成都 610059)

由于地质构造复杂，气候多变，我国多发滑坡、泥石流之类的地质灾害。在山区，大部分可利用的土地是坡体边缘相对平缓的部分，如果坡体分布有不稳定的斜坡，发生滑坡后坡体前缘的构筑物区域正是堆积区，易造成严重的损失和伤亡事件。如：2015年12月20日，广东省深圳市发生一起渣土堆填物特大滑坡事故，造成33栋建筑物被掩埋或不同程度受损，73人遇难，4人失联，覆盖面积约38.5×104m2，填筑体积约627.24×104m3[1]；2016年7月1日，贵州省毕节市偏坡村金星组发生山体滑坡，掩埋金星组居民11户30人，共造成23人死亡，7人受伤[2]。

为了应对碎屑流的冲击灾害，可以通过建设拦挡结构来减小滑坡致灾范围、减弱致灾强度，为此国内外学者开展了滑坡碎屑流运动能量和冲击机理的研究。运动能量方面，有学者根据相似原理设计制作滑体下滑及冲击模型，研究滑坡与受灾体之间的作用[3]，根据冲量定律和能量守恒定律换算得到滑体冲击能[4]，并对受灾体易损性定量评估[5-6]。同时考虑滑坡体运动过程中耗散规律[3,7]，依据离散元能量转化和能量守恒计算法则[8]，建立用于评估滑速、滑距及冲击能滑体运动模型[7]。冲击机理方面，有学者采用模型试验与数值模拟相结合的方法，考虑颗粒粒径大小的影响[9]，研究不同位置[10]和高度[11]的挡墙对滑坡冲击的影响以及易损性和动力响应[12-13]。以冲击力为研究对象，分解为碎屑流拖曳力、重力和摩擦力以及被动土压力[14-16]，基于深度平均冲击理论，构建了力学模型[17-18]。然而冲击力分布的复杂性和非线性特征，与碎屑流冲击挡墙的土拱效应关系密切[19]，通过分析微观结构的异质性在流动的冲击和沉积阶段，得出土拱效应对颗粒的影响[20]，可以进一步探讨土拱效应的形成机理、冲击力拱的形成特征以及对冲击力分布变化[19]。

结合分析已有的研究成果，发现以往的碎屑流冲击灾害研究主要以运动特征和冲击机理为切入点，侧重碎屑流正面与刚性或柔性防护结构[9,21]和桥墩[22]等承灾体之间的作用效应，对碎屑流的冲击考虑偏转效应的研究较少。从海岸防护堤的导流结构使得海水返回到大海得到启示，导流结构的研究应用到了碎屑流当中[23]，但并未考虑碎屑流侧面与结构之间的作用改变其运动路径[24]，使得这些研究不够全面。

为了丰富碎屑流侧向冲击作用的研究，更加清晰、全面地说明碎屑流与不同导引结构之间的冲击效应和堆积特性，本文利用三维离散元软件，结合室内模型试验数据，分析了不同导引结构对碎屑流颗粒运动的导向作用，并为滑坡灾害防治提供建议。

1 三维离散元模型

1.1 模型建立

本模型由斜坡、侧板、底板、拦挡结构和碎屑物源组成(图1)。其中，碎屑物源距离底面垂直高度固定为1.25 m，滑槽宽度为0.5 m，滑槽高度为0.5 m，滑坡体纵向宽度2.5 m，斜坡坡度为45°。导引结构设置在坡脚处。堆积区纵向宽2.5 m，横向宽1.5 m。

图1 滑槽模型示意图(单位：mm)Fig.1 Schematic diagram showing the chute model diagram (mm)

本文所用的导引结构，长(l)、宽(d)、高(h)都相同，为250 mm，考虑了三种典型结构，凹型圆弧结构(B1)、直线型结构(B1)、凸型圆弧结构(B3)(图2)。

图2 导引结构示意图Fig.2 Schematic diagram showing the guidance structure diagram

已有学者对颗粒单元进行校准之后发现，使用多个球体黏结作为基础颗粒与单个球形颗粒相比，前者能够反映实际颗粒含有棱角的特点，数值模拟结果与试验结果更相符[20,25]。本文通过将四个球形颗粒黏结作为基础颗粒(图3)，模型中最小颗粒单元为刚性体，颗粒密度2 100 kg/m3，总质量52.5 kg。

图3 单个颗粒大样[11]Fig.3 Single particle sample[11]

1.2 数值模型参数

根据滑槽模型试验具体数据设置相应仿真模型参数。离散元数值模拟中，泊松比、剪切模量和密度为材料的本征参数，表示各种材料的自身属性，与外界条件无关。改变滑体粒径或级配，不会影响材料属性。基于模型试验实测及计算结果，确定颗粒与几何体的材料属性见表1。

表1 材料参数

图4 颗粒级配曲线[11]Fig.4 Particle gradation curve[11]

图5 碎屑流离散元颗粒模型[11]Fig.5 Discrete element model for debris flow[11]

休止角是颗粒间流动特性和摩擦特性的宏观描述，是颗粒间静摩擦系数和滚动摩擦系数标定常用的指标[27]。因此，以颗粒间的静摩擦系数、滚动摩擦系数作为试验因素，采用堆积法进行休止角试验。根据休止角试验结果分别拟合得到适用于堆积法的颗粒间静摩擦系数、滚动摩擦系数与颗粒休止角正切值的回归方程。

利用有机玻璃建立边长为0.3 m的立方体模型，模型填充0.008 m3碎屑流颗粒，在刚性底板上开展室内堆积法试验(图6)，重复进行3次，并采集堆积形态的图像，获取现场室内试验休止角的正切值(表2)。

图6 现场休止角试验Fig.6 Rest angle field test

试验次数123平均值休止角正切值0.5660.5630.5710.567

根据实测得到的岩土体的休止角正切值来求解颗粒间静摩擦系数和滚动摩擦系数。最后通过引入标定的参数进行休止角仿真试验(图7)，以休止角基本吻合为条件，完成参数标定，最终接触参数见表3。

图7 模拟休止角试验Fig.7 Simulated angle field test

物理参数颗粒-颗粒钢材-颗粒有机玻璃-颗粒静摩擦系数0.570.500.01滚动摩擦系数0.100.100.01碰撞恢复系数0.200.280.25

2 冲击力分析

通过碎屑流与结构之间的冲击力分析，可以对结构优化提出有益的建议。法向力是颗粒与接触面垂直的作用力，法向力越小，对结构的损毁越小；切向力是指颗粒作用在与接触面平行的力，切向力越大，对颗粒的导引作用越大。

如图8(a)所示，导引结构B1所受法向力随时间的变化可以分为四个阶段：线性增加、恒力阶段、线性减小和恒力阶段。从碎屑流与结构作用开始(t=0.65 s)，法向力随时间增加而增加，持续至0.9 s；在0.9～1.1 s之间法向力大小基本保持不变，1.1～1.3 s法向力逐渐减小，最后对结构作用恒定的静压力。导引结构B2和B3所受法向力的趋势相同：随时间的增加先增加后减小，最后保持不变。两种结构法向力的峰值出现时间有所不同，导引结构B3峰值法向力的时间(1.02 s)早于导引结构B2(1.09 s)。其中，导引结构B1所受的法向力最小，其峰值约为(50 N)，而导引结构B2 、B3法向力峰值分别为180 N、250 N。最后由于碎屑流颗粒堆积在导引结构后方，对结构作用的静压力使得法向力曲线最后没有归零。按照静压力的大小排序，分别为导引结构B3>导引结构B2>导引结构B1。对比三种结构法向力的曲线，导引结构B2和B3的法向力曲线波动比较大，而导引结构B1的法向力曲线波动比较小。

图8 冲击力时程曲线变化Fig.8 History curve of impact force vs. time(a)法向力；(b)切向力

而从切向力时程曲线(图8b)发现，切向力时程曲线可以分为三个阶段：线性增加、线性减小和恒力阶段。碎屑流颗粒开始与导引结构接触，切向力逐渐增大，三种结构分别达到切向力峰值。按照切向力峰值出现的时间先后进行排序，分别是结构B2、结构B3、结构B1。其中，导引结构B2和B3切向力峰值大小相同，约为70 N，而导引结构B1的切向力峰值约为160 N，是三种结构切向力峰值最大的。当碎屑流颗粒静止时，颗粒对导引结构的作用力使得切向力维持恒定的数值，结构B1、B2、B3的切向力分别为40 N、11 N、13 N。

出现以上现象的原因在于：拦挡结构在受到碎屑流冲击后，拦挡结构上的受力随着时间增加逐渐增加，导引结构后方堆积颗粒的体积也随之增加，对后缘下滑的颗粒的作用力越大。当法向力和切向力达到峰值后，拦挡结构前堆积着一定体积的颗粒，而后缘的颗粒不断下滑，与堆积的颗粒碰撞、摩擦，间接作用在拦挡结构上，法向力和切向力不断减小，最终碎屑流颗粒堆积完成。堆积在拦挡结构前方的颗粒的静压力最终使得法向力和切向力维持在一定数值。对比三种导引结构，结构B1可以有效地将碎屑流颗粒的作用力转换，结构所受的法向力最小，结构的损毁最小；切向力最大，对颗粒的导引作用最大。

3 能量分析

导引结构的设置可以一定程度上改变碎屑流颗粒的运动方向。为进一步研究导引结构的设置对碎屑流颗粒运动能量的影响，根据颗粒的速度、质量和高度，通过式(1)、(2)可以计算出每个颗粒的动能和势能：

(1)

(2)

式中：mi——颗粒i的质量；

vi——颗粒i的速度；

g——重力加速度；

hi——颗粒i与水槽基部的垂直距离。

如图9(a)所示，碎屑流颗粒的动能随时间的增加先逐渐增加，在0.8 s左右达到峰值，随后逐渐减小；无导引结构时，碎屑流颗粒的动能随时间变化趋势与有导引结构的相同，先增加后减小。从图9(b)可以发现，导引结构的存在与否，对结构的势能变化影响不显著，势能随时间增加不断减小，而最终有导引结构设置时，颗粒堆积在坡脚处的颗粒越多，堆积高度会增加，导致势能的终值大于无导引结构时。

当t=0.65 s时，碎屑流颗粒经过加速区加速后，势能与动能之间进行了部分转换，此时前缘颗粒与导引结构发生碰撞，而大部分颗粒仍在加速区滑槽上，具有较大的势能，由图9中可得颗粒的动能EK=117 J，势能Ep=297 J。对比有无导引结构时颗粒的动能大小，可以发现，导引结构的设置对颗粒动能有一定程度的耗散，减小了约20 J。经过三种不同导引结构后，颗粒的动能时程曲线基本重合，说明三个导引结构对颗粒动能耗散无显著差异。

图9 碎屑流颗粒能量时程曲线变化Fig.9 Variation in energy time history curve of debris flow particles (a)碎屑流颗粒的动能；(b)碎屑流颗粒的势能

碎屑流颗粒动能的变化，归因于颗粒与颗粒(图10a)、颗粒与滑槽(图10b)、颗粒与导引结构之间(图10c)的碰撞和摩擦。如图10(a)所示，颗粒之间的接触数量随时间增加而不断增加。在t=1.25 s之前，导引结构的设置，对颗粒之间的碰撞影响并不显著，碎屑流颗粒之间接触数量的所有时程曲线基本重合，颗粒之间碰撞和摩擦所消耗的能量基本相同；1.25 s之后，颗粒继续运动，此时仍有部分颗粒在加速区滑槽上，后缘颗粒推挤前缘颗粒向前运动。从图9(a)中可以看到：有导引结构时，颗粒的动能明显小于无导引结构时，说明无导引结构时后缘颗粒对前缘颗粒推挤作用更显著，颗粒运动得越分散；而有导引结构时，后缘颗粒的运动堆积得更紧密，颗粒之间的碰撞数量大于无导引结构时，颗粒之间碰撞和摩擦消耗的能量更大，动能减小得更快，运动时间更短。从图9(a)中可以看出，导引结构使碎屑流颗粒的运动时间减小了0.1 s。从图10(b)中可以看出碎屑流颗粒与滑槽之间的接触数量存在先增加后减小再增加的趋势；导引结构的设置大大增加了碎屑流颗粒与滑槽之间的接触数量。由于导引结构的设置，碎屑流颗粒的运动路径被阻挡，前缘的颗粒大部分堆积在坡脚处，而后缘颗粒受到前缘堆积颗粒的作用，在滑槽上与滑槽之间的碰撞时间增加。而从图10(c)可以发现，导引结构的设置，颗粒与结构之间的接触数量随着时间的增加不断增加，说明堆积在导引结构后方的颗粒越多。虽然颗粒与结构之间的接触逐渐增加，但颗粒与结构之间的接触数量明显小于颗粒之间和颗粒与滑槽的数量，说明颗粒与结构之间的碰撞和摩擦对能量的影响很小。因此，对碎屑流颗粒动能变化的影响，主要是由于颗粒与滑槽之间的碰撞和摩擦，这一发现与Ashwood等[30]的试验观测一致。

图10 颗粒接触数量变化图Fig.10 Changes in number of contact particles(a)颗粒之间；(b)颗粒与滑槽；(c)颗粒与结构

4 堆积体分析

通过对碎屑流颗粒的堆积形态进行提取，绘制颗粒堆积形态图(图11)，图中(1.25,0)为坡脚位置，(1.5,0)为导引结构放置的位置。可以看到：三种导引结构作用下，碎屑流颗粒的堆积形态相似；从最大水平运动距离上看，约为0.75 m。由此可以看出，三种不同导引结构作用对颗粒的最终运动形态影响不大。

图11 颗粒堆积形态图Fig.11 Particle packing pattern

沿着滑坡碎屑流运动方向，从坡脚处开始划分，每隔10 mm划分一个140 mm宽的区域(图12)，绘制直方图(图13)对堆积体体积分布进行研究。

图12 堆积体划分图Fig.12 Partition diagram of the accumulation body

图13 体积分布直方图Fig.13 Volume distribution histogram

如图13所示，碎屑流颗粒主要分布在坡脚处，堆积在A1区的颗粒所占百分比最大，导引结构B1作用下A1区堆积了25%，导引结构B2作用下A1区堆积了24%，导引结构B3作用下A1区堆积了19%。A2区中颗粒的百分比在导引结构B1、B2、B3的作用下，分别为20%、16%、13%。B1与B2在A1区时颗粒堆积所占的百分比相差不大；在A2区时，由于B1结构是凹型，所能容纳颗粒堆积的空间大于直线型结构B2，而结构B3是凸型，其空间最小，所以在A2区时，颗粒堆积所占的百分比出现了显著的差异。在A3区之后，三种导引结构作用下碎屑流颗粒堆积所占的百分比基本相同。碎屑流颗粒堆积在水平滑槽上的百分比随距离的增加而不断减小。通过对堆积在水平滑槽上的颗粒所占百分比进行汇总，导引结构B1、B1、B3颗粒所占百分比分别为76%、70%、64%。由此得出，导引结构的作用对于颗粒堆积体积分布会有显著的影响，主要影响区是靠近坡脚处，对导引结构之后的堆积区域的颗粒体积分布影响不显著。

5 结论

(1)对比三种导引结构，凹型圆弧结构B1可以有效地将碎屑流颗粒的作用力进行转化，结构所受的法向力最小，结构的损毁最小；切向力最大，对颗粒的导引作用最大。

(2)导引结构的存在与否，对颗粒的势能变化影响不显著，势能随时间增加不断减小，而最终有导引结构设置时，颗粒堆积在坡脚处的颗粒越多，堆积高度会增加，导致颗粒静止时的势能大于无导引结构时。经过三种不同导引结构后，颗粒与滑槽之间的碰撞和摩擦是导致颗粒动能减小的主要原因，但三种不同导引结构对颗粒动能耗散效果无显著差异。

(3)三种不同导引结构作用下，对颗粒的最终运动形态影响不大。导引结构的作用对于颗粒堆积体积分布会有显著的影响，主要影响区是靠近坡脚处，对导引结构之后的堆积区域的颗粒体积分布影响不显著。