近年来，中国西南地区地震及其次生灾害频发，流域地表及局部沟床汇集了大量的松散堆积物，在地震、降雨等因素的诱发下，大量的松散堆积物质常以崩塌、滑坡及泥石流等形式侧向汇入河道，进而导致山区河流的泥沙补给条件发生突变^[1–4]。大量破碎的山体侧向进入主河后，会挤压甚至堵断河流流路，改变河道水沙运动规律，严重威胁山区河流防洪安全^[5–7]。如2008年5月12日四川省汶川县发生了8.0级特大地震，地震造成唐家山大量山体崩塌、滑坡，大量的堆积体堵塞湔江河道，形成巨大的堰塞湖—唐家山堰塞湖，一度威胁湔江下游30多万人民的生命财产安全。

国内外众多学者通过野外观测、室内试验等研究方法对地震作用下的松散体的变形及运动特征做了大量的研究工作。在野外观测和室内振动试验方面，如黄河清等^[8]通过对文家沟滑坡进行野外调查指出滑坡的高启动速度是地震荷载作用导致的结果，与山体的碰撞是造成滑坡体破碎的主要原因。詹威威等^[9]基于汶川地震诱发的38个沟道堆积体运动的野外实测数据，建立了具有较高精度的堆积体最大水平运动距离的预测模型。刘婧雯等^[10]借助系列振动台模型试验发现基岩与边坡中的堆积体存在动力响应差异，一旦出现裂缝便会快速扩展，坡肩处的块体会在动力作用下脱离抛出。孙志亮等^[11]通过大型振动台对连续多级地震荷载作用下松散体的永久变形以及不同类型地震波对于堆积体边坡变形的影响进行了研究。Guo等^[12]通过输入正弦地震波和汶川地震波发现在垂直方向上，斜坡对于水平加速度的动力响应呈现出非线性的放大特征。到目前为止，基于连续介质与非连续介质模型的数值模拟已经成为研究松散体变形及运动过程的重要研究手段。连续介质模型把滑坡体假设为连续体，采用连续介质理论模拟滑坡运动，例如Savage等^[13]通过拉格朗日摩擦模型对干沙土的运动过程进行了数值模拟。O’Brien等^[14]基于二次流变模型利用FLO-2D计算了碎屑流、洪水和泥石流的运动过程。虽然它具有理论发展成熟、应用范围广和计算效率高等特点，但无法获取滑坡运动过程中岩体之间相互挤压断裂等细节。与连续介质模型相比，非连续介质模型既能模拟滑坡的整体运动过程，还能计算滑坡体在移动过程中内部岩体破坏随时间的变化情况，当前非连续介质模型中以离散单元法（discrete element method）应用最为广泛。He等^[15]基于2维离散单元程序对汶川地震水平剪切波作用下堆积体变坡的动力响应过程进行了计算，研究总结出低速表面堆积体是控制基岩上覆盖层运动的重要因素。孙新坡等^[16]采用离散单元法对“5·12”汶川地震引起的牛圈沟滑坡运动过程进行了数值模拟，并对基底摩擦系数对滑坡堆积形态、运动距离的影响进行了分析。

作者在众多学者的研究基础上，采用离散单元法（discrete element method）对不同地震荷载作用下松散体的变形过程及堆积特点进行了数值模拟，重点讨论了地震波频率、振幅等因素对于松散体运动速度、单元接触数量以及补给质量的影响，为深入研究地震作用下松散体破坏变形、产沙堆积及河道泥沙超量补给致灾提供了理论依据。

1）颗粒运动控制方程。颗粒的运动由平动和转动两部分所组成，根据牛顿第二定律，颗粒l的运动方程形式如下：

