考虑阻力作用的碎石土滑坡涌浪数值模拟方法研究

薛宏程 彭杏瑶 马倩 卿正阳

摘要：

滑坡涌浪是库区水灾害的风险来源之一，开展涌浪预测研究对提升航道安全，减小近岸洪水风险具有重要的工程应用价值，但滑坡体入水形成涌浪的动力过程十分复杂，目前对其进行准确模拟仍是棘手的问题。基于等效流体假设，将非均质的碎石土滑坡视为连续的特殊流体，通过自定义函数编写滑坡体沿滑动面运动的摩擦阻力程序，利用Realizable k-ε湍流模型与VOF法建立了碎石土滑坡涌浪数值模拟模型。研究结果表明：该方法可较好地模拟碎石土滑坡体水下运动与堆积形态，且获得的涌浪波幅精度较高；阻力模型中摩擦系数与内摩擦系数同滑坡体的变形速率分别成反比关系，而湍流阻力系数增大仅能减缓滑坡后缘的变形速率。研究成果可为定量评估碎石土滑坡涌浪灾害风险、制定库区涌浪灾害紧急避险预案提供参考。

关 键 词：

滑坡涌浪； 数值模拟； 等效流体； k-ε湍流模型； VOF模型

中图法分类号： TV139

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2023.07.022

0 引 言

碎石土滑坡是山区水库滑坡的主要类型之一，在库区沿岸分布十分广泛，主要由碎石与土的混合物组成，结构松散，在强降雨作用后极易失稳，一旦高速滑入水库、湖泊或河道后，形成的涌浪会使灾害影响区域急剧扩大，严重威胁航行船舶、水工建筑物及沿岸城镇的安全。2003年7月13日，三峡库区支流千将坪滑坡及其形成的涌浪造成13人死亡，经济损失高达8 000 万元［1］；2008年11月23日，三峡库区龚家方发生崩滑，对岸涌浪高度达13 m，涌浪传播至4.2 km外的巫山新县城时浪高仍超过1 m，使停靠码头的多只船舶受损［2］；2015年6月24日，三峡库区红岩子滑坡形成的涌浪打翻船舶13艘，造成2人死亡，4人受伤［3］。为了对此类次生灾害进行预测和早期防范，准确获得涌浪传播特征与影响范围是工程安全运行中亟需解决的问题。

数值模拟方法具有可重复性强、数据监测灵活等优势，已成为涌浪灾害预警的重要研究手段之一［4-5］。在模拟岩质滑坡冲击水体形成的涌浪时，常将滑坡体简化为刚体，再通过经验公式计算给出刚体下滑速度［6-7］。但是碎石土滑坡在运动过程中不断变形，滑坡下滑速度无法准确模拟，且滑坡入水后与水体发生强烈的掺混作用，极大增加了涌浪模拟难度。一些学者尝试采用离散元法（Discrete Element Method，DEM）将土石滑坡划分为若干个离散元单元阵来模拟其运动过程，该方法可以较真实地还原滑坡体中的非线性大变形特征。徐文杰［8］通过构建SPH-DEM流固耦合算法，再现了意大利瓦依昂（Vajont）库区滑坡涌浪灾害链的全过程；谭海等［9］引入浮力和拖曳力实现了DEM模块与SPH模块的耦合，发现此模型更适合模拟水上散粒体滑坡形成的涌浪特征。但是离散元法在提高精度的同时将面临控制参数较多、粒子分散度不能太大、计算量巨大等难点。Hungr等［10-11］基于连续介质理论，将碎石土滑坡视为一种特殊流体，认为滑坡运动形态与其内部黏性和底部摩擦阻力有关。在此基础上，朱圻等［12］、袁小一等［13］、高杨等［14］、夏式伟等［15］基于上述方法分别对牛圈沟、文家沟、鸡尾山、汤家沟发生的滑坡事件进行了模拟验证，计算结果与实际情况吻合较好。虽然上述方法能够较好地模拟碎石土滑坡的运动过程与堆积形态，但是滑坡涌浪问题的复杂性在于滑坡运动过程中除受气相影响外，还要充分考虑滑坡与水体的相互作用。针对这类典型的三相流问题如何使计算快速收敛，且能够准确模拟涌浪的形态特征，还有待进一步研究和探索。

