泥流在实际运动过程中，由于受到拦挡坝的拦挡作用，泥流运动行为得到控制、流速减少、侵蚀能力和流量减少，使其偏离流动方向，可以减少灾害体致灾范围和致灾规模（Chen et al., 2000, 2003）。因此，关于拦挡坝拦挡作用的研究是泥流防治中的一项有科学意义和现实意义的议题。拦挡坝作为一种被动保护措施抵御灾害体破坏，其作用不容忽视，体现在数值模拟中，即为基底高程变化的体现。伴随着科学理论与计算技术的日益发展与进步，数值模拟方法已成为解决自然界和科学界复杂动力学难题的一种必要方案。采用数值模拟技术可以对灾害事件进行反演计算，同时也能够对未来灾害事件发生的过程进行预测。尽管数值模拟方法不能取代物理模型试验，但与物理模型试验相互结合可以取得更好的机理分析结果，因此得以广泛应用（Chen et al., 2006; 马建全等，2022）。近年来，由于计算机技术的飞速进步、微分方程的改进和人们对于灾害动态特性的更加清晰的认识，泥石流、滑坡等自然灾害的影响范围和动态特性的模拟技术也在不断地进步。这些技术的应用，不仅使泥石流、滑坡等自然灾害的影响更加明显，而且也使自然灾害的动力机制和影响因素的研究变得更加有效（Hungr, 1990; Erlichson, 1991; Fannin et al., 2001; Crosta et al., 2004）。

工程实践中的地质灾害危险性评估主要涉及两个步骤。第一步是边坡失稳概率和易发性的评估；其次需要建立边坡失稳模型，以确定不同失稳概率下致灾强度的潜在分布。其中，地质灾害被动措施防治效果评估的一个重要工具是对灾害体运动行为进行数值模拟，可用于定量阐明灾害体运动过程中的危险区域。泥流数值模拟一般分为质量集中模型、独立单元模型和连续介质模型3类（Chen et al., 2000, 2003）。

独立单元法虽然可以应用于香港滑坡的研究，它利用二维颗粒流程序PARTI-2D和UDEC来模拟颗粒流的流态，以及粒子之间的交互，但却没能充分考虑到粒流、滑坡的宏观结构，比如流场的流向、流量的变化、流域的影响以及流场的不稳定性，因此，独立单元法存在着明显的局限性（Chen et al.,2000, 2003）。体现“雪橇”的质量集中模型以能量准则为依据，将泥浆的运动特性表现得更加清晰，它可以更好地捕捉到固有的内摩擦角以及孔隙水的压力，从而更好地体现出它的综合影响，诸多学者都对此做了深入的研究（Hungr, 1990; Erlichson, 1991; Fannin et al., 2001; Crosta et al., 2004）。由于该模型假设物体的形心和重心重合，导致不能考虑失稳形态下的复杂模式，仅仅将下垫面简化为光滑表面，从而无法准确反映基底高程的变化，以及灾害体在运动过程中的侵蚀和沉积作用，而这些现象都是自然界实际存在的，因此，该方法仍有一定的局限性。而遵守质量、动量守恒的连续介质力学模型，以表征灾害体动态过程，可以更为精细化的对关键动力参数进行模拟计算，再现灾害体的运动过程，因而适用性更强（Koch et al.,1994; Jakob et al., 2005）。

鉴于以上分析，本研究构建坡面泥流三维运动概化模型，引入侵蚀速率概念，基于连续介质运动模型及求解方法，模拟拦挡坝对泥流运动过程、动力行为的影响，揭示拦挡坝防护对泥流动力过程的作用机制，以期为泥流防治提供技术参考。

1.   实验材料与方法

1.1   实验设计与模型几何结构

通过数值模拟，本研究选取甘肃天水2013年群发性坡面泥流的动力特征，并分析拦挡坝的设置如何影响其动力特征。为此，笔者建立了一个基础泥流地形模型，并对其进行了拦挡坝动力机制模拟分析。模型几何结构为一倾斜的35°坡面和一个水平面，之间采用5°线性过渡区连接。x轴方向与泥流运动方向一致（图1）。试验模型坡面在y方向上的宽度为8 m （y：−4，4），x水平方向长11 m（x：−2，9），过渡区长2 m（x：9，11），水平区域长9 m（x: 11，20）。坡型采用限制坡面，坡面限制部分宽4 m（−2，2），中间部分下凹其最深（x=0），深度h=0.25 m。物源区设定为一高度（H）为0.25 m，投影半径（R）为2 m的球冠，其体积为0.3987 m³。