式中： $m$ 为颗粒 $l$ 的质量； ${{u}}$ 为颗粒 $l$ 的平动速度； $ {{g}}$ 为重力加速度； ${{\omega }}$ 为颗粒 $l$ 的转动速度； $I$ 为颗粒的转动惯量； ${n_{\rm c}}$ 为与颗粒 $l$ 接触的总颗粒数; ${{{f}}_{n,lk}}$ 、 ${{{f}}_{t,lk}}$ 与 ${{{T}}_{lk}}$ 分别为颗粒 $l$ 与颗粒 $k$ 的法向碰撞力、切向碰撞力以及碰撞力矩，碰撞力矩 ${{{T}}_{lk}} = {{{R}}_l} \times ({{{f}}_{n,lk}} + {{{f}}_{t,lk}})$ ，其中 ${{{R}}_l}$ 为碰撞点与颗粒质心的距离矢量。

2）颗粒碰撞模型。采用软球模型来描述颗粒间的相互碰撞过程。基于线性弹簧阻尼模型^[11]，颗粒间的法向碰撞力与切向碰撞力可以分别表示为：

式中， ${k_{\rm s}}$ 为刚度系数， $\eta $ 为阻尼系数， ${\mu _{\rm s}}$ 为摩擦系数， ${{\delta }}$ 为颗粒间的变形量， ${{v}}$ 为发生碰撞时颗粒间的相对速度：

式中： ${{t}}$ 为单位切向量， ${{t}} = {{{v}}_{t,lk}}/\left| {{{{v}}_{t,lk}}} \right|$ ； ${{n}}$ 为由颗粒l的质心指向颗粒 $k$ 质心的单位向量， ${{n}} = {{{R}}_l}/\left| {{{{R}}_l}} \right|$ 。

为了分析和讨论地震波频率、振幅对松散体变形及运动过程的影响，采用正弦波来模拟地震波，在计算过程中所选取的物理模型如图1所示，该物理模型由滑动面、堆积面和堆积槽3部分组成，堆积面的长度为20 m，坡高为10 m，坡比为1:2，堆积槽的几何形状是三棱主体，用于在仿真开始前存放颗粒物质，堆积槽的尺寸为20 m $ \times $ 3 m $ \times $ 3 m。在模拟过程中不考虑颗粒之间的粘结力，在设置计算参数时作者参考了文献[17]，具体的计算参数见表1，地震参数结合文献[18–19]，计算工况如表2所示。

图1(Fig. 1)

图1 物理模型示意图 Fig. 1 Sketch of physical model

表1(Tab. 1)

表1 DEM中相关计算参数 Tab. 1 Parameters used for the DEM calculation

表1 DEM中相关计算参数 Tab. 1 Parameters used for the DEM calculation 模型部件 颗粒 坡面、堆积面 密度/(kg·m–3) 2 500 — 剪切模量/MPa 100 1 200 泊松比 0.28 0.26 碰撞恢复系数 0.35 0.35 滚动摩擦系数 0.05 0.05 滑动摩擦系数 0.4 0.42

表2(Tab. 2)

表2 不同地震荷载作用下的计算工况 Tab. 2 Parameters of calculation conditions under different seismic loads

表2 不同地震荷载作用下的计算工况 Tab. 2 Parameters of calculation conditions under different seismic loads 计算工况 频率/Hz 振幅 工况1 — — 工况2 20 0.025g 工况3 20 0.05g 工况4 20 0.075g 工况5 10 0.1g 工况6 15 0.1g 工况7 20 0.1g 注： $g$ 为重力加速度， $g$ =9.8 m/s2。

为探究地震及其次生灾害对于松散体运动变形特点的影响，基于3维颗粒离散单元法（discrete eement method）对松散堆积体在坡面上的运动堆积过程进行了数值模拟。为了验证本文计算的合理性与准确性，本文首先对文献[20]中松散体的室内物理试验的两种工况进行了计算，并将计算结果与试验结果进行了比对，如图2所示。由图2可以发现计算结果与试验数据吻合较好，这也说明了DEM可以较好模拟松散体的运动过程。

图2(Fig. 2)

图2 河道中松散堆积体累积体积随时间变化关系 Fig. 2 Cumulative volume of loose deposit body varies with time in the river

图3(Fig. 3)

