考虑坝后淤积的泥石流冲击拦挡坝动力响应研究

王东坡，瞿华南，沈 伟，何思明

(1.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059；2.博洛尼亚大学 生物地质与环境科学系，博洛尼亚 40126；3.中国科学院 成都山地灾害与环境研究所，四川 成都 610041)

“5·12”汶川特大地震诱发数以万计的滑坡崩塌地质灾害，并在地震灾区新增约1×1010m3的松散固体物[1]，为震后泥石流的孕灾提供了丰富的物源，泥石流特征较震前发生了显著变化，其趋势由震前的低频中小规模泥石流向高频大规模泥石流转化[2-3]。为此，修砌各类防治工程，如拦挡坝、停淤厂、排导槽等，然而，近年来泥石流防治工程的失效案例在西部地区时有发生[4-6]。经过对泥石流防治工程失效案例的现场调查，发现其中重要原因之一是泥石流拦挡坝、谷坊坝、排导槽等防治工程在泥石流反复作用下被淤积填满，失去部分或全部减灾能力。因此，研究坝后淤积对泥石流动力特征的影响具有重要意义。

目前，针对泥石流坝后淤积的研究主要集中在泥石流沟道淤积形态以及淤积坡度的规律方面。赵静静等[7]通过模型实验，得到了坝后回淤比降与沟床比降及坝高的相关关系；韩文兵等[8]探讨了梳子坝淤砂的形态特征、淤砂长度和回淤坡度；黄海等[9]基于冲淤总量和颗粒参数建立了不同沟道段泥石流容重经验计算公式；Yang等[10]在实地调查结果和理论分析的基础上，建立了淤积坡度，淤积长度和大坝储蓄量的计算方法；孙昊等[11]在拦沙坝荷载组合与稳定性分析的基础上，推导空库工况和半库工况黏性泥石流过流时坝体稳定性系数的表达式；Gao等[12]模拟泥石流冲击城市建筑物，发现建筑物前碎屑淤积会增大冲击力；Liu等[13]基于Kanako对泥石流防治工程进行评估，发现已被先前泥石流沉积物填满的坝体，仍然具有控制泥石流的功能；贾世涛等[14]通过改变一次过程体积总量，进行半库和满库的泥石流冲击物理试验；陈洪凯等[15]以被泥石流淤埋的桥梁为研究对象，推导了泥石流沉积物流变固结的力学公式；何思明等[16]假设沟床以圆弧面侵蚀，研究了黏性泥石流运动对淤积体的侵蚀启动机制。上述众多研究成果并未考虑坝后淤积对泥石流动力过程的影响，然而，经过泥石流频繁作用后其防治工程部分或全部淤满，它们仍具有多少减灾效果？迫切需要一套科学、合理的量化评价方法或平台。

鉴于网格方法的局限性，非连续介质模型的无网格方法已逐渐应用于泥石流的数值模拟。其中SPH是一种纯拉格朗日方法，在处理自由表面流动和大变形问题方面被广泛认为优于传统的数值方法[17]。该方法不需要繁琐的网格划分，克服了传统网格方法中局部网格畸变及网格重划分的问题。基于上述优点，近年来，国内外学者采用SPH的方法针对洪水、滑坡和泥石流等自然灾害的动力学问题开展了深入研究。韩征等[18]基于SPH研究HBP本构模型在泥石流模拟中的应用。Chen等[19]采用SPH-FEM耦合方法研究考虑坝体破坏对泥石流的运动的影响。Cuomo等[20]基于SPH分析了坝体位置数量对流体流速、厚度的影响。Dai等[21]应用SPH方法探讨了泥石流运动过程中坝体冲击力的演变。胡嫚等[22]基于弹塑性理论将SPH应用于滑坡模拟。SPH方法的快速发展，对本项目的研究具有重要意义，因此，作者拟在光滑粒子流体动力学的基础上，采用HBP模型描述泥石流流变特性，基于DP准则判断坝后淤积物的侵蚀，开展坝后淤积条件下泥石流动力过程的研究，并分析坝后淤积体对拦挡坝动力响应的影响。

1 基于SPH算法的泥石流冲击拦挡坝理论计算模型

1.1 泥石流控制方程

SPH是拉格朗日无网格方法，各个粒子上承载多种物理量，包括质量、速度等，通过求解粒子的动力学方程和跟踪每个粒子的运动轨道，求得整个系统的力学行为。泥石流流体的运动可采用纳维-斯托克斯（N-S）方程描述，在SPH中N-S方程的离散形式[23]为：