图  1  数值模拟实验结构示意图

Figure  1.  Set-up of the numerical experiments

下载: 全尺寸图片 幻灯片

在实际的工程应用中，拦挡坝遵循具体的地形和地貌布置。根据地形特点，本研究将拦河坝底座设置于x=12.3 m的地方。由于地形的复杂性，将拦挡坝高度调整至1 m，底座的宽度0.6 m，拦挡坝顶宽度为0.4 m，并且把迎水面与背水面的斜度调整至1:0.2。

1.2   流变模型和参数选取

非牛顿流体，如泥流、高速远程滑坡，其运行机制主要依赖于其自身的速度梯度，从而产生了摩擦力，进一步导致了剪切力的产生，这种剪切率的大小可以用来表征流体的变形，即其对剪切力的响应（Fannin et al., 2001）。通常采用摩擦模型、Voellmy模型和Bingham模型（Hungr, 1990; Erlichson, 1991; Fannin et al., 2001）描述流变特性。在这项研究中，我们使用了数值模拟方法（Friction model）来进行分析。

在摩擦模型中，更多是表征多为粗颗粒存在时基底侵蚀、摩擦、夹带的力学行为。基底有效正应力与摩擦力正相关，其关系为：

在太沙基公式中：$ \tau $表示底部的摩擦力；$ {\sigma }_{Z} $表示底部的有效正应力；$ {\gamma }_{u} $表示孔隙水压力的比率；$ \delta $表示底部的摩擦角。根据太沙基的计算，当底部的位置保持不变的情况下，孔隙水压力的加大将高效地减少底部的有效正应力位，进而减少底部的磨擦力。$ {\gamma }_{u} $和$ \delta $的值通常用$ {\varphi }_{b} $来描述，公式如下：

本研究中选取摩擦模型的流变参数为$\delta $=20°，$ {\varphi }_{b}=30° $。

1.3   侵蚀速率计算与选取

侵蚀是泥流运动过程中常见的现象。由于复杂的下垫面地形和土壤异质性，量化泥流侵蚀作用仅能事件反分析获取。学者对此开展了大量的研究工作，提出了侵蚀速率这一概念，并经过的实际的工程验证（Chen et al.,2000, 2003; 马建全等，2022）。在此引入侵蚀速率这一概念，定量表征泥流侵蚀能力。

侵蚀作用强烈时，灾害体沉积体积为初始物质体积的数倍，尽管下垫面地形变化不大。侵蚀速率（E）是一种衡量法向沟床流失和切向滑动速率之间关系的指标，它的取值范围在10⁻³之间，具体表达式如下：

式中：$ {V}_{eroded} $为侵蚀物质总体积；$ {A}_{effect} $为侵蚀面积；$ {d}_{center} $为灾害体中心的移动距离。以上参数均通过天水2013年娘娘坝坡面泥流案例测量估算（于国强等，2014）。通过对多个地质灾害的反复演算，确定修正系数$ \alpha $，本研究采用$ \alpha $=2.0（Hungr, 1990; Koch et al., 1994; Jakob et al., 2005），最终确定侵蚀速率$ E\cong 0.007 $。

2.   结果与分析

2.1   拦挡坝对泥流运动过程的影响

根据图2的数据显示，发现坡面的地貌发生了显著的改变，而且泥沙的运动也发生了改变。此外，通过云图还能够清楚的观察到最终的淤积层的厚度。通过对比图2b，泥沙的移动速率为3.0 s，并且它们的前部接触到了阻碍作用的坡底。随着时间的推移，泥流的堆积越来越厚，使得大量的淤积物无法通过堤坝的抵抗。甚至当遇到严重的泥流时，堤坝的存在仍是必须考虑的因素。根据图2，即使泥流的堆积层只占整个主体的50%，堤坝的抵抗力仍不足，仍会导致泥流的扩散。

图  2  有无拦挡坝泥流运动模拟结果

a. 无拦挡坝；b. 有拦挡坝

Figure  2.  Simulation results of motion of debris flow with and without dam

下载: 全尺寸图片 幻灯片