图3 颗粒分区示意图 Fig. 3 Schematic diagram of particle zoning

为了分析和研究松散体不同区域在运动过程中的差异，在初始时刻，将松散体划分成了5个区域，分别为前端上部、前端下部、中间上部、中间下部、后部并分别用Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ对5个区域进行标号，如图3所示。不同区域颗粒运动速度随时间的变化关系如图4所示，图5给出了不同区域颗粒在不同时刻的运动情况。

由于颗粒在运动过程中横向运动并不显著，因此这里仅给出颗粒水平 $x$ 方向速度和垂直 ${\textit{z}}$ 方向速度随时间变化情况。在松散体刚开始运动时，位于Ⅰ和Ⅱ区的颗粒最先起动，从图4（b）中可以发现，区域Ⅰ与区域Ⅱ中的颗粒在运动初期前0.5 s内， ${\textit{z}}$ 方向的运动速度迅速增加，随后曲线斜率趋于稳定，这主要是由于颗粒突然释放时，在重力的作用下松散物料前段的颗粒在运动初期存在明显的掉落现象所引起的。随着时间的推移，Ⅲ、Ⅳ、Ⅴ依次开始滑动并在重力的作用下，颗粒的运动速度不断增加，当前段颗粒运动到堆积面并与其发生碰撞后，颗粒的动能迅速减小，颗粒的运动速度快速降低，并在坡脚位置逐渐堆积。通过图4可以发现在运动过程中颗粒 $x$ 方向的速度最大（峰值速度约为4 m/s）、 ${\textit{z}}$ 方向速度次之(峰值速度越为2 m/s)，各区域颗粒 $x$ 方向和 ${\textit{z}}$ 方向的速度随着时间的推移整体上呈现先增大后减小的趋势，且各区域颗粒 $x$ 方向与 ${\textit{z}}$ 方向的峰值速度几乎相等。从图5（f）给出的是颗粒运动后的最终堆积形态，可以发现Ⅰ和Ⅱ区的颗粒几乎完全堆积在坡脚下，Ⅲ和Ⅳ的颗粒大多数堆积在坡面上，Ⅴ区的颗粒全部堆积在坡面上，并且离坡脚越近，颗粒的堆积厚度越厚。

图4(Fig. 4)

图4 不同区域颗粒运动速度随时间变化关系 Fig. 4 Relationship between particle velocity and time in different regions

图5(Fig. 5)

图5 不同时刻各区域颗粒运动情况 Fig. 5 Motion of particles in each region at different time

天然情况下松散体大多都是由宽级配的卵砾石所组成的，为了探究颗粒粒径对于运动速度的影响，本文计算时所采用的松散物料由4种粒径（5、10、15、20 cm）的颗粒所组成，且每种粒径颗粒所占的百分比相等，均为25%，图6给出了不同粒径下颗粒运动速度随时间的变化关系，从图6中可以发现随着颗粒粒径的增大，颗粒的运动速度也随之增加，且达到峰值速度的时间随之减小。图7（a）给出的是颗粒在河道中的最终的堆积形态，为了探究堆积体不同区域内的颗粒级配分布特点，这里将颗粒的最终堆积形态沿 $x$ 方向平均分成了3个区域，分别为 ${\rm x}$ 1、 ${\rm x}$ 2、 ${\rm x}$ 3，各区域内的颗粒级配如图7（b）所示。可以发现，堆积体颗粒的分布呈现出大颗粒位于堆积体表面前段，细小颗粒位于堆积体下后部的特点，这主要是因为颗粒在运动过程中会不断的发生摩擦碰撞，颗粒间存在着一定的空隙，细小颗粒运动速度较慢并且能够沿着颗粒间的空隙进入到堆积体的底部，这与自然情况下松散体的堆积规律是吻合的。

图6(Fig. 6)

图6 颗粒粒径对运动速度的影响 Fig. 6 Effect of particle size on particle velocity

图7(Fig. 7)

图7 不同区域颗粒级配情况 Fig. 7 Particle size distribution in different