式中：ρi为第i个粒子的密度；mj为第j个粒子的质量；να为速度分量；P为正应力；g 为重力；Wij为核函数（Wij=W（x-x′，h）），核函数有3种选择：2次函数、3次样条函数和Wendland5次核函数。通常情况下，SPH中积分近似的精度随核函数次数增加而更好，但其计算效率则反之。Gomez-Gesteira等[24]对SPH方法中比较常用的光滑函数特性进行了对比分析，指出Wendland 5次函数在计算精度和计算效率上有一个很好的折衷。因此，作者使用了Wendland核函数。τ为与应变张量相关的剪应力张量，其表达式如下：

式中：μeff为有效黏滞系数，用以描述流体的流变特性；ε为应变张量。

目前学者提出了多种本构模型描述泥石流流变特性，其中Herschel-Bulkley模型和双线性流变模型较为复杂，在实际中很少使用，Bingham模型对其流变特性进行描述并取得一定满意的结果，但它不能近似描述冲刷过程中遇到的应力状态范围。为避免这一问题，采用Herschel-Bulkley-Papanastasiou（HBP）模型[25]对泥石流的流变特性进行了建模。HBP模型及其等效黏性系数表示为：

式中：γ为剪切应变率；μ为表观动态黏度；τy表示屈服应力；m控制应力的指数增长，n是一个幂律指数，可以模拟剪切变薄或剪切增稠行为。当m=0和n=1时，模型简化为牛顿模型，而当m→∞和n=1时，模型简化为Bingham模型。

最终SPH动量方程表示如下：

1.2 淤积体侵蚀准则

坝后淤积体受到重力、黏聚力和相互间的内摩擦力，具有一定的抗冲刷能力。需要屈服准则来提供淤积体剪应力的临界值。只有当上游泥石流达到其临界值，沟床中的物质才能开始流动并被泥石流带走，如图1所示。

图 1 泥石流与坝后淤积体相互作用示意图Fig. 1 Schematic diagram of the interaction between debris flow and deposition behind the dam

根据Fourtakas等[26]对Mohr-Coulomb和Drucker-Prager两种侵蚀启动条件的研究，发现由前者得到的发生侵蚀的底床颗粒层等效黏度偏大，从而导致侵蚀层厚度比实验结果小。因此，采用Drucker-Prager屈服准则。

综上，将HBP本构模型与DP准则在SPH框架下集成，实现考虑坝后淤积下泥石流动力过程的模拟。

2 物理模型试验验证

为验证上述数值模型的可靠性，结合室内泥石流冲击实验结果进行校验，试验装置如图2所示。泥石流物理模型试验可模拟泥石流从启动、运移、冲击拦挡坝的全过程。水槽顶部为供料池，实验浆体颗粒级配参考红椿沟泥石流颗粒级配进行配制。拦挡坝安装在水槽的末端，坝体上安装有压力传感器，便于测量坝体受到的泥石流冲击力。在水槽的下方设置回收池，方便回收循环利用泥石流浆体。

图 2 试验水槽示意图Fig. 2 Sketch of laboratory flume and measuring equipment

数值计算按照实际3维尺寸进行建模，试验中对密度ρ，抗剪强度参数 φ、黏性系数μ等关键物理量进行测定。SPH流体粒子初始间距dp=0.005 m，光滑长度按h=(3(dp)2)1/2确定，采用Verlet时间积分算法，模拟流体4 s内的流动状态，模拟关键参数见表1。

表1 SPH方法模拟关键参数Tab.1 Parameters used in the simulation

图3为不同时刻泥石流运动过程，结果较好地反映了泥石流触坝、爬高、最终淤积的过程。

图3 泥石流冲击过程模拟与实验对比结果Fig.3 Comparison of simulation and experimental results of debris flow impact process

图4提取了a0、a1和a2处的冲击压力与试验数据进行对比，从图4可以看出，HBP模型可以较好地捕获峰值压力和压力的剩余峰后值，测试结果和模拟结果的峰值和波形吻合度较好。由此可见，基于SPH所建立的数值模型具有可行性，能够较好模拟泥石流的运动过程并预测冲击力。

图4 冲击力试验结果与数值模拟结果Fig.4 Impact force-time test and numerical simulation results

3 实例应用

在第2节模型可靠度验证的基础上，进一步对汶川地震灾区的红椿沟泥石流动力过程开展数值模拟研究。总共模拟了3种工况，包括空库过坝、半库过坝和满库过坝，从而分析在坝后不同淤积条件下泥石流冲击拦挡坝的动力响应。