通过对比图3和图4，可以清楚地观察到，当泥流碰到拦挡坝时，堤岸上的沉淀物会迅速扩散，其中沉积物的最高深度高达0.2 m，比未被阻碍的时候高出0.1 m，而沉积物的分散范围也会更广，这表明拦挡坝此时已被泥流的冲刷所覆盖，阻止了不断上升的沉积物。通过速度矢量的变化可知，即便没有拦挡坝的阻拦，泥流仍会沿着x轴方向增加，这表示它的速度仍然相对较快。但是，如果泥流漫过拦挡坝时，它的速度会急剧下降，且速度矢量十分稀疏，这表示泥流的运动路径发生变化。

图  3  有无拦挡坝泥流运动速度矢量及堆积区厚度分布（撞击时刻）

a. 无拦挡坝；b. 有拦挡坝

Figure  3.  Distributions of moving speed vector and thickness of accumulation zone of debris flow with and without dam (impact moment)

下载: 全尺寸图片 幻灯片

图  4  有无拦挡坝泥流运动速度矢量及堆积区厚度分布（漫坝时刻）

a. 无拦挡坝；b. 有拦挡坝

Figure  4.  Distributions of moving speed vector and thickness of accumulation zone of debris flow with and without dam (overtopping moment)

下载: 全尺寸图片 幻灯片

当没有拦挡坝的情况下，图5和图6显示了泥流进入水平区域撞击拦挡坝后，堆积体的形状发生了巨大的变化，由于受到拦挡坝的阻挡，泥流的横向扩展增强，从而导致堆积体的厚度和分布范围明显增加，从而使得坝前和坝后的堆积区域明显扩散。随着拦挡坝的出现，泥流的纵向延伸受到了明显的抑制，堆积区的中心和最大深度位置处的滑移距离显著缩短，从而使得灾害的影响范围也大大缩小。同时，从图6的坡面泥流运动轮廓模拟结果也可以看出，在有拦挡坝的作用下，坡面泥流前缘已经明显后退，泥流的主体部分已经大部分后退至过渡区位置，仅有少部分停留在坝后位置，拦挡坝已经拦截了泥流灾害体的大部分物质，致灾范围和致灾物质已大部分减少。且由于拦挡坝的拦挡作用，泥流堆积体横向扩展均堆积在坝前，使得泥流的运动方向发生了改变。

图  5  有无拦挡坝泥流堆积区厚度分布情况（堆积时刻）

a. 无拦挡坝；b. 有拦挡坝

Figure  5.  Thickness of accumulation zone of debris flow with and without dam (accumulation moment)

下载: 全尺寸图片 幻灯片

图  6  坡面泥流运动轮廓模拟结果

a. 无拦挡坝；b. 有拦挡坝

Figure  6.  Projection of debris profiles on the horizontal plane

下载: 全尺寸图片 幻灯片

2.2   拦挡坝对泥流动力参数的影响

图7和图8显示，当没有拦挡坝时，泥流动态特性会发生显著改变，其中，前端速度、平均速率和动能都会达到峰值，拦挡坝的作用会迅速起效。在有无拦挡坝作用下，泥流运动过程中的各个动力参数的变化趋势基本保持一致（图7）。泥流沿坡面运动过程中，其各个动力参数（平均速度、前端速度、总动能）均表现出迅速增加趋势，在坡脚位置处达到峰值；当泥流进入水平面后，动力参数呈现下降态势。同时，前端速度在下降过程中呈现出一定的波动态势，这是由于泥流前端在接触到坡脚时，受到坡脚的挤压与惯性作用的双重影响所致，而且在有拦挡坝的情况下，波动趋势更为明显。总势能则表现出平稳下降，然后保持稳定的态势，在有拦挡作用下，部分泥流物质堆积于拦挡坝上，总势能最终略高于无坝状况。

图  7  有无拦挡坝泥流动力参数随时间变化模拟结果

a. 平均速度；b. 前端速度；c. 总动能；d. 总势能

Figure  7.  The changes in dynamic parameters of debris flow with time with and without dam

下载: 全尺寸图片 幻灯片

图  8  有无拦挡坝泥流面积、侵蚀物质体积模拟结果

a. 泥流面积；b. 泥流体积

Figure  8.  The changes in area and volume of debris flow with time with and without dam

下载: 全尺寸图片 幻灯片

从有无拦挡坝的平均速度、前端速度、总动能的对比分析可以看出，在拦挡坝加持的作用下，3个动力参数的运动过程加快，均表现出压缩状态，表明泥流的运动过程收到限制。同时，由于拦挡作用，使得泥流整体的平均速度、前端速度和总动能均有所减少。平均速度平均降低10.62%，最高可减少15.65%；总动能平均削减16.17%，最高可削减22.89%，可在一定程度上减少泥流致灾规模。