本文基于Hungr的理论，通过自定义函数编写阻力模型程序控制滑坡运动过程，结合Realizable k-ε湍流模型与流体体积法（Volume of Fluid，VOF），对滑坡、空气与水三相共同作用下涌浪形成过程进行模拟，同时还讨论了滑坡阻力模型中摩擦系数、内摩擦系数、湍流阻力系数3个参数对滑坡形态特征的影响。

2 模型验证与结果分析

2.1 模型构建

验证算例源于Fritz在苏黎世联邦理工学院开展的散粒体滑坡涌浪模型试验［16］。该试验在长11 m、深1 m、宽0.5 m的水槽中进行，即涌浪在宽度方向上的传播可忽略不计。水槽的前端设置了长3 m且坡度固定为45°的滑坡体滑动面，滑坡散粒体通过气动装置启动，该装置可以给滑坡提供初始速度。如图1所示，试验以水面与滑动面交界点为原点O，在滑动面1号（距原点0.4 m）與2号（距原点0.07 m）断面位置安装了激光测距传感器，以获得入水前滑坡厚度随时间变化的情况，同时在①与②区域用工业相机捕捉涌浪形态。

根据上述试验建立二维几何模型如图2（a）所示，模型采用四边形网格，并对整个计算区域各相交界面附近的网格进行了加密处理（见图2（b）），网格总数约3万个。值得说明的是，本文计算前对计算区域内的网格开展了敏感性分析，发现进一步加密网格对计算精度影响较小，说明目前的网格密度足以满足计算要求。模型边界条件包括速度进口、压力进口、自由出流与固壁边界。现选取Fritz试验中的两组工况进行模拟验证，其中滑坡体密度ρs均为1 620 kg/m3，水体密度ρw为998.2 kg/m3；壁面摩擦系数f′根据Fritz的滑坡材料（PP-BaSO4）取0.445，其余参数设置见表1。

2.2 计算结果分析

2.2.1 滑坡滑动形态对比

图3（a）和图3（b）分别给出了1号与2号断面处激光测距传感器在两组试验中捕获的滑坡厚度s随时间t的变化情况，其中滑坡前缘与水面接触的时刻记为t=0 s。从图3中可以看出，工况一中数模计算得到的滑坡厚度随时间的变化情况与试验结果基本吻合，较好地反映了散粒体滑动的真实状态；工况二中滑坡体较厚且滑动速度快，滑坡体表层的变形主要受剪切应力影响，使得计算得到的滑坡变形速率略大于试验结果。不过总体来看，滑坡厚度的微小差别对涌浪形态的影响可以忽略不计。

2.2.2 近场涌浪形态对比

图4为工况一数模计算得到的涌浪形成过程与Fritz的试验结果对比（上方为实际形态，下方为模拟形态，下同），数模结果中深蓝色代表空气相，青色代表水相，红色代表滑坡相。从图4（a）中可以看出在相对时间T=0.93时涌浪形态存在一定差异，主要是由于滑坡在水下的滑动速度比试验略快；从图4（b）中可以看出计算得到的滑坡堆积形态与试验基本一致，但滑坡冲入水中使两相掺混区域的动力黏度会迅速增大，导致空气很难掺混进水相中，故无法较好地模拟出试验中白色的水气掺混现象；图4（c）中，速度梯度减小促使滑坡在水下逐渐停止运动，滑坡与水体两相边界附近的动力黏度则逐渐增大，且远远大于水体的动力黏度值，导致部分卷入水中且靠近滑坡边界的空气无法逸出，出现了一层很薄的空气滞留在滑坡表面。

