CFD–DEM Coupling Numerical Simulation of the Seyu Landslide–Surge Disaster Chain in Houziyan Reservoir
-
摘要: 滑坡涌浪是一种常见的地质灾害现象，由于滑坡体与水体之间存在复杂的流固耦合作用，使得传统的单一介质模型无法进行准确求解。为此，介绍一种基于计算流体力学方法（CFD）与离散单元法（DEM）的流固耦合模型CFD–DEM，采用计算流体力学方法（CFD）求解水体流动，采用离散单元法（DEM）模拟散粒体滑坡运动，充分利用不同计算模型的优势，对滑坡及涌浪演进过程进行数值模拟分析。首先，利用该耦合模型对Robbe–Saule开展的颗粒堆积体坍塌–涌浪试验进行了相同工况下的数值计算，从颗粒坍塌运动过程、涌浪高度演化过程等方面进行对比，结果表明模拟结果与试验结果吻合很好，验证了CFD–DEM流固耦合模型的有效性。然后，将该方法应用于四川省猴子岩水库色玉滑坡–涌浪灾害的演进过程分析，重现了该事件滑坡失稳运动、涌浪产生及传播、涌浪爬升、涌浪回流的全过程，计算结果显示：计算得到的电站进水口处涌浪高度与实测数据较为接近；色玉滑坡从失稳运动至静止堆积的持续时间约为20 s，颗粒平均速度最大达到16.12 m/s；滑坡引起的涌浪约在滑坡失稳10 s后传播到对岸，之后开始沿坡面向上爬升，最大爬升高度达到27.32 m。研究表明CFD–DEM流固耦合模型能够很好地应用于模拟山区河谷大规模滑坡涌浪灾害，可为库区防灾减灾提供高效的技术支持。
-
关键词:
- 滑坡涌浪 /
- 计算流体力学（CFD） /
- 离散单元法（DEM） /
- 耦合模型 /
- 数值模拟
Abstract: Landslide surge is a common geological disaster phenomenon. Due to the complex fluid-solid interaction between landslides and water, the traditional single-medium model cannot be solved accurately. Therefore, based on a coupling model of computational fluid dynamics (CFD) and discrete element method (DEM), the CFD method was used to simulate the water flow, and the DEM was used to simulate the motion of a granular landslide. Taking full advantage of the advantages of different computational models, the evolution process of landslide and surge was numerically simulated and analyzed. Firstly, this method was used to numerically simulate the particle accumulation collapse-surge experiment carried out by Robbe-Saule under the same working conditions. The particle motion process and the surge propagation process showed that the simulation results were in good agreement with the experimental results, which verified the effectiveness of the CFD–DEM fluid-solid coupling model. Then, the method was applied to analyze the evolution process of Seyu landslide-surge disaster in Houziyan reservoir of Sichuan Province, and the whole process of landslide movement, surge generation and propagation, surge climbing on the opposite bank and surge backflow was reproduced. The calculated surge height at the inlet of the power station was close to the measured data. The calculation results showed that the duration from unstable movement to static accumulation of Seyu landslide was about 20 s, and the maximum average particle velocity reached 16.12 m/s. The surge caused by the landslide spread to the opposite bank after about 10 s and then began to climb along the slope, with the maximum climbing height reaching 27.32 m. The research showed that the CFD–DEM fluid-solid coupling model could be well applied to the simulation of large-scale landslide and surge disasters in mountainous valleys, which could provide efficient technical support for disaster prevention and mitigation in the reservoir area. -
位于河道、水库周边的大型滑坡、崩塌等地质灾害一旦发生，可能造成堵江并产生巨大的涌浪，对人民的生命财产安全、船舶安全及水电站的运行等构成威胁[1-3]。例如：2003年7月，位于三峡库区的千将坪地区突然发生了一起方量为1.5×107 m3的滑坡，该滑坡形成的堰塞坝将青干河堵塞，造成24名人员死亡，直接经济损失约700万美元[4]。2013年7月，云南省永善县黄华镇黄坪村发生山体滑坡，滑坡体冲入金沙江激起15～20 m高的涌浪，造成12人失踪，多部车辆损毁[5]。2015年6月，重庆市巫山县红岩子发生滑坡，其方量约为2.3×105 m3，坡体入水后产生5 m多高的涌浪，造成2人死亡，6人受伤，13条船只毁坏[6-7]。
数值模拟方法是国内外开展滑坡–涌浪灾害相关研究的重要手段[8-9]，可较为快速、准确、全面地分析滑坡–涌浪动力学过程。然而，由于滑坡入水产生涌浪的过程伴随滑坡体与水体之间的流固耦合作用，较为复杂，且空间尺度较大，很难用单一介质模型进行求解。考虑到流体和滑坡的不同力学性质，应采用不同的数值模拟方法。在处理离散颗粒时，离散单元法（discrete element method，DEM）具有很好的灵活性，能够很好地反映滑坡体滑动、碎裂的过程，因而被广泛应用于滑坡灾害动力过程的数值模拟[10-12]；计算流体力学方法（computational fluid dynamics，CFD）在水体演进过程的计算精度方面有其独特优势[13]。通过两种方法间的耦合求解，能改善传统模拟方法在研究流固耦合现象时的不足，考虑滑坡运动过程中库水与滑坡体间的水土相互作用，还原水库滑坡在水环境影响下的真实演进过程，在模拟时间和空间尺度都较大的真实滑坡涌浪灾害时具有一定可靠性。
CFD–DEM流固耦合模型已被广泛应用于化学分析、工业制造等领域，目前在地质灾害领域也初步开展了一些应用。Li等[14]使用具有远端相互作用的键合粒子网络模拟柔性防护网，利用计算流体力学和离散元耦合研究泥石流对柔性防护网冲击的动力响应。Park[15]针对溃坝过程中的固液气混合流动，研究了颗粒密度对自由表面和颗粒运动的影响，并对运动特征进行了分类。Shi等[16]应用CFD–DEM的耦合方法研究滑坡坝的渗流特性，定量分析了粗颗粒与细颗粒不同的渗流破坏机理。Zhao等[17]采用CFD–DEM耦合程序对平面应变条件下Vajont滑坡过程进行了2维分析研究，并将模拟结果与ALE–FEM方法获得的结果进行了比较。
本文利用开源程序CFDEM、OpenFOAM及LIGGGHTS，基于计算流体力学方法与离散单元法构建流固耦合模型[18-20]；将该耦合模型应用于真实情况下3维大尺度滑坡–涌浪灾害动力过程演进模拟；以四川省猴子岩水库色玉滑坡为例，对滑坡涌浪灾害的动力学过程进行分析，一次解决滑坡变形失稳及涌浪形成与传播的问题，为山区大区域尺度大型滑坡–涌浪复合灾害的分析研究提供新的技术手段。
1. CFD–DEM耦合模型
1.1 计算流体力学控制方程
CFD–DEM耦合模型中，流体包含水体与空气两相，采用VOF模型求解多相流的界面问题。在VOF模型中，流体界面的演化通过求解以下方程来描述：
$$ \frac{{\partial {\alpha _i}}}{{\partial t}}{\text{ + }}\nabla \cdot ({\alpha _i}{{\boldsymbol{u}}_{\text{f}}}){\text{ = }}0 $$ (1) 式中，
$ {\alpha _i} $ 为流体i的体积分数，$ {{\boldsymbol{u}}_{\text{f}}} $ 为流体速度，$t$ 为时间。考虑颗粒相的影响，添加人工对流项以保证界面尖锐后，方程如下：$$ {\;\;\;\;\;\;\frac{{\partial {\alpha _{\text{p}}}{\alpha _i}}}{{\partial t}}{\text{ + }}\nabla \cdot ({\alpha _{\text{p}}}{\alpha _i}{{\boldsymbol{u}}_{\text{f}}}) - \nabla \cdot ({{\boldsymbol{u}}_{\text{c}}}{\alpha _i}(1 - {\alpha _i})){\text{ = }}0 }$$ (2) 式中，αp为孔隙率，
$ {{\boldsymbol{u}}_{\text{c}}} $ 为压缩速度。网格中的流体密度
$\; {\rho _{\text{f}}} $ 、黏度$\; {\mu _{\text{f}}} $ 可以分别表示为：$$ {\rho _{\rm{f}}} = \sum\limits_{i = 1}^n {{\alpha _i}{\rho _i}} $$ (3) $$ {\mu _{\rm{f}}} = \sum\limits_{i = 1}^n {{\alpha _i}{\mu _i}} $$ (4) 式中，i为网格数，ρi、µi分别为第i个网格的流体密度和黏度。
在CFD–DEM耦合模型中，流体的控制方程为Navier–Stokes方程，如式（5）～（6）所示：
$$ \frac{{\partial {\alpha _{\text{p}}}}}{{\partial t}}{\text{ + }}\nabla \cdot ({\alpha _{\text{p}}}{{\boldsymbol{u}}_{\text{f}}}){\text{ = }}0 $$ (5) $$ \begin{aligned}[b] & \frac{{\partial ({\alpha _{\text{p}}}{\rho _{\text{f}}}{{\boldsymbol{u}}_{\text{f}}})}}{{\partial t}}{\text{ + }}\nabla \cdot ({\alpha _{\text{p}}}{\rho _{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}{{\boldsymbol{u}}_{\text{f}}}) = - {\alpha _{\text{p}}}\nabla p - \\&\qquad {{\boldsymbol{R}}_{\text{pf}}} + {\alpha _{\text{p}}}\nabla \cdot {\boldsymbol{\tau}} + {\text{ }}{\alpha _{\text{p}}}{\rho _{\text{f}}}{\boldsymbol{g}} \end{aligned} $$ (6) 式（5）～（6）中：
${\boldsymbol{ \tau}} $ 为流体的应力张量；$ {\boldsymbol{g}} $ 为重力加速度；αp为孔隙率，以解释颗粒相占用的体积；$ {{\boldsymbol{R}}_{{\rm{p}}{{{\rm{f}}}}}} $ 为计算每个网格单元内基于拖拽力模型的动量交换项，代表了流体与颗粒的相互作用：$$ {{\boldsymbol{R}}_{{\rm{p}}{{{\rm{f}}}}}}{\text{ = }} {K_{{\rm{p}}{{{\rm{f}}}}}}[{{\boldsymbol{u}}_{\rm{f}}} - \left\langle {{{\boldsymbol{u}}_{\rm{p}}}} \right\rangle ] $$ (7) $$ {K_{{\rm{p}}{{{\rm{f}}}}}}{\text{ = }}\frac{{\displaystyle\sum\limits_{i = 1}^{{n_{\rm{p}}}} {{\boldsymbol{F}}_{{{i,{\rm{drag}}}}}^{}} }}{{{V_{{{{\rm{cell}}}}}} \cdot \left| {{{\boldsymbol{u}}_{{{\rm{f}}}}} - \left\langle {{{\boldsymbol{u}}_{\rm{p}}}} \right\rangle } \right|}} $$ (8) 式中，
$ \left\langle {{{\boldsymbol{u}}_{\rm{p}}}} \right\rangle $ 为基于流体网格的平均颗粒速度，$ {K_{{\rm{p}}{{{\rm{f}}}}}} $ 为动量交换系数，$ {\boldsymbol{F}}_{{{i,{\rm{drag}}}}}^{} $ 为单个颗粒受到的拖拽力，$ {V_{{{{\rm{cell}}}}}} $ 为网格体积，$ {n_{\rm{p}}} $ 为网格内的颗粒数量。1.2 离散单元法控制方程
在CFD–DEM耦合模型中，颗粒的运动状态由牛顿第二定律进行描述，包括平动及转动两部分。除颗粒之间的相互作用力和力矩外，平动及转动两部分控制方程中还分别添加了流体作用于颗粒上的力和力矩，方程如下：
$$ {m_i}\frac{{{\rm{d}}{{\boldsymbol{u}}_i}}}{{{\rm{d}}t}} = {m_i}{\boldsymbol{g}} + \sum\limits_{j = 1}^{{n_j}} {{\boldsymbol{F}}_{ij}^{}} + {\boldsymbol{F}}_i^{\rm{f}} $$ (9) $$ {{\boldsymbol{I}}_i}\frac{{{\rm{d}}{{\boldsymbol{\theta}} _i}}}{{{\rm{d}}t}}{\text{ = }}\sum\limits_{j = 1}^{{n_j}} {{\boldsymbol{M}}_{ij}^{} + {\boldsymbol{M}}_i^{\rm{f}}} $$ (10) 式（9）～（10）中，
$ {m_i} $ 为颗粒i的质量，$ {{\boldsymbol{u}}_i} $ 为颗粒i的平动速度，nj为与颗粒j接触的颗粒总数，$ {\boldsymbol{F}}_{ij}^{} $ 为颗粒i与颗粒j的相互作用力，$ {\boldsymbol{F}}_i^{\rm{f}} $ 为流体作用在颗粒i上的合力，$ {{\boldsymbol{I}}_i} $ 为转动惯量，$ {{\boldsymbol{\theta}} _i} $ 为颗粒i的旋转角速度，$ {\boldsymbol{M}}_{ij}^{} $ 为颗粒j作用于颗粒i的力矩，$ {\boldsymbol{M}}_i^{\rm{f}} $ 为流体作用于颗粒i的力矩，g为重力加速度。颗粒间相互作用力的计算采用基于Hertz接触理论的Hertz–Mindlin模型[21]。当两个颗粒的接触距离小于两者半径之和时，将产生相互作用力，法向包括弹簧力和阻尼力两项，切向包括剪切力和阻尼力两项。该模型计算公式如式（11）所示：
$$ {\;\;\;\;\;\;\;{\boldsymbol{F}}_{ij}^{} = \left( {{k_{\rm{n}}}{{\boldsymbol{\delta}} _{\rm{n}}} - {\gamma _{\rm{n}}}{\boldsymbol{v}}_{\rm{n}}^{}} \right) + \left( {{k_{\rm{t}}}{{\boldsymbol{\delta}} _{\rm{t}}} - {\gamma _{\rm{t}}}{\boldsymbol{v}}_{\rm{t}}^{}} \right) }$$ (11) 式中，
$ {k_{\rm{n}}} $ 、$ {k_{\text{t}}} $ 分别为法向、切向弹性系数，$ {\boldsymbol{\delta }}_{\rm{n}}^{} $ 为法向重叠量，$ {{\boldsymbol{\delta}} _{\rm{t}}} $ 为切向位移，$ {\gamma _{\rm{n}}} $ 、$ {\gamma _{\text{t}}} $ 分别为法向、切向阻尼系数，$ {\boldsymbol{v}}_{\rm{n}}^{} $ 为法向相对速度，$ {\boldsymbol{v}}_{\text{t}}^{} $ 为切向相对速度。1.3 流体–颗粒相互作用
在CFD–DEM耦合模型中，流体–颗粒间的相互作用通过在控制方程中交换相互作用力来考虑。本文主要考虑拖拽力和浮力的影响[21]，即：
$$ {{\boldsymbol{F}}_{\rm{f}}}{\text{ = }}{{\boldsymbol{f}}_{\rm{d}}} + {{\boldsymbol{f}}_{\rm{b}}} $$ (12) 式中，
$ {{\boldsymbol{F}}_{\rm{f}}} $ 为流体作用于颗粒上的合力，$ {{\boldsymbol{f}}_{\rm{d}}} $ 为流体与颗粒之间的拖拽力，$ {{\boldsymbol{f}}_{\rm{b}}} $ 为浮力。作用于颗粒上的浮力的表达式为：
$$ {{\boldsymbol{f}}_{\rm{b}}}{\text{ = }}\frac{4}{3}{\text{π}} {r^3}{\rho _{\rm{f}}}{\boldsymbol{g}} $$ (13) 式中，
$ r $ 为颗粒半径。作用于颗粒上的拖拽力的表达式采用di Felice关系式[22-24]，拖拽力的方向与颗粒运动方向相反，表达式为：
$$ {\;\;\;\;\;\;\;{{\boldsymbol{f}}_{\rm{d}}}{\text{ = }} \frac{1}{2}{C_{\rm{d}}}{\rho _{\rm{f}}}{\text{π}} {r^2}\left| {{{\boldsymbol{u}}_{\rm{f}}} - {{\boldsymbol{u}}_{\rm{p}}}} \right|({{\boldsymbol{u}}_{\rm{f}}} - {{\boldsymbol{u}}_{\rm{p}}}){\alpha _{\text{p}}^{ - \chi }}} $$ (14) 式中：Cd为拖拽力系数；
$ {{\boldsymbol{u}}_{\rm{p}}} $ 为颗粒速度；$ \chi $ 为孔隙率的一个经验修正系数，此系数使得该关系式能够适用于更大范围的雷诺数的情况。1.4 CFD–DEM耦合流程
在CFD–DEM耦合框架下，应用OpenFOAM程序求解流体运动方程，应用LIGGGHTS程序计算散粒体滑坡颗粒运动方程；将颗粒与流体的相互作用力作为耦合纽带，通过CFDEM耦合模块进行动量交互传递。
CFD–DEM耦合流程为：1）建立计算模型，给定初始流场、初始颗粒位置、初始颗粒速度；2）利用 DEM求解器计算作用于颗粒上的各个力，求解颗粒动量方程，更新每个颗粒的位置和速度；3）完成DEM循环后，计算CFD网格单元中的孔隙率与流体–颗粒相互作用力，并将其传递到CFD循环；4）求解流体动量方程和连续性控制方程，得到流体速度场与压力场；5）流场收敛后，计算作用于颗粒上的流体作用力，代回DEM求解器；6）CFD–DEM循环计算，直至达到计算结束时间。
2. 算例验证
2.1 颗粒堆积体坍塌–涌浪试验
为验证CFD–DEM流固耦合模型模拟滑坡–涌浪灾害的有效性，对Robbe–Saule等[25]所做的室内水槽试验进行数值模拟。该试验过程包含了颗粒堆积体失稳、涌浪形成及传播等类似滑坡入水的关键过程，研究了颗粒堆积体的高宽比和体积、颗粒的直径和密度等因素对涌浪高度的影响情况。因此，被广泛用作标准检验模型。
试验装置左侧为颗粒堆积区，高宽比为H0/L0；右侧为一个长2.00 m、高0.30 m、宽0.15 m的玻璃水槽，其中充满深度为0.05 m的水。通过提拉挡板释放颗粒，通过改变颗粒堆积体的高度和宽度控制体积变量V。数值模拟使用的计算参数见表1。
表 1 颗粒堆积体坍塌–涌浪试验计算参数Table 1 Calculation parameters of the particle accumulation collapse–surge experiment模块 计算参数 取值 DEM模块 颗粒半径/m 0.002 5 颗粒密度/(kg·m–3) 2 500 泊松比 0.25 杨氏模量/Pa 2.5×107 碰撞恢复系数 0.5 静摩擦系数 0.5 滚动摩擦系数 0.01 CFD模块 水体密度/(kg·m–3) 1 000 运动黏度/(m2·s–1) 1.0×10–6 表面张力系数/(N·m–1) 0.07 耦合模块 耦合步数 10 时间步长/s 0.000 1 计算时长/s 1 2.2 对比验证结果
图1为在H0/L0=2.5，V=10.2 dm3工况下，颗粒堆积体坍塌演进过程中不同时刻的水槽试验结果（图1左）与数值模拟结果（图1右）。对比图1中颗粒堆积体坍塌过程及水面涌浪运动状态可以清晰看出，在整个演进过程中，数值模拟结果与水槽试验结果基本保持一致。当颗粒撞击水面时会形成初始涌浪，随着颗粒体大量侵入水体，涌浪高度随之增加；而后，涌浪在传播过程中逐渐变陡，波峰在重力影响下溅落到水体表面。除此之外，颗粒堆积体坍塌到水中产生涌浪的整个过程，涉及到颗粒、水体及空气的复杂相互作用。由图1试验结果可以观察到，水体在短时间内并不能完全充满颗粒间隙，在颗粒堆积体坍塌过程中将形成一个明显的凹形干湿边界；对比数值模拟结果可知，CFD–DEM流固耦合模型能够较为准确地捕捉这一现象，初步验证了该模型解决流固耦合问题的可靠性。
图 1 颗粒堆积体坍塌演进过程对比Fig. 1 Comparison of particle accumulation collapse evolution process
下载:
全尺寸图片
通过提取水槽试验及数值模拟中各时刻涌浪高度，并进行对比分析，结果如图2所示。从图2中的涌浪峰值高度看，试验中测量的最大涌浪高度为8.3 cm，模拟计算得到的最大涌浪高度为8.1 cm，结果较为一致。从图2中涌浪高度演进过程来看，模拟结果稍微超前于试验结果，但计算模拟的整体趋势及高差范围与试验结果较为吻合。分析造成误差的原因有两点：一是，实际试验中需提拉挡板，造成时间误差，导致涌浪延后产生，在提升过程中对颗粒堆积体也会造成扰动，一定程度上影响涌浪产生过程；二是，试验中颗粒大多是不规则的，为简化模拟采用粒径单一的球形颗粒代替，会低估颗粒间的相互作用。
图 2 数值模拟与物理试验得到的涌浪高度对比Fig. 2 Comparison of surge height between numerical and experimental results下载:
全尺寸图片
初始颗粒堆积体的体积对涌浪的形成起着重要作用，因此，Robbe–Saule等[25]进行了一系列不同颗粒堆积体体积工况下的试验。根据对应的试验参数，本文进行颗粒堆积体体积分别为3.71、5.69、7.42、10.15 dm3的4组数值模拟，以进一步验证CFD–DEM流固耦合模型在不同体积工况下的适用性。各组工况下试验与模拟结果的最大涌浪高度如图3所示。由图3可知：在一定范围内，最大涌浪高度随初始颗粒堆积体体积的增大而增大；在各组颗粒堆积体体积工况下，模拟与试验结果的最大涌浪高度拟合度较高，误差均在8%以内。
图 3 不同颗粒体积工况下最大涌浪高度对比Fig. 3 Comparison of maximum surge height under different particle volume conditions下载:
全尺寸图片
通过对颗粒坍塌运动过程、涌浪高度演化过程等进行对比，发现颗粒运动及涌浪传播模拟结果与试验结果吻合很好，可证明本文的流固耦合模型能有效还原颗粒运动状态及涌浪的波动变化。
3. 色玉滑坡–涌浪灾害演进过程分析
3.1 色玉滑坡–涌浪灾害概况
基于上述CFD–DEM流固耦合模型，对四川省猴子岩水库色玉堆积体前缘滑坡进行数值模拟研究，分析滑坡失稳、涌浪形成和传播的灾害演进全过程。
猴子岩水电站位于四川省甘孜藏族自治州康定县境内大渡河干流上；色玉堆积体位于电站库区右岸，距下游大坝约4.3 km。色玉堆积体顺河长约650 m，横河宽约350 m，前、后缘高程分别为1 820、2 100 m，高差约为280 m。高程约1 930～2 100 m处为较宽缓的台状地形，台地前缘至坡脚地形坡度约为50～70°。堆积台地主要由冰碛堆积物（glQ3）组成，台地后缘覆盖少量崩坡积堆积物（col+dlQ4）。堆积体厚度较大，一般为83～100 m；下伏基岩主要为泥盆系中上统河心组（D2–3h）灰、灰白色中厚层白云岩、白云质灰岩夹灰岩。
2017年10月24日，色玉堆积体前缘发生滑坡，此时水库蓄水位为1 832 m。滑坡失稳后的地貌特征如图4所示。由图4可知，失稳方量约80×104 m3，失稳高度最大达160 m，失稳后在下游3.7 km处的电站进水口形成1.7 m高的涌浪。截止目前，该堆积体又多次发生垮塌，并形成多处凹槽地形。由于该堆积体距离大坝相对较近，其失稳涌浪对工程和周边人类活动安全的不利影响需重点关注和研究。
图 4 色玉滑坡全貌Fig. 4 Panorama of the Seyu landslide下载:
全尺寸图片
3.2 计算模型及参数设置
通过对色玉堆积体失稳前后地形数据进行差值分析，确定模拟的物源。计算使用了2万多个半径为2 m的球形颗粒，参考Hu等[26]所做的色玉滑坡附近区域土工试验测试结果，设置模型中颗粒间的泊松比和杨氏模量等参数。色玉滑坡数值计算模型见图5，计算域网格大小为40 m×40 m×40 m。为减小涌浪高度计算误差，将库水位上下20 m范围内的网格进行加密，网格大小为10 m×10 m×10 m。计算参数选取见表2。
图 5 色玉滑坡3维模型Fig. 5 3D model of Seyu landslide下载:
全尺寸图片
表 2 色玉滑坡计算参数Table 2 Calculation parameters of Seyu landslide模块 计算参数 取值 DEM模块 颗粒半径/m 2.0 颗粒密度/(kg·m–3) 2 150 泊松比 0.35 杨氏模量/Pa 2.5×107 碰撞恢复系数 0.8 静摩擦系数 0.5 滚动摩擦系数 0.01 CFD模块 水体密度/(kg·m–3) 1 000 运动黏度/(m2·s–1) 1.0×10–6 表面张力系数/(N·m–1) 0.07 耦合模块 耦合步数 10 时间步长/s 0.01 计算时长/s 200 3.3 模拟结果
为验证数值模拟的准确性，将涌浪计算结果与实测数据进行对比。图5中，P点（225，–2 700）所示位置为色玉堆积体下游3.7 km的电站进水口处，堆积体失稳后，在该处实测到了高1.7 m的涌浪。提取P点的涌浪高度计算结果如图6所示。
图 6 P点涌浪高度变化Fig. 6 Changes of surge height at point P下载:
全尺寸图片
由图6可知：在滑坡发生142 s后，首浪传播到P点，高度为2.12 m，这与实测的数据较为接近。当涌浪还没有传播到P点时，水面在±0.5 m范围内波动，这是由于流体初始条件设置及部分网格不规则造成的误差。通过对该处涌浪计算结果与实测数据进行对比，可以证明色玉滑坡–涌浪灾害的数值模拟结果较为准确。
色玉滑坡运动及涌浪传播动态演进过程模拟结果如图7所示，显示了0～30 s时间内滑坡及涌浪形态。由图7可知：从滑坡体失稳运动至静止堆积的持续时间约为20 s，从滑坡入水处涌浪形成到涌浪传播至对岸的持续时间约为10 s。t=5 s左右时，在滑坡颗粒冲击作用下，涌浪开始逐渐形成；t=10 s，大部分滑坡颗粒滑下，涌浪高度进一步增加，沿滑坡发生点呈圆弧状向外扩散；t=15 s，滑坡颗粒全部进入水库，近岸水体沿滑坡颗粒运动方向下陷，涌浪在对岸已经开始爬升；t=20 s，滑坡颗粒不再运动，堆积于水库底部，此刻对岸涌浪爬升到最大高度；t=20～30 s，涌浪开始衰退，并继续向上下游传播。
图 7 色玉滑坡运动及涌浪传播动态演进过程Fig. 7 Dynamic evolution process of Seyu landslide movement and surge propagation下载:
全尺寸图片
图8展示了滑坡颗粒速度和涌浪速度在滑坡从失稳运动到静止堆积这段时间内的变化情况。图9展示了所有滑坡颗粒的最大速度及平均速度随时间的变化曲线。
图 8 滑坡颗粒及涌浪速度Fig. 8 Velocity of the landslide particles and surge下载:
全尺寸图片
图 9 滑坡颗粒速度随时间变化曲线Fig. 9 Variation curves of landslide particles velocity with time下载:
全尺寸图片
由图8、9可知：随着滑坡运动不断发展，重力势能转换为动能，滑坡颗粒速度开始处于快速增加阶段，平均速度在6.0 s时达到最大值，为16.12 m/s；滑坡中部运动速度明显高于两侧运动速度，滑坡前部运动速度高于滑坡后部运动速度。之后，由于滑坡前缘颗粒持续滑入水库底部，受到地形影响，颗粒在水库底部逐渐减速并开始堆积，滑坡体的平均速度也开始减小。直至滑坡发生约20 s后，平均速度减到0.04 m/s，可以认为此时滑坡完全停止运动。虽然滑坡整体运动速度在6 s后减小，但仍有部分颗粒在激烈碰撞下达到较高的运动速度，单个颗粒的最大速度达到44.29 m/s，发生在10.0 s左右。
采用横向剖面进一步分析涌浪沿滑坡运动方向的传播过程及涌浪的爬升高度，剖面线位置如图10所示。为更加精确地显示涌浪高度变化，在剖面上选取3个监测点对各时刻的涌浪高度进行提取分析，监测点分别为滑坡入水处A点、河道中部B点及滑坡对岸处C点，3点坐标分别为（472，1 018）、（561，946）、（658，868）。
图 10 剖面线及监测点位置Fig. 10 Location of profile lines and monitoring points下载:
全尺寸图片
图11为剖面位置涌浪演进过程模拟分析结果，可清楚观察到滑坡发生后短时间内产生巨大涌浪的全过程。
图 11 涌浪演进过程剖面Fig. 11 Profile of surge evolution process下载:
全尺寸图片
结合图11可以看出，涌浪运动过程主要分为3个阶段，具体如下：
1）涌浪产生及传播阶段：滑坡失稳后在重力作用下冲击水库，受滑坡前缘颗粒挤压，水面沿滑坡方向向下凹陷，涌浪高度快速增加并向对岸传播。
2）涌浪爬升阶段：色玉滑坡方量较大，且由于河谷狭窄，涌浪在10 s后即传播到对岸，在惯性力作用下开始沿坡面向上爬升，并在20 s左右爬升到最大高度27.32 m；在这一阶段，滑坡入水处水体产生的凹陷开始收缩。
3）涌浪回流阶段：在20 s之后，由于滑坡的运动趋于静止，冲向对岸的涌浪在重力的作用下，开始回流。
剖面上监测点的涌浪高度变化如图12所示。
图 12 沿滑坡运动方向监测点涌浪高度变化Fig. 12 Changes of surge height at monitoring points along the landslide movement direction下载:
全尺寸图片
由图12可知：A点显示了滑坡发生后产生的初始涌浪高度，在3.5 s时达到6.84 m，随后近岸水体在大量滑坡颗粒冲击下迅速下陷；涌浪波峰在10.5 s时传播至B点，此时河道中部的最大涌浪高度为16.62 m；C点为滑坡对岸附近，可显示涌浪传播到对岸边时的高度，该点涌浪高度在15.5 s时达到最大值19.04 m。滑坡产生的涌浪经过A、B、C点之后，开始在对岸沿坡面爬升，最大爬升高度达到27.32 m。由此可见，在山区狭窄河谷地形条件下，水库岸坡失稳将产生巨大涌浪，尤其在对岸爬高最大。
为分析滑坡涌浪沿河道方向的传播过程，在滑坡区与电站进水口P点之间，等距离选取了D（525，200）、E（725，–525）、F（700，–1 272）、G（425，–1 975）、P（225，–2 700）5个监测点进行涌浪高度监测，具体位置如图5所示。沿河道方向各监测点涌浪高度变化如图13所示。由图13可知：相较于滑坡涌浪近场区域，远场区域涌浪高度明显降低。当滑坡体高速冲击水体时，能量主要沿滑动方向传递，造成沿滑动方向上的涌浪高度、传播速度和能量大，沿河道方向较小。色玉滑坡失稳入水后，涌浪传播到各监测点的时间分别为22、46、80、110、142 s。滑坡涌浪在河流弯道处EF段的传播时间最长，在直道段传播时间稍短。以上结果表明滑坡涌浪在沿河道传播过程中，其传播特征受到库区地形条件的影响较大。各监测点最大涌浪高度分别为4.26、3.75、3.15、3.09、2.12 m，均出现在第1个波峰，后续波峰峰值衰减较大。沿程涌浪是一个逐渐衰减过程，最大涌浪高度随着传播距离的增加而减小。
图 13 沿河道方向监测点涌浪高度变化Fig. 13 Changes of surge height at monitoring points along the river下载:
全尺寸图片
4. 结 论
1）介绍了一种CFD–DEM耦合算法，流体的控制方程为连续性方程及动量方程，颗粒的运动由牛顿第二定律控制，通过基于拖拽力模型的动量交换项实现流体与固体之间的耦合。
2）基于颗粒堆积体坍塌试验，进行了不同颗粒体积工况下的模型验证，模拟结果与试验结果吻合较好，证明了耦合模型能准确模拟颗粒与流体的相互作用过程。
3）将耦合模型应用于实际发生的大规模滑坡入水–涌浪灾害链的工程四川省猴子岩水库色玉滑坡–涌浪的模拟。结果表明：滑坡体从失稳运动至静止堆积的持续时间约20 s，颗粒平均速度最大达到16.12 m/s；滑坡引起的涌浪在约10 s后传播到对岸，之后开始沿坡面向上爬升，最大爬升高度达到27.32 m。
CFD–DEM耦合模型不仅对室内坍塌–涌浪试验的模拟具有高效性和准确性，也可以应用于山区河谷大规模滑坡入水–涌浪灾害链这一复杂流固耦合问题的模拟，实现滑坡失稳运动及涌浪形成传播全过程分析，为山区河谷防灾减灾提供新的技术手段。
-
图 1 颗粒堆积体坍塌演进过程对比
Fig. 1 Comparison of particle accumulation collapse evolution process
下载:
全尺寸图片
图 2 数值模拟与物理试验得到的涌浪高度对比
Fig. 2 Comparison of surge height between numerical and experimental results
下载:
全尺寸图片
图 3 不同颗粒体积工况下最大涌浪高度对比
Fig. 3 Comparison of maximum surge height under different particle volume conditions
下载:
全尺寸图片
图 4 色玉滑坡全貌
Fig. 4 Panorama of the Seyu landslide
下载:
全尺寸图片
图 5 色玉滑坡3维模型
Fig. 5 3D model of Seyu landslide
下载:
全尺寸图片
图 6 P点涌浪高度变化
Fig. 6 Changes of surge height at point P
下载:
全尺寸图片
图 7 色玉滑坡运动及涌浪传播动态演进过程
Fig. 7 Dynamic evolution process of Seyu landslide movement and surge propagation
下载:
全尺寸图片
图 8 滑坡颗粒及涌浪速度
Fig. 8 Velocity of the landslide particles and surge
下载:
全尺寸图片
图 9 滑坡颗粒速度随时间变化曲线
Fig. 9 Variation curves of landslide particles velocity with time
下载:
全尺寸图片
图 10 剖面线及监测点位置
Fig. 10 Location of profile lines and monitoring points
下载:
全尺寸图片
图 11 涌浪演进过程剖面
Fig. 11 Profile of surge evolution process
下载:
全尺寸图片
图 12 沿滑坡运动方向监测点涌浪高度变化
Fig. 12 Changes of surge height at monitoring points along the landslide movement direction
下载:
全尺寸图片
图 13 沿河道方向监测点涌浪高度变化
Fig. 13 Changes of surge height at monitoring points along the river
下载:
全尺寸图片
表 1 颗粒堆积体坍塌–涌浪试验计算参数
Table 1 Calculation parameters of the particle accumulation collapse–surge experiment
模块 计算参数 取值 DEM模块 颗粒半径/m 0.002 5 颗粒密度/(kg·m–3) 2 500 泊松比 0.25 杨氏模量/Pa 2.5×107 碰撞恢复系数 0.5 静摩擦系数 0.5 滚动摩擦系数 0.01 CFD模块 水体密度/(kg·m–3) 1 000 运动黏度/(m2·s–1) 1.0×10–6 表面张力系数/(N·m–1) 0.07 耦合模块 耦合步数 10 时间步长/s 0.000 1 计算时长/s 1 表 2 色玉滑坡计算参数
Table 2 Calculation parameters of Seyu landslide
模块 计算参数 取值 DEM模块 颗粒半径/m 2.0 颗粒密度/(kg·m–3) 2 150 泊松比 0.35 杨氏模量/Pa 2.5×107 碰撞恢复系数 0.8 静摩擦系数 0.5 滚动摩擦系数 0.01 CFD模块 水体密度/(kg·m–3) 1 000 运动黏度/(m2·s–1) 1.0×10–6 表面张力系数/(N·m–1) 0.07 耦合模块 耦合步数 10 时间步长/s 0.01 计算时长/s 200 -
[1] 黄波林,殷跃平.水库区滑坡涌浪风险评估技术研究[J].岩石力学与工程学报,2018,37(3):621–629. doi: 10.13722/j.cnki.jrme.2017.1047 Huang Bolin,Yin Yueping.Risk assessment research on impulse wave generated by landslide in reservoir[J].Chinese Journal of Rock Mechanics and Engineering,2018,37(3):621–629 doi: 10.13722/j.cnki.jrme.2017.1047 [2] 王芳,殷坤龙,桂蕾,等.单体库岸滑坡及其次生涌浪灾害风险分析[J].地球科学,2018,43(3):899–909. Wang Fang,Yin Kunlong,Gui Lei,et al.Risk analysis on individual reservoir bank landslide and its generated wave[J].Earth Science,2018,43(3):899–909 [3] Fritz H M,Mohammed F,Yoo J.Lituya bay landslide impact generated mega-tsunami 50th anniversary[J].Pure and Applied Geophysics,2009,166(1/2):153–175. doi: 10.1007/s00024-008-0435-4 [4] Jian Wenxing,Xu Qiang,Yang Hufeng,et al.Mechanism and failure process of Qianjiangping landslide in the Three Gorges Reservoir,China[J].Environmental Earth Sciences,2014,72(8):2999–3013. doi: 10.1007/s12665-014-3205-x [5] 蒋权,陈希良,肖江剑,等.云南黄坪库区滑坡运动及其失稳模式的离散元模拟[J].中国地质灾害与防治学报,2018,29(3):53–59. doi: 10.16031/j.cnki.issn.1003-8035.2018.03.08 Jiang Quan,Chen Xiliang,Xiao Jiangjian,et al.Discrete element numerical simulation and analysis of Yunnan Huangping Reservoir areas landslide and its failure mode[J].The Chinese Journal of Geological Hazard and Control,2018,29(3):53–59 doi: 10.16031/j.cnki.issn.1003-8035.2018.03.08 [6] Zhou Jiawen,Xu Fugang,Yang Xingguo,et al.Comprehensive analyses of the initiation and landslide-generated wave processes of the 24 June 2015 Hongyanzi landslide at the Three Gorges Reservoir,China[J].Landslides,2016,13(3):589–601. doi: 10.1007/s10346-016-0704-8 [7] Xiao Lili,Wang Jiajia,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. doi: 10.1007/s10346-018-1057-2 [8] 徐文杰.滑坡涌浪流–固耦合分析方法与应用[J].岩石力学与工程学报,2020,39(7):1420–1433. doi: 10.13722/j.cnki.jrme.2020.0061 Xu Wenjie.Fluid-solid coupling method of landslide tsunamis and its application[J].Chinese Journal of Rock Mechanics and Engineering,2020,39(7):1420–1433 doi: 10.13722/j.cnki.jrme.2020.0061 [9] 刘威,杨宗佶,游勇,等.滑坡涌浪诱发冰湖溃决灾害链过程分析与模拟[J].工程科学与技术,2022,54(2):41–48. doi: 10.15961/j.jsuese.202100577 Liu Wei,Yang Zongji,You Yong,et al.Dynamic analysis and numerical simulation of disaster chain from landslide,surge wave to glacier lake outburst[J].Advanced Engineering Sciences,2022,54(2):41–48 doi: 10.15961/j.jsuese.202100577 [10] Katz O,Morgan J K,Aharonov E,et al.Controls on the size and geometry of landslides:Insights from discrete element numerical simulations[J].Geomorphology,2014,220:104–113. doi: 10.1016/j.geomorph.2014.05.021 [11] Mead S R,Cleary P W.Validation of DEM prediction for granular avalanches on irregular terrain[J].Journal of Geophysical Research(Earth Surface),2015,120(9):1724–1742. doi: 10.1002/2014JF003331 [12] 徐寅,陈胜宏.基于离散单元法的滑坡堆积及其涌浪计算[J].岩土力学,2012,33(9):2850–2856. doi: 10.16285/j.rsm.2012.09.020 Xu Yin,Chen Shenghong.Calculation of heap shape of landslide and its surge based on discrete element method[J].Rock and Soil Mechanics,2012,33(9):2850–2856 doi: 10.16285/j.rsm.2012.09.020 [13] Abadie S,Morichon D,Grilli S,et al.Numerical simulation of waves generated by landslides using a multiple-fluid Navier–Stokes model[J].Coastal Engineering,2010,57(9):779–794. doi: 10.1016/j.coastaleng.2010.03.003 [14] Li Xingyue,Zhao Jidong.A unified CFD–DEM approach for modeling of debris flow impacts on flexible barriers[J].International Journal for Numerical and Analytical Methods in Geomechanics,2018,42(14):1643–1670. doi: 10.1002/nag.2806 [15] Park K M.CFD–DEM based numerical simulation of liquid-gas-particle mixture flow in dam break[J].Communications in Nonlinear Science and Numerical Simulation,2018,59:105–121. doi: 10.1016/j.cnsns.2017.11.010 [16] Shi Zhenming,Zheng Hongchao,Yu Songbo,et al.Application of CFD–DEM to investigate seepage characteristics of landslide dam materials[J].Computers and Geotechnics,2018,101:23–33. doi: 10.1016/j.compgeo.2018.04.020 [17] Zhao T,Utili S,Crosta G B.Rockslide and impulse wave modelling in the vajont reservoir by DEM–CFD analyses[J].Rock Mechanics and Rock Engineering,2016,49(6):2437–2456. doi: 10.1007/s00603-015-0731-0 [18] Vångö M,Pirker S,Lichtenegger T.Unresolved CFD–DEM modeling of multiphase flow in densely packed particle beds[J].Applied Mathematical Modelling,2018,56:501–516. doi: 10.1016/j.apm.2017.12.008 [19] Jasak H.OpenFOAM:Open source CFD in research and industry[J].International Journal of Naval Architecture and Ocean Engineering,2009,1(2):89–94. doi: 10.2478/IJNAOE-2013-0011 [20] Berger R.Hybrid parallelization of the LIGGGHTS open-source DEM code[J].Powder Technology,2015,278:234–247. doi: 10.1016/j.powtec.2015.03.019 [21] Potyondy D O.A bonded-particle model for rock[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(8):1329–1364. doi: 10.1016/j.ijrmms.2004.09.011 [22] 王胤,艾军,杨庆.考虑粒间滚动阻力的CFD–DEM流–固耦合数值模拟方法[J].岩土力学,2017,38(6):1771–1780. doi: 10.16285/j.rsm.2017.06.027 Wang Yin,Ai Jun,Yang Qing.A CFD–DEM coupled method incorporating soil inter-particle rolling resistance[J].Rock and Soil Mechanics,2017,38(6):1771–1780 doi: 10.16285/j.rsm.2017.06.027 [23] Felice R D.The voidage function for fluid-particle interaction systems[J].International Journal of Multiphase Flow,1994,20(1):153–159. doi: 10.1016/0301-9322(94)90011-6 [24] Zhao Jidong.Coupled CFD–DEM simulation of fluid-particle interaction in geomechanics[J].Powder Technology,2013,239:248–258. doi: 10.1016/j.powtec.2013.02.003 [25] Robbe–Saule M,Morize C,Henaff R,et al.Experimental investigation of tsunami waves generated by granular collapse into water[J].Journal of Fluid Mechanics,2021,907:A11. doi: 10.1017/jfm.2020.807 [26] Hu Yuxiang,Zhu Yongguo,Li Haibo,et al.Numerical estimation of landslide-generated waves at Kaiding Slopes,Houziyan Reservoir,China,using a coupled DEM–SPH method[J].Landslides,2021,18(10):3435–3448. doi: 10.1007/s10346-021-01718-1