通过图8中泥流灾害体面积和体积的变化可知，有无拦挡坝时变化趋势一致。泥流灾害体面积随运动过程逐渐增加，然后稍有回落，这是由于灾害体后缘的堆积和挤压，使泥流面积稍微减少。泥流灾害体体积随运动过程逐渐平稳增加，并维持稳定。但在拦挡坝的作用下，泥流灾害的面积和体积较无拦挡坝时均有一定的下降，平均下降2.48%和3.63%，说明拦挡坝在经受泥流冲击的同时，拦挡作用降低了泥流物质的增加。

3.   讨论

当泥流漫坝时，泥流主体迅速在坝前横向扩散，大部分泥流堆积于坝前；相反，泥流的纵向延伸逐渐减弱，泥流的厚度也会逐渐下降，致灾范围得以减小。

可以说，时间和空间的关系决定了泥流前端轮廓，图4展示了泥流前端碰撞拦挡坝时刻（t=3.0 s）的瞬时速度分布，此时，泥流前端的速度处于3～5 m/s之间；但当泥流漫坝时，其速度会急剧减小，范围为0.5～1 m/s。拦挡坝的作用是极其重要的，它能够有效地控制泥流的高度和前端速度，从而使泥流的运动变得更加平稳，从而显著地提升了泥流的运动速度。同时，结合泥流运动过程中撞击时刻和漫坝时刻的速度矢量分布（图3、图4）以及泥流平均速度和前端速度可以看出，有拦挡坝时的前端速度与无拦挡坝时的变化不大，表明拦挡坝此时正经受泥流的冲击，尽管总体数值变化不大，但此时拦挡坝的作用已经使得泥流前端速度方向发生了改变，使其偏离主流动方向，变为向上、向下或向两端扩散（图4、图6），使得泥流在纵向（x方向）上的前端速度在很大程度上已有所降低，可以在很大程度上减少了泥流的冲击力（陆鹏源等，2016），进一步降低了泥流的致灾强度和致灾规模，这与之前诸多学者数值模拟的研究结果基本一致（Hungr, 1990; Koch et al., 1994; Gray et al., 1999; Fannin et al., 2001; Jakob et al., 2005）。

由于拦挡坝的拦截作用，降低了泥流侵蚀夹带作用，使得泥流灾害体的面积和体积有一定程度的减少，拦挡作用降低了泥流物质的增加（图8），在一定程度上减少了泥流灾害体质量。在泥流灾害体运动速度（平均速度减少10.62%）和灾害体质量（面积和体积平均减少2.48%和3.63%）减少的双重叠加作用之下，导致泥流灾害体的总动能（减少16.17%）进一步降低，最终削减了泥流的致灾能量，这与Chen等（2000, 2003）采用数值模型计算结果相类似（王玉峰等，2021），也验证了此次数学模型的准确性。

经过计算，笔者发现拦挡坝的设计能够一定程度上阻止泥流的运动。但是，即使这种工程能够有效地阻止大规模泥流的运动，未来仍然可能出现漫坝的情况。因为，即使采用相同的触发机制和沟道条件，也可能会出现极端的情况。

4.   结论

（1）当泥流撞击拦挡坝堤造成漫坝时，会对拦挡坝造成重要的影响。堆积体由于拦挡坝的作用，它会在横向上不断膨胀，并且会逐渐朝着拦挡坝的方向变薄，从而使淤积层的厚度逐渐变小，限制淤积层的破坏能力。

（2）当没有拦挡坝的情况下，泥流的前端速度几乎没有变化，但是当坝前坝后的泥流增加时，速度矢量明显减小，而且分布也比较稀疏，这说明拦挡坝正在承受着泥流的持续冲击。由于拦挡坝的拦截和拦挡作用，降低了泥流侵蚀夹带作用，减少了泥流灾害体质量，改变了运动方向，减少了10.62%平均速度，使得泥流在纵向上的前端速度在很大程度上已有所降低；在双重作用叠加之下，降低了泥流运动的总动能，最终削减了灾害体16.17%的致灾能量，从而减少了泥流的冲击力，进一步降低了泥流的致灾强度和致灾规模。

（3）通过数值模拟实验，可以更好地研究实际的工程情况，特别是在研究过去的泥流灾害时，如何应对地形、水力特性、水位变化以及水土流失的复杂性。通过数值模拟方法，可以获得更准确的模拟结果，从而有助于更有效的识别、预测、控制、管理地质灾害。