作者利用正弦波来模拟地震波，重点研究地震荷载中频率和振幅对于松散体运动过程的影响。在模拟时仅考虑颗粒相互碰撞作用，不考虑颗粒间的粘结力。图8是不同地震荷载作用下松散体运动过程中单元接触数量随时间的变化关系，图8（a）给出的是频率20 Hz时不同振幅对于单元接触数量的影响，图8（b）是振幅为0.1g时不同频率对于颗粒单元接触数量的影响，从图8中可以发现在相同频率和相同振幅条件下，在颗粒运动的初始时刻，颗粒的堆积较为紧密，单元接触数量大约有37 000，随着颗粒运动的持续颗粒单元接触数量逐渐减小，在5 s附近单元接触数量减小到最小值，而后随着颗粒前段与堆积面接触逐渐开始堆积，单元接触数量逐渐增加。地震荷载导致松散体中颗粒单元的接触数量大幅度减少，当地震频率一定时，松散坡体中的颗粒的接触数量随着地震波的振幅的增加而逐渐减少；当地震振幅一定时，松散坡体中颗粒的接触数量随着地震波频率的增加而增加。图9给出的是不同地震荷载作用下松散坡体的运动速度随时间的变化规律，从图9中可以看出随着地震波振幅的增加，颗粒的运动速度逐渐增大，而随着地震波频率的增加颗粒的速度是逐渐降低的。对比图8与9可以发现在颗粒运动过程中单元接触数量随时间的变化曲线与颗粒合速度随时间的变化曲线正好相反。这说明堆积体中单元接触数量越少，颗粒越稀疏，颗粒间的相互碰撞越少，颗粒的能量损失越小，颗粒的运动速度越大。

图8(Fig. 8)

图8 颗粒接触数量随时间的变化过程 Fig. 8 Number of contacts changes with time

图9(Fig. 9)

图9 地震作用下颗粒运动速度随时间的变化 Fig. 9 Relationship between particle velocity and time under seismic loads

图10(Fig. 10)

图10 河道中松散堆积物累积质量随时间变化 Fig. 10 Cumulative mass of loose deposit body varies with time in river

图10（a）和10（b）分别表示的是地震波振幅与地震波频率对于山区河道泥沙补给质量的影响，根据图中曲线的变化可以发现地震荷载作用会导致进入山区河道的泥沙质量会大幅度的增加，并且随着时间的推移，进入河道的泥沙质量的增加速率呈现出先增大后减小的特征，地震波振幅的增加与地震波频率的降低均会导致泥沙补给速度加快，泥沙补给的必要时间减小，这主要是由于在地震荷载的作用会导致松散体的单元接触数量大幅度降低，颗粒的能量损失减小，运动速度加快，从而使堆积体可以运动到较远的位置，这使得原来停留在沟道中的松散堆积物质进入河道，河道中的泥沙补给质量陡增，这为山区河流的泥沙补给提供充足的来源，这也与地震发生后天然河流中的输沙量会大幅度增加所印证。

基于3维离散单元法，分别对不同地震荷载作用下松散体的变形及运动过程进行了数值计算，研究了地震波频率、振幅等因素对松散体运动速度、单元接触数量、补给质量的影响，主要成果如下：

1）松散体的运动过程由起动加速、高速运动以及减速堆积3个阶段所组成，不同区域颗粒速度的峰值基本相同，但达到最大值的时间不同，这主要与颗粒与堆积面的接触时间有关。颗粒的粒径越小，颗粒的运动速度越慢，颗粒进入堆积面发生堆积后大颗粒堆积在表面前段，细小颗粒堆积在堆积体的下后部。

2）地震荷载使松散体的运动速度和运动距离大幅度提高，地震波的频率越低、振幅越大，堆积体中单元接触数量越小，颗粒越稀疏，颗粒间的相互碰撞越少，能量损失越小，松散体的运动速度与运动距离越大，这会使原本未进入河道的松散堆积物进入河道，进而增大河流的泥沙补给质量。