由于工况二滑坡厚度与滑动速度较工况一显著增加，其形成的涌浪高度也更高，滑坡体与水体接触后呈现出明显的分离状态［17］。图5（b）与图5（c）中涌浪自由面与试验结果相比略有差异，主要体现为掺气效果不佳，出现该问题的原因依然是水体表面黏度过大导致空气难以掺入，无法模拟出水花飞溅的形态。但总体来看，掺气效果不佳和空气滞留的影响仅局限于滑坡相的边界附近，不会对涌浪的总体形态和传播衰减产生较大影响。值得说明的是，对滑坡体边界层区域的网格进行适当加密，可以使滞留气体的范围减小。

3 讨 论

改变阻力模型中的摩擦系数f′、湍流阻力系数ζ，及内摩擦系数K均会对滑坡运动特征及涌浪形态产生影响，故基于上述算例对3个参数进行单因素影响分析。以工况一为例，当湍流阻力系数ζ为定值2 000 m/s2，内摩擦系数K为定值0.5，摩擦系数f′分别取值0.1，0.3，0.5时，可得到2号断面处滑坡厚度随时间变化的情况。从图6可以看出，摩擦系数f′越大，对滑坡阻力作用越明显，减缓了滑坡的变形速率。

当摩擦系数f′取定值0.5，内摩擦系数K取定值0.1，湍流阻力系数ζ分别取100，500，1 500 m/s2时，滑坡体在2号断面处随时间变化的情况如图7所示，从图7中可以看出改变湍流阻力系数ζ对滑坡厚度变化的影响较小，但湍流阻力系数ζ增大会减缓滑坡后缘的变形速率。

当摩擦系数f′取定值0.5，湍流阻力系数ζ取定值100 m/s2，内摩擦系数K分别取0.1，0.5，0.9，从图8可以看出内摩擦系数K对滑坡变形速率影响较大。由于增大内摩擦系数K值会使滑坡内部流速梯度减小，故内摩擦系数K越大，滑坡体的变形速率就越慢。

4 结 论

（1） 本文的数值模拟方法能够较好地模拟碎石土滑坡冲击水体形成涌浪的全过程，通过获得涌浪波高、流场等信息，可以为山区水库滑坡涌浪灾害预测提供依据。

（2） 在提高涌浪模拟精度的同时，该模拟方法相比离散元法具有计算量适中、输入控制参数较少等优势。

（3） 模型中的摩擦系数f′、湍流阻力系数ζ及内摩擦系数K均会对滑坡运动特征产生影响，其中摩擦系数与内摩擦系数同滑坡体的变形速率成反比关系，而湍流阻力系数ζ的增大仅能减缓滑坡后缘的变形速率。

（4） 值得特别说明的是，阻力模型会使滑坡与水、气交界区域的动力黏度迅速增大，导致在模拟水花飞溅等强掺气现象时效果不佳，同时会使小部分气体滞留在滑坡与水的交界面处难以逸出水面，针对上述问题还需进一步优化阻力模型。此外，在实际工程中，库岸碎石土滑坡体特征千变万化，还应对碎石土滑坡的流变特性开展进一步研究，为模擬计算提供更加准确的控制参数取值范围。

参考文献：

［1］ YIN Y P，HUANG B L，CHEN X T，et al.Numerical analysis on wave generated by the Qianjiangping landslide in Three Gorges Reservoir，China［J］.Landslides，2015，12（2）：355-364.

［2］ HUANG B L，CHEN L D，PENG X M，et al.Assessment of the risk of rockfalls in Wu Gorge，Three Gorges，China［J］.Landslides，2010，7（1）：1-11.

［3］ XIAO L L，WANG J J，WARD S N，et al.Numerical modeling of the June 24，2015，Hongyanzi landslide generated impulse waves in Three Gorges Reservoir，China［J］.Landslides，2018，15（12）：2385-2398.

［4］ ALEXANDRE P，PHILIPPE H，STPHANE A.Landslide tsunamis：Comparison between depth-averaged and Navier-Stokes models［J］.Coastal Engineering，2021，170：104022.