红椿沟（图5）距离映秀镇东北方向约500 m，该沟位于“5·12”汶川大地震的震中区。红椿沟流域的总面积为5.35 km2，主沟长度约3.6 km，沟道高程在880 m至1 700 m之间，平均坡度约36°。红椿沟的上游包括3个分支：从西到东分别是甘溪铺沟、大水沟和新店子沟。2010年8月12日至14日，该地区经历一次强降雨事件，降雨最终引发了该次泥石流事件。根据野外调查[27]该次泥石流启动量约700 000～800 000 m3，其中约一半的泥石流物质（350 000～400 000 m3）冲出了沟口并阻塞了岷江，导致岷江改变流向，对灾后重建的映秀镇造成了极大破坏，导致17人死亡或失踪。为防止灾害再次对人民生命财产造成损失，红椿沟泥石流治理工程[28]于2011年5月竣工，主要布置工程为上游3条支沟修建4座拦挡坝、15座谷坊坝，中游修建4座较高的格拦坝。红椿沟防治工程竣工至今，红椿沟发生2次中型泥石流、2次小型泥石流以及多次洪水。在工程防治下，泥石流未形成规模、冲出沟口，防治工程取得一定的减灾成效。然而，多次的泥石流灾害使得红椿沟拦挡坝坝后淤积严重，减灾效果已然不如建成时期。主沟中下段布置1#、2#格拦坝，利用其库容优势主要起拦挡作用，从而在坝后产生淤积，因此本文结合原主沟1#、2#及排导槽处设置拦挡坝旨在研究淤积作用下拦挡坝的动力响应过程，坝体高度结合实际工程统一为20 m。

图5 红椿沟渠流域全景Fig.5 Panorama of Hongchungou drainage basin

从无人机获得红椿沟流域数字高程模型（digital elevation model，DEM），由原始DEM数据导出的三角地形网格，并将地形网格转换为一系列边界粒子，而流动材料也被离散为一系列具有特定性质的粒子。模拟泥石流初始体积约700 000 m3，淤积体布置在S1拦挡坝坝后，满库工况下淤积总量基于周必凡等[29]提出的经验公式确定：

式中：A为坝址处沟道的横断面积；ls为回淤段长度；K为经验系数，一般取值为0.3～0.5，本文取0.4。

粒子距离设置为3 m，最终生成684 824个边界粒子和31 349个流体粒子，满库工况下生成2 488个淤积粒子，半库工况生成1 200个淤积粒子。参考黄雨[21]、Ouyang等[30]开展红椿沟泥石流相关数值模拟，本次研究采用的主要参数见表2。

表2 模拟关键参数Tab.2 Parametersused in the simulation

图6为泥石流在空库和满库工况下的运动过程，粒子的颜色代表泥石流的流速。在t=60 s时泥石流运动至第一座拦挡坝，此时流速约25 m/s，随后泥石流与拦挡坝相互作用流速减小。在t=150 s时可以看出满库工况下泥石流运动距离（约150 m）相比空库工况更远，表明坝后淤积体将减小泥石流与拦挡坝相互作用时间。泥石流最终淤积的范围如图7所示。泥石流主要堆积在大坝的后面，并且向沟道上游延伸，此时空库工况下冲出沟口体积约198 000 m3，满库工况下冲出沟口的体积约330 000 m3。

图6 泥石流的运动过程Fig.6 Runout processof the debris flow

图7 不同工况下泥石流最终淤积Fig.7 Simulated deposition processes with different conditions.

为进一步分析在淤积条件下泥石流与拦挡坝的相互作用过程，图8展示了泥石流冲击拦挡坝S1的详细过程。对于空库工况（图8（a）），泥石流在t=61 s到达拦挡坝，此时泥石流具有较大的速度并在较短的时间内与拦挡坝坝体产生强烈的相互作用，部分泥石流被拦挡坝拦截在坝后形成死区，部分流体反弹与后面的流体相互碰撞；在t=69 s泥石流爬升到最高点并越过大坝继续行进，拦挡坝所受到的冲击力峰值为1 690 kPa主要由动压力决定。对于半库工况（图8（b）），t=61 s时泥石流侵蚀淤积体并夹带淤积体向前运动，并在t=64 s时泥石流发生垂直喷射现象。图8（c）展示了满库工况泥石流与坝体相互作用的过程，主要表现为泥石流夹带淤积物质向前运动，在t=68 s时泥石流已翻过坝体。此外从t=80 s可以看出，同一时刻下随着坝后淤积体增大过坝后流量明显增加。

图8 泥石流与拦挡坝相互作用Fig.8 Interaction between debris flow and check dam

泥石流流速是泥石流的主要特征之一，图9显示了泥石流到达第一座坝后流速（v1）和坝前出口的流速（v2）的变化。由图9可以看出，在空库工况下，泥石流达到拦挡坝流速峰值明显减小。为评价泥石流流速的减小程度，以泥石流到达坝后流速和坝前出口的流速之差与坝后流速的比值作为该结构对泥石流流速削减程度的指标，即定义速度减小率T=(v1-v2)/v1×100%。计算结果显示，当泥石流空库、半库、满库过坝时，速度的减少率分别为31.6%、16.4%、6.5%，表明随坝后淤积体的增加拦挡坝的削峰效应明显减低。

图9 过坝前后速度变化Fig.9 Velocity changesbefore and after thedam

图10通过比较不同工况下坝体不同位置的冲击力变化，分析坝后淤积体对泥石流拦挡坝冲击动力响应的影响。

图10 不同位置处冲击力时程曲线Fig.10 Relationship of impact force-time at different positions

根据图10可直观看出，3种工况下拦挡坝受到的冲击力均表现为随着时间增大迅速增大，随后冲击力趋于平缓，造成这个现象的原因有两个：1）泥石流流量不再随着时间继续增加；2）拦挡坝后的冲击力不能直接作用于坝体，主要为淤积在坝后的土压力。根据图10（c）在空库工况下拦挡坝坝顶受到的冲击力为323 kPa，没有明显的峰值；而在满库工况下峰值冲击力为890 kPa。空库和满库工况下拦挡坝顶部冲击力的变化体现了坝体动力响应过程明显的区别，在空库过坝时坝顶主要受到土压力控制，泥石流冲击力方向主要为竖向方向，坝体的垂直分量较小；而满库过流冲击力先达到峰值再逐渐下降，这是刚开始坝顶受到泥石流直接作用，随后由于顶部淤积物质厚度较小，造成土压力没有刚开始冲击力大，所以呈现下降趋势。值得注意的是在图10（a）中，在半库工况下坝体底部受到的冲击合力最小为1 466 k Pa，有两个因素可能导致此结果：一是当坝体后面有淤积体时，由于淤积体与泥石流的相互作用，增大泥石流运动的阻力，使得泥石流到达坝体前的速度相对空库时较小；另一个原因是淤积体的存在导致泥石流的冲击力存在分量，因此半库过坝时拦挡坝受到的冲击力最小。满库过坝时拦挡坝所受到的冲击合力最大为1 913 kPa，这是因此初始的主动土压力较大所导致的，空库过坝时的冲击力介于两者之间为1 690 kPa。

4 结 论

以红椿沟灾害为背景，通过光滑粒子流体动力方法（smoothed particle hydrodynamics，SPH）研究了考虑坝后淤积条件下泥石流冲击拦挡坝的动力响应机制，得出下列结论：

1）依托水槽物理模型试验，利用SPH方法对试验过程进行数值模拟。通过理论结果和试验结果的对比分析表明，基于HBP模型的SPH方法能够较好的模拟泥石流动力过程。

2）坝后淤积显著改变了泥石流的流动特性和拦挡效果。在考虑坝后淤积条件下泥石流运动距离更远，满库工况下冲出沟口的泥石流体积为空库工况的1.67倍。随着坝后淤积体增加，空库、半库、满库工况的速度减小率分别为31.6%、16.7%、6.5%，拦挡坝削峰效应逐渐降低；但在满库工况下，坝后淤积物质使得沟床坡度减小，拦挡坝对泥石流仍具有一定的调控作用。

3）坝后淤积通过影响泥石流与拦挡坝相互作用过程，从而改变坝体动力响应过程。在满库工况下拦挡坝坝顶冲击力受泥石流动压力控制，其峰值相比空库工况增大1.75倍，从而增加拦挡坝坝体发生倾覆破坏的风险；在半库工况下，拦挡坝底部受到的冲击合力峰值相比空库工况减小13.25%，满库工况相比空库工况底部冲击合力峰值增大13.19%。考虑实际中拦挡坝坝底容易发生破坏，可采取在坝底放置块石的措施保护坝体，同时当坝体淤满应及时进行清淤工作。