［5］ ZHANG G，CHEN J，QI Y，et al.A WCSPH two-phase mixture model for tsunami waves generated by granular landslides［J］.Computers and Geotechnics，2022，144：104657

［6］ MAGDALENA I，FIRDAUS K，JAYADI D.Analytical and numerical studies for wave generated by submarine landslide［J］.Alexandria Engineering Journal，2022，61：303-7313

［7］ 鄧成进，党发宁，陈兴周.库区滑坡涌浪传播及其与大坝相互作用机理研究［J］.水利学报，2019，50（7）：815-823.

［8］ 徐文杰.滑坡涌浪流-固耦合分析方法与应用［J］.岩石力学与工程学报，2020，39（7）：1420-1433.

［9］ 谭海，徐青，陈胜宏，等.基于DEM-SPH耦合模型的散粒体滑坡涌浪仿真模拟［J］.岩土力学，2020（增2）：1-11.

［10］ HUNGR O，EVANS S G，BOVIS M J，et al.A review of the classification of landslides of the flow type［J］.Environmental & Engineering Geoscience，2001，7（3）：221-238.

［11］ HUNGR O，MCDOUGALL S.Two numerical models for landslide dynamic analysis［J］.Computers & Geosciences，2009，35（5）：978-992.

［12］ 朱圻，程谦恭，王玉峰，等.高速远程滑坡超前冲击气浪三维动力学分析［J］.岩土力学，2014，35（10）：2909-2926.

［13］ 袁小一，许强，程谦恭，等.高速远程滑坡-碎屑流超前冲击气浪分析［J］.水文地质工程地质，2016，43（6）：113-119.

［14］ 高杨，李滨，王国章.鸡尾山高速远程滑坡运动特征及数值模拟分析［J］.工程地质学报，2016，24（3）：425-434.

［15］ 夏式伟，郑昭炀，袁小一，等.芦山地震汤家沟滑坡-碎屑流过程模拟［J］.山地学报，2017（4）：102-109.

［16］ FRITZ H M.Initial phase of landslide generated impulse waves［D］.Zurich：ETH Zürich，2002.

［17］ XUE H C，MA Q，DIAO M J，et al.Propagation characteristics of subaerial landslide-generated impulse waves［J］.Environmental Fluid Mechanics，2019，19（1）：203-230.

（编辑：胡旭东）

Numerical simulation method of gravel soil landslide surge considering resistance effect

XUE Hongcheng1，2，PENG Xingyao1，MA Qian3，QING Zhengyang1

（1.College of River and Ocean Engineering，Chongqing Jiaotong University，Chongqing 400074，China； 2.National Engineering Research Center for Inland Waterway Regulation，Chongqing Jiaotong University，Chongqing 400074，China； 3.Southwest Waterway Engineering Institute，Chongqing Jiaotong University，Chongqing 400016，China）

Abstract：

Landslide surge is one of the risk sources of water disasters in a reservoir area.It is of great engineering application value to carry out surge prediction research to improve channel safety and reduce near-shore flood risk.However，the dynamic process of landslide body entering water to form surge is very complicated.At present，it is still a difficult problem to accurately simulate it.Based on the assumption of equivalent fluid，the heterogeneous gravel soil landslide was regarded as a continuous special fluid，and a friction resistance program of the landslide body moving along the sliding surface was written by a custom function.The Realizable k-ε turbulence model and VOF method were used to establish the numerical simulation model of gravel soil landslide surge.The results showed that this method can better simulate the underwater movement and accumulation form of gravel soil landslide，and the accuracy of surge amplitude was high.In the resistance model，the friction coefficient and the internal friction coefficient were inversely proportional to the deformation rate of the landslide body，while the increasing of the turbulent flow resistance coefficient could only slow down the deformation rate of the trailing edge of the landslide.The research results can provide a reference for quantitatively evaluating the risk of surge disaster of gravel soil landslide and formulating the emergency plan of surge disaster in a reservoir area.

Key words：

landslide surge；numerical simulation；equivalent fluid；k-ε turbulence model；VOF model