Coupled DEM–CFD Method for Landslide–river Blockage–impulse Wave Disaster Chain Simulation and its Application
-
摘要: 滑坡–堵江–涌浪灾害链的演化过程涉及复杂的滑坡–河流相互作用。为了真实再现这一演进过程，本文提出了一种扩展的离散元法（DEM）与计算流体动力学（CFD）耦合数值方法，通过在局部平均DEM–CFD耦合方法中引入流体体积法（VOF）建立了一套能够描述河流自由水面演化的扩展DEM–CFD耦合数值框架，实现了对滑坡涌浪形成及传播过程的追踪；并提出了一种新的局部孔隙率计算方法即虚拟球模型，来克服网格颗粒临界尺寸比的限制，从而实现复杂地形建模要求；进一步通过多个典型算例验证其准确性和有效性。在此基础上，开展了滑坡–堵江–涌浪灾害链演化模拟，并应用至白格滑坡坝形成过程复演，深入分析了滑坡–河流动态演化过程及耦合作用机理。结果表明：所提出的方法能够准确模拟滑坡–河流相互作用；滑坡–河流–河谷地形的耦合作用驱动了灾害演化，滑坡体驱动了涌浪的传播，而河流则增加了滑坡体的动能耗散，且对其沉积过程在河流流向上产生相反影响；模拟得到的白格滑坡运移路径、沉积形态及涌浪侵蚀面积与实地调查结果吻合度较高，成功再现了2018年10月11日白格滑坡堵江事件。本文所提出的扩展DEM–CFD耦合数值方法为深入研究滑坡–堵江–涌浪灾害链演化机制提供了有力的数值工具，对于灾害预测及防灾减灾策略的制定具有一定的参考价值。Abstract: The intrusion of landslides into a river can block the river and induce surge waves. The evolution of the landslide–river blockage–surge wave disaster chain involves complex landslide–river interactions. To reproduce it realistically, an extended coupled discrete element method (DEM) and computational fluid dynamics (CFD) numerical method is proposed in this research: the VOF model is introduced to track the evolution of river free surface; a virtual sphere model is proposed to overcome the limitation of the critical size ratio of the mesh and particle; and the model is verified by several typical cases. Based on this, the simulation of the landslide–river blockage–surge wave disaster chain was carried out to reproduce the formation process of the Baige landslide-induced river damming, and the dynamic evolution process and interaction mechanism of landslide–river dynamics are investigated comprehensively. The results show that the proposed method can accurately simulate the landslide–river interaction during the evolution of landslide–river blockage–surge wave disaster chain; the landslide mass drives the propagation of surge waves, while the river significantly increases the kinetic energy dissipation of the landslide mass and has opposites effect on the deposition of the landslide along the flow directions; the simulated migration path, deposit morphology and wave erosion area of the Baige landslide are in good agreement with the field survey results, and the method successfully reproduces the October 11th, 2018 Baige landslide. The extended coupled DEM–CFD method proposed in this research provides a robust numerical tool for the understanding of the river blockage and impulse wave evolution mechanism, as well as the prediction of the disaster evolution process, which is of reference value for the of disaster emergency response and disaster mitigation strategies.
-
中国西南地区，尤其是青藏高原地区，沟壑林立且河流发育，一直是滑坡类灾害的高发区[1-2]。滑坡体阻塞河流所形成的堰塞坝可诱发如上游洪涝、溃坝洪水及泥石流等一系列次生灾害[3-4]。同时，滑坡体侵入河流造成的涌浪会导致灾害影响范围显著扩大[5-6]。滑坡–堵江–涌浪灾害链具有致灾机理复杂、致灾范围广及破坏性强等显著特点，可对沿岸居民生命财产安全及流域内水利设施安全造成严重威胁。
滑坡–堵江–涌浪灾害链演化过程涉及复杂的滑坡–河流多相动力演化及耦合作用，从数值角度对其演化过程进行模拟是一个极大的挑战[7]。一些研究基于连续力学理论的双层深度平均方法将滑坡体及河流用两相流体建模[7-9]。然而，上述研究假定滑坡体为连续介质流体，对固–液耦合作用的描述过于简化，且忽略了滑坡体具有孔隙的非连续特性[10]。相较而言，基于非连续介质理论的离散单元法（discrete element method，DEM）通过颗粒对土体建模，在处理土体大变形问题中具有突出优势[11]。DEM与计算流体动力学（computational fluid dynamics，CFD）的耦合方法能够有效处理土–水耦合作用问题，已成功应用于水下滑坡模拟[12-13]。然而，传统的固–液两相DEM–CFD耦合方法并不足以描述滑坡涌浪的形成及传播过程。此外，局部平均的DEM–CFD耦合方法受限于网格颗粒临界尺寸比，难以满足滑坡堵江及涌浪过程模拟中滑坡体粗颗粒材料的建模及高分辨率真实地形网格的要求，这些都限制了DEM–CFD耦合数值方法在滑坡–堵江–涌浪灾害链模拟中的应用。
本文旨在通过数值模拟方法研究滑坡–堵江–涌浪灾害链演化过程。首先，通过在局部平均DEM–CFD耦合方法中引入流体体积法（volume of fluid，VOF）建立了一套能够描述河流自由水面演化的扩展DEM–CFD耦合数值框架，实现了对滑坡涌浪形成及传播过程的追踪。然后，提出了一种新的局部孔隙率计算方法，即虚拟球模型，来克服网格颗粒临界尺寸比的限制。在对该方法的准确性和有效性进行验证后，在模型尺度上模拟了滑坡堵江及涌浪演化过程，并分析其滑坡–河流耦合作用机制。最后，将其应用至2018年10月11日金沙江白格滑坡堵江事件的数值反演。
1. 扩展DEM–CFD流固耦合数值方法
1.1 控制方程
本研究所提出的DEM–CFD耦合数值方法中，滑坡体通过DEM建模，并通过求解牛顿第二定律来模拟滑坡体的运移过程[14]：
$$ {m_{{i}}}\frac{{{\rm{d}}{{\boldsymbol{v}}_{{i}}}}}{{{\rm{d}}t}} = {m_{{i}}}{\boldsymbol{g}} + \sum\limits_{{{j}} = 1}^{{n_{{i}}}} {{{\boldsymbol{F}}_{{{ij}}}}} + {{\boldsymbol{F}}_{{{\rm{f}}}}} $$ (1) $$ {I_{{i}}} = \frac{{{\rm{d}}{\text{ω}_{{i}}}}}{{{\rm{d}}t}} = \sum\limits_{{{j}} = 1}^{{{{n}}_{{i}}}} {{{\boldsymbol{M}}_{{{ij}}}}} $$ (2) 式（1）～（2）中，mi、vi和
${\text{ω}_{{i}}} $ 分别为颗粒质量、速度及角速度，g为重力加速度，Fij和Mij分别为颗粒的接触作用力及扭矩，Ff 为流体对颗粒的耦合作用力。对于流体相，通过CFD求解颗粒相体积分数修正后的局部平均Navier–Stokes方程[15]：
$$ \frac{{\partial \varepsilon }}{{\partial t}} + \nabla \cdot (\varepsilon {\boldsymbol{u}}) = 0 $$ (3) $$\begin{aligned}[b]& \frac{{\partial ({\rho _{\text{f}}}\varepsilon {\boldsymbol{u}})}}{{\partial t}} + \nabla \cdot ({\rho _{\text{f}}}\varepsilon {\boldsymbol{uu}}) = \\ &\qquad - \varepsilon \nabla p + \varepsilon \nabla \cdot ({\mu _{\text{f}}}(\nabla {\boldsymbol{u}} + \nabla {{\boldsymbol{u}}^{\text{T}}})) + \varepsilon {\rho _{\text{f}}}g + {f_{{\text{pf}}}} \end{aligned}$$ (4) 式（3）～（4）中，
$ \varepsilon $ 为局部孔隙率，$ {\boldsymbol{u}} $ 、$ \;{\rho _{\text{f}}} $ 、$\; {\mu _{\text{f}}} $ 分别为流体速度、密度及黏度，p为流体压力，$ {f_{{\text{pf}}}} $ 为表征颗粒对流体阻力的动量源项，式（4）右侧其余3项则依次为压力梯度项、应力项及重力项。1.2 VOF模型
在本研究中，将VOF法引入传统的局部平均DEM–CFD耦合方法中，从而允许对多相流体介质建模，涌浪的形成及演化过程通过对多相流体交界面的追踪来实现。本研究的扩展DEM–CFD耦合方法中，每个流体单元被分为水相和空气相，并用
$ {\alpha _{\text{w}}} $ 和$ {\alpha _{\text{a}}} $ 分别表示相应相的体积分数。其中，水相的体积分数满足$0 \le {\alpha _{\text{w}}} \le 1$ ，则空气相的体积分数可表示为[16]：$$ {\alpha _{\text{a}}} = 1 - {\alpha _{\text{w}}} $$ (5) 因此，
$ 0 \lt {\alpha _{\text{w}}} \lt 1 $ 被定义为自由液面的存在。式（4）中的流体密度$ \;{\rho _{\text{f}}} $ 及黏度$\; {\mu _{\text{f}}} $ ，参考两相流体性质，通过线性插值计算。通过求解基于流体速度的输运方程，自由液面运动可计算如下[17]：$${\;\;\;\; \frac{{\partial (\varepsilon {\alpha _{\text{w}}})}}{{\partial t}} + \nabla \cdot (\varepsilon {\alpha _{\text{w}}}{\boldsymbol{u}}) + \nabla \cdot (\varepsilon {\alpha _{\text{w}}}(1 - {\alpha _{\text{w}}}){{\boldsymbol{u}}_{\text{r}}}) = 0} $$ (6) 式中，
${\boldsymbol{u}}_{\rm{r}} $ 为两相介质相对速度，$ {{\boldsymbol{u}}_{\text{r}}} = {{\boldsymbol{u}}_1} - {{\boldsymbol{u}}_2} $ 。值得注意的是，方程中已经考虑颗粒相的存在。左侧第1项中$ \varepsilon $ 的时间导数考虑到了体积分数变化对液相体积分数的影响。式（6）的前两项与式（4）具有相同的平流形式，第3项为人工项，用于减少数值耗散，以达到使自由液面更为清晰的目的。1.3 耦合作用力计算
本研究中，流体与颗粒之间的耦合作用力
$ {F_{\text{f}}} $ 包括浮力$ {F_{\text{b}}} $ 和拖曳力$ {F_{\text{d}}} $ ，其中，浮力$ {F_{\text{b}}} $ 基于压力梯度求解，可表示为：$$ {F_{\text{b}}} = - {V_{\text{p}}}\nabla p $$ (7) 式中，
$ {V_{\text{p}}} $ 为颗粒体积。拖曳力是由流体域颗粒间的相对运动引起的，作用于颗粒质心，其方向与相对运动方向相反。通过Di Felice模型来计算拖曳力[18]：$$ {\;\;\;\;\;{F_{\text{d}}} = \frac{1}{2}{C_{\text{d}}}{\rho _{\text{f}}}\frac{{{\text{π}} {d^2}}}{4}\left| {{\boldsymbol{u}} - {\boldsymbol{v}}} \right|({\boldsymbol{u}} - {\boldsymbol{v}}){\varepsilon ^{ - \chi + 1}} }$$ (8) 式中：
${C_{\rm{d}}} $ 为拖曳力系数，${C_{\rm{d}}} = \left( {0.63 + \dfrac{{4.8}}{{\sqrt {{{{{{Re}}} }_{\text{p}}}} }}} \right)$ ；d为颗粒直径；v为颗粒速度矢量；$\;\chi $ 为修正系数，$\;\chi = 3.7 - 0.65\cdot \exp \left[ { - \dfrac{{{{(1.5 - {{\lg }}\;{{{{{Re}}} }_{\text{p}}})}^2}}}{2}} \right]$ 。颗粒雷诺数$ {{{{Re}}} _{\text{p}}} $ 表示为：$$ {{Re} _{\text{p}}} = \frac{{\varepsilon d{\rho _{\text{f}}}\left| {{\boldsymbol{u}} - {\boldsymbol{v}}} \right|}}{{{\mu _{\text{f}}}}} $$ (9) 式（4）中的颗粒对流体阻力项
$ {f_{{\text{pf}}}} $ 可以表示为：$$ {f_{{\text{pf}}}}{\text{ = }} - \frac{{\displaystyle\sum\limits_{{{i}} = 1}^{{n}} {{F_{{\text{d}i}}}} }}{{{V_{{\text{cell}}}}}} $$ (10) 1.4 虚拟球模型
局部孔隙率
$ \varepsilon $ 是颗粒相在流体控制方程中的重要表征，对于计算流体对颗粒的拖曳力至关重要。然而，在局部平均DEM–CFD耦合数值方法中，颗粒尺寸接近或大于最小流体网格尺寸时，会导致颗粒相占据的体积分数接近甚至大于网格体积，从而导致流体域孔隙度场不连续并使数值收敛困难[19]。因此，对于最小网格尺寸与最大颗粒粒径的比值，建议保持在4倍以上[20]。但对于真实案例模拟而言，复杂地形建模显然要求高分辨率的地形网格，因此限制了传统DEM–CFD耦合方法在滑坡–堵江–涌浪灾害链演化模拟中的应用。为此，本研究提出了一种基于虚拟球模型的局部孔隙率计算的改进方法，如图1所示。
图 1 虚拟球模型示意图Fig. 1 Schematic of the virtual sphere model
下载:
全尺寸图片
在扩展DEM–CFD耦合系统中，孔隙率计算的网格范围通过虚拟球判定，其质心与真实颗粒相同，直径为4倍真实颗粒直径。真实颗粒对孔隙率计算的贡献体积被均匀划分至虚拟球确定的相关网格。此外，同一网格被允许关联至不同虚拟球。因此，某一网格单元的局部孔隙率可通过式（1）计算得出[21]：
$$ \varepsilon = 1{{ - }}\sum\limits_{{i}} {\frac{{{V_{{\text{p,}i}}}}}{{{V_{{\text{vs,}i}}}}}} $$ (11) 式中，
${V_{{\text{p,}i}}}$ 为颗粒i的体积，${V_{{\text{vs,}i}}}$ 为颗粒i的所有相关网格总体积。1.5 扩展DEM–CFD耦合方法的计算流程
本研究中，扩展的DEM–CFD耦合方法基于离散元软件EDEM及CFD软件Fluent实现，其中耦合模块基于C语言开发。其计算流程如图2所示。
图 2 DEM–CFD耦合计算流程Fig. 2 Numerical solution strategy for the extended coupled DEM–CFD model下载:
全尺寸图片
通过DEM对滑坡体建模，并求解其运移过程。而在CFD模块中，首先通过引入的VOF模型追踪河流自由液面，其次通过压力隐式算子分裂（pressure-implicit with splitting of operators，PISO）算法迭代求解流体相控制方程，进而求解计算域内流体演化。DEM及CFD的时间步长分别满足Rayleigh时间及Courant–Friedrichs–Lewy（CFL）条件的约束。在耦合模块中，基于虚拟球模型计算局部孔隙率，并进一步计算包括浮力及拖曳力在内的耦合作用力。随后，将其传输至DEM模块进行接触力计算，而颗粒对流体的作用力则通过显式及隐式动量源项代入流体控制方程，从而实现双向耦合计算。
2. 模型验证
2.1 单颗粒沉降验证
为了验证本研究中扩展DEM–CFD耦合方法的有效性及准确性，采用单颗粒从空气入水沉降模型，并将颗粒沉降速度的数值解与理论解进行对比。单颗粒沉降验证案例如图3所示。从图3（a）可以看出：计算域为0.1 m×0.1 m×0.2 m的立方体，并采用4 mm尺寸的正六面体网格均匀离散。其中，上半部分为空气，下半部分为水。直径为 1 mm 的颗粒在距水面0.05 m处生成并从静止状态开始下落，模拟所需具体参数见表1。
图 3 单颗粒沉降验证Fig. 3 Verification of single particle sedmentation下载:
全尺寸图片
表 1 单颗粒沉降模拟参数Table 1 Parameters for single particle sedimentation simulation参数 取值 参数 取值 气相密度/(kg.m–3) 1 颗粒直径/mm 1 气相黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2500 水相密度/(kg·m–3) 1000 重力/(m·s–2) –9.81 水相黏度/(Pa∙s) 1×10–3 DEM时间步长/s 1×10–5 CFD时间步长/s 1×10–4 在此过程中，颗粒的运动可以表示为：
$$ {\;\;\;\;\frac{{4{\text{π}} {r^3}}}{3}{\rho _{\text{p}}}\frac{{{\rm{d}}{\boldsymbol{v}}}}{{{\rm{d}}t}} + \frac{{{\text{π}} {r^2}}}{2}{\rho _{\text{f}}}{C_{\text{d}}}{{\boldsymbol{v}}^2} = \frac{{{\text{π}} {r^3}}}{3}({\rho _{\text{p}}} - {\rho _{\text{f}}}){{g}} }$$ (12) 式中，
$ r $ 为颗粒半径，$ \;{\rho _{\text{p}}} $ 为颗粒密度。根据式（12），采用显式正向差分法分别计算不同时刻颗粒的沉降速度理论解。颗粒沉降速度–时间曲线如图3（b）所示，由图3（b）可见，颗粒沉降过程的理论解及数值解表现出良好的一致性，并都达到了0.134 m/s的稳定沉降速度。
2.2 不同网格尺寸下单颗粒沉降
保持颗粒直径不变，将网格尺寸设定为1～5 mm以改变网格尺寸与颗粒直径比，分别采用传统的中心点模型[22]和虚拟球模型计算孔隙率，并进行颗粒沉降过程模拟，结果如图4所示。
图 4 中心点模型和虚拟球模型中尺寸效应对比Fig. 4 Comparison of size ratio in central model and virtual sphere model下载:
全尺寸图片
由图4可知：在中心点模型中，随着尺寸比的不断减小，颗粒最终沉降速度的误差不断增大，这一结果反映了局部平均DEM–CFD方法的临界尺寸效应。与之相反，虚拟球模型在不同网格尺寸下的模拟结果相对稳定，其结果与理论解的相对误差在1%以内，且计算的颗粒最终沉降速度的精度最大能提高 40%。结果表明，虚拟球模型有效地消除了局部孔隙率计算中网格尺寸与颗粒直径比的限制。
3. 滑坡–堵江–涌浪灾害链模拟
3.1 模型设置及参数
滑坡体延坡面快速运移，随后侵入并阻塞河流，同时激发涌浪延河道及对岸边坡传播。在堵江过程中，碎屑物质进入河谷，随后到达对岸并沉积下来形成堰塞坝。滑坡–堵江–涌浪演化过程模拟模型如图5所示。由图5可知，滑坡–堵江–涌浪灾害链是一种典型的斜坡运动–河流系统的耦合作用现象，从模型尺度上建立灾害链演化的简化模拟模型，有助于理解滑坡–河流相互作用过程。
图 5 滑坡–堵江–涌浪演化过程模拟模型Fig. 5 Numerical modelling of the landslide river blocking and impulse wave evolution下载:
全尺寸图片
野外调查及统计数据表明，坡度在30°～45°的V型峡谷为滑坡堵江灾害的高发区域[23]。因此，模型设置如图5所示，研究区域为8 m×4 m的V型峡谷，其两侧岸坡均设置为35°，滑坡运动方向设定为与河道垂直。共计520 kg的滑坡体，初始位置设置在模型中段顶部，其初始速度为7 m/s。河流水深为0.3 m，流速为0.5 m/s。在 DEM 和 CFD 模块中，岸坡表面均施加无滑移边界条件，CFD模块中顶部边界设为开放边界。计算所需参数见表2，其中DEM参数已在之前的研究中经过标定和验证[24]。
表 2 滑坡坝形成过程模拟参数Table 2 Simulation parameters for landslide dam formation process参数 取值 参数 取值 空气密度/(kg∙m–3) 1 颗粒直径/mm [10, 40] 空气黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2700 水密度/(kg∙m–3) 1000 泊松比 0.25 水黏度/(Pa∙s) 0.001 剪切模量/Pa 1×109 重力/(m∙s–2) –9.81 恢复系数 0.5 CFD时间步长/s 1×10–4 摩擦系数 0.81 DEM时间步长/s 1×10–5 滚动摩擦系数 0.275 3.2 滑坡堵江及涌浪演化过程
图6为滑坡–堵江–涌浪演化的模拟结果。由图6可知：滑坡体侵入河流后，与河水、河谷之间的剧烈动能交换导致涌浪被激发。此外，碰撞下的能量耗散导致滑坡体前缘速度骤降，随后对中后部滑坡体产生阻碍作用致使其开始沿河道方向扩散。与此同时，滑坡体前缘受到中后部滑坡体的推挤作用，在河谷地形控制下转向，并混合涌浪一起沿对岸山谷爬升。当滑坡体–涌浪的混合体沿对侧山谷坡面爬升达到最大高度后便开始回落，滑坡体在河谷中稳定并沉积形成滑坡坝。在此过程中，河流被完全阻塞。随着上游来水持续被阻，坝前水位下降，而坝后水位不断上升并形成堰塞湖。
图 6 滑坡–堵江–涌浪演化模拟结果Fig. 6 Simulation results of the river blockage and impulse wave evolution下载:
全尺寸图片
图7显示了滑坡–堵江–涌浪灾害链过程模拟中的耦合作用力演化过程。由图7可以观察到：在滑坡体侵入河流的瞬间，颗粒–流体耦合作用力出现显著的短暂下降，这是由于河流对滑坡体的制动作用导致前缘速度降低。随着滑坡体持续侵入河流，耦合作用力继续上升，直到滑坡体到达谷底。滑坡体与谷底之间强烈碰撞产生的能量耗散致使耦合作用力不断减小。此后，随着堵江的完成及蓄水的增加，耦合作用力缓慢增长。
图 7 耦合作用力演化过程Fig. 7 Evolution process of the coupling forces下载:
全尺寸图片
图8为滑坡坝几何形态及涌浪在对岸最大爬升高度的演化过程。由图8可知：值得注意的是，涌浪在对岸的传播要略微晚于坝体，体现了滑坡体对涌浪传播的驱动作用。此外，从坝体沉积长度的演化过程可以发现，坝体上游的沉积长度显著小于坝体下游方向，上游一侧坝体的延伸速度也明显低于下游一侧。可见河流对滑坡体的运移沉积有明显的干扰作用，即河水能显著增强滑坡体顺流方向的运移，同时抑制逆流方向的运移。
图 8 坝体形态及涌浪高度演化过程Fig. 8 Evolution process of the dam geometry and the wave height下载:
全尺寸图片
4. 白格滑坡堵江案例模拟
2018年10月11日凌晨，西藏江达县与四川省白玉县交界处的白格村金沙江右岸边坡（东经98°42′17.98″，北纬31°4′56.41″）突发大规模滑动[25]。超过2.3×107 m3的失稳高位滑坡体高速滑动并侵入金沙江，激发了最大高度约160 m的涌浪沿对岸快速攀爬侵蚀，在对岸留下了最大标高达3 040 m的涌浪侵蚀区域，明显扩大了灾害影响范围[26-27]。随后，滑坡体阻塞金沙江并形成堰塞湖，形成了顺河长度达1500～2 000 m，横河宽度约510 m的堰塞坝，并在溃决前阻水达2.9×108 m3 [26]。
本研究以四川省测绘地理信息局提供的2017年实测数据的10 m分辨率灾前数字高程模型和2018年10月12日高性能无人机航空摄影测量的1 m分辨率灾后数字高程模型为基础建立模拟模型。“10·11”白格滑坡堵江模拟模型如图9所示，滑坡体被离散为2～5 m的颗粒，并通过DEM建模。CFD计算域内最小网格尺寸为2.5 m，考虑到灾害发生在金沙江汛期，参考调查信息设置金沙江初始水位高程为2890 m，上部由空气填充，顶部为开放边界。根据水文数据，确定入流量为1 680 m3/s[28]。在上游设置质量流入入口，并提前进行充分时间的模拟，以在耦合计算前达到稳定的初始河流条件。DEM参数则通过经验公式、文献[27,29-30]数据及参数标定获得，相关计算参数见表3。其中，摩擦系数可通过滑坡的高程落差及最大水平运动距离之比来确定，对于白格滑坡，其值可设定为0.45。此外，滑坡体的密度、泊松比及剪切模量等材料参数可参考前人研究的结果或结论。进一步地，基于坝体沉积形态开展参数反演，从而确定模拟参数取值。
图 9 “10·11”白格滑坡堵江模拟模型Fig. 9 Numerical modelling for the “10·11” Baige landslide simulation下载:
全尺寸图片
表 3 白格滑坡案例模拟参数Table 3 Simulation parameters for Baige landslide event参数 取值 参数 取值 空气密度/(kg·m–3) 1 颗粒半径/m [2,5] 空气黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2600 水密度/(kg·m–3) 1000 泊松比 0.2 水黏度/(Pa∙s) 0.001 剪切模量/Pa 1×109 重力/(m·s–2) –9.81 恢复系数 0.5 CFD时间步长/s 1×10–3 摩擦系数 0.45 DEM时间步长/s 1×10–4 滚动摩擦系数 0.05 图10为“10·11”白格滑坡–堵江–涌浪演化过程。由图10可见：在11 s左右，滑坡体前缘已抵达金沙江面，并有少量侵入。随后，滑坡体高速冲击河道，在滑坡体与河床的剧烈碰撞下，滑坡体前缘速度骤然下降，并沿对岸向上攀升。河流在滑坡体的剧烈扰动下激发涌浪，涌浪随后沿对岸岸坡向上传播。后部滑坡体持续侵入河道，一方面推动滑坡体前缘继续向上攀升，另一方面在狭窄山谷地形及前部滑坡体的制约下，沿着河道方向向上下游扩展。在33 s时，可以观察到河道中的滑坡体与坡面上滑坡体之间存在明显的速度差异。此时，涌浪波已经到达对岸，并保持在滑坡体的前方。值得注意的是，在上游方向沿山谷的初始延伸明显大于下游方向，上游侧山谷夹角小于下游侧，因此可知，地形约束加剧了滑坡体延河谷方向的运动。在44 s左右，涌浪波到达对岸最大高程，并开始回落到河道中。同时，可以观察到涌浪波也沿河道向上游和下游方向传播（图10（e））。随后，大部分滑坡体沉积在山谷中，形成了滑坡坝，而两侧的坝体前缘继续沿山谷方向缓慢延伸（图10（f））。
图 10 “10·11”白格滑坡–堵江–涌浪演化过程Fig. 10 River blockage and impulse wave evolution of the “10·11” Baige landslide下载:
全尺寸图片
“10·11”白格滑坡灾害的滑坡速度演化过程及最大涌浪高度演化过程如图11所示。由图11可以发现：其速度演化包括启动、加速、减速及沉积4个阶段，最大平均速度为47 m/s。此外，其涌浪高度演化曲线峰值平缓，可见受地形影响，涌浪爬高在对侧岸坡上不同位置先后达到峰值。最大涌浪爬高为海拔3 033 m，超过水面高程153 m，十分接近调查所得的160 m。
图 11 滑坡速度及涌浪高度演化过程Fig. 11 Evolution of the landslide velocity and the wave height下载:
全尺寸图片
“10·11”白格滑坡事件模拟结果如图12所示。通过观察图12中的累积滑坡运移路径及涌浪侵蚀面积可以发现，其模拟得到的滑坡最终沉积形态整体呈弓形，顺河方向沉积长度为1 759.01 m，横河方向沉积长度为544.14 m，均与现场调查值一致。模拟得到的最大坝顶高程为海拔2 990 m，与海拔2 985 m的实际最大坝顶高程十分接近。此外，通过与实际灾害影响边界对比可进一步发现，整体模拟灾害影响范围与实地调查获得的灾害边界高度吻合，进一步说明了模拟结果的可靠性。
图 12 “10·11”白格滑坡事件模型结果Fig. 12 Simulation result of the “10·11” Baige landslide event下载:
全尺寸图片
5. 结 论
本文提出了一种扩展DEM–CFD耦合数值方法，并将其应用至滑坡堵江及涌浪过程模拟和演化机制研究，进一步模拟了“10·11”金沙江白格滑坡堵江事件。具体结论如下：
1）通过引入VOF模型及虚拟球模型，所提出的扩展DEM–CFD耦合方法实现了对河流自由液面演化的追踪，并克服了网格颗粒临界尺寸的限制，实现了耦合作用力的精确求解，能够满足滑坡体粗颗粒材料的建模及高分辨率地形网格的要求。
2）模拟结果表明，滑坡坝形成过程涉及强烈的滑坡–河流相互作用现象。滑坡体驱动涌浪在对岸爬升及河道内传播；而河流显著增加了滑坡体的动能耗散，且在上下游方向上对滑坡体的沉积产生相反的影响。
3）所提出的扩展DEM–CFD耦合数值方法准确再现了“10·11”金沙江白格滑坡堵江及涌浪演化过程，模拟得到的滑坡运移路径、沉积形态及涌浪侵蚀面积与实地调查结果保持了较高的一致性，从而进一步证明了方法的有效性。
-
图 1 虚拟球模型示意图
Fig. 1 Schematic of the virtual sphere model
下载:
全尺寸图片
图 2 DEM–CFD耦合计算流程
Fig. 2 Numerical solution strategy for the extended coupled DEM–CFD model
下载:
全尺寸图片
图 3 单颗粒沉降验证
Fig. 3 Verification of single particle sedmentation
下载:
全尺寸图片
图 4 中心点模型和虚拟球模型中尺寸效应对比
Fig. 4 Comparison of size ratio in central model and virtual sphere model
下载:
全尺寸图片
图 5 滑坡–堵江–涌浪演化过程模拟模型
Fig. 5 Numerical modelling of the landslide river blocking and impulse wave evolution
下载:
全尺寸图片
图 6 滑坡–堵江–涌浪演化模拟结果
Fig. 6 Simulation results of the river blockage and impulse wave evolution
下载:
全尺寸图片
图 7 耦合作用力演化过程
Fig. 7 Evolution process of the coupling forces
下载:
全尺寸图片
图 8 坝体形态及涌浪高度演化过程
Fig. 8 Evolution process of the dam geometry and the wave height
下载:
全尺寸图片
图 9 “10·11”白格滑坡堵江模拟模型
Fig. 9 Numerical modelling for the “10·11” Baige landslide simulation
下载:
全尺寸图片
图 10 “10·11”白格滑坡–堵江–涌浪演化过程
Fig. 10 River blockage and impulse wave evolution of the “10·11” Baige landslide
下载:
全尺寸图片
图 11 滑坡速度及涌浪高度演化过程
Fig. 11 Evolution of the landslide velocity and the wave height
下载:
全尺寸图片
图 12 “10·11”白格滑坡事件模型结果
Fig. 12 Simulation result of the “10·11” Baige landslide event
下载:
全尺寸图片
表 1 单颗粒沉降模拟参数
Table 1 Parameters for single particle sedimentation simulation
参数 取值 参数 取值 气相密度/(kg.m–3) 1 颗粒直径/mm 1 气相黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2500 水相密度/(kg·m–3) 1000 重力/(m·s–2) –9.81 水相黏度/(Pa∙s) 1×10–3 DEM时间步长/s 1×10–5 CFD时间步长/s 1×10–4 表 2 滑坡坝形成过程模拟参数
Table 2 Simulation parameters for landslide dam formation process
参数 取值 参数 取值 空气密度/(kg∙m–3) 1 颗粒直径/mm [10, 40] 空气黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2700 水密度/(kg∙m–3) 1000 泊松比 0.25 水黏度/(Pa∙s) 0.001 剪切模量/Pa 1×109 重力/(m∙s–2) –9.81 恢复系数 0.5 CFD时间步长/s 1×10–4 摩擦系数 0.81 DEM时间步长/s 1×10–5 滚动摩擦系数 0.275 表 3 白格滑坡案例模拟参数
Table 3 Simulation parameters for Baige landslide event
参数 取值 参数 取值 空气密度/(kg·m–3) 1 颗粒半径/m [2,5] 空气黏度/(Pa∙s) 1×10–5 颗粒密度/(kg·m–3) 2600 水密度/(kg·m–3) 1000 泊松比 0.2 水黏度/(Pa∙s) 0.001 剪切模量/Pa 1×109 重力/(m·s–2) –9.81 恢复系数 0.5 CFD时间步长/s 1×10–3 摩擦系数 0.45 DEM时间步长/s 1×10–4 滚动摩擦系数 0.05 -
[1] 黄润秋,许强.中国典型灾难性滑坡[M].北京:科学出版社,2008:1-4. [2] 年廷凯,吴昊,陈光齐,等.堰塞坝稳定性评价方法及灾害链效应研究进展[J].岩石力学与工程学报,2018,37(8):1796–1812. doi: 10.13722/j.cnki.jrme.2017.1655 Nian Tingkai,Wu Hao,Chen Guangqi,et al.Research progress on stability evaluation method and disaster chain effect of landslide dam[J].Chinese Journal of Rock Mechanics and Engineering,2018,37(8):1796–1812 doi: 10.13722/j.cnki.jrme.2017.1655 [3] 郑鸿超,石振明,彭铭,等.崩滑碎屑体堵江成坝研究综述与展望[J].工程科学与技术,2020,52(2):19–28. doi: 10.15961/j.jsuese.201900946 Zheng Hongchao,Shi Zhenming,Peng Ming,et al.Review and prospect of the formation mechanism of landslide dams caused by landslide and avalanche debris[J].Advanced Engineering Sciences,2020,52(2):19–28 doi: 10.15961/j.jsuese.201900946 [4] 戴兴建,殷跃平,邢爱国.易贡滑坡–碎屑流–堰塞坝溃坝链生灾害全过程模拟与动态特征分析[J].中国地质灾害与防治学报,2019,30(5):1–8. doi: 10.16031/j.cnki.issn.1003-8035.2019.05.01 Dai Xingjian,Yin Yueping,Xing Aiguo.Simulation and dynamic analysis of Yigong rockslide-debris avalanche-dam breaking disaster chain[J].The Chinese Journal of Geological Hazard and Control,2019,30(5):1–8 doi: 10.16031/j.cnki.issn.1003-8035.2019.05.01 [5] Crosta G B,Imposimato S,Roddeman D.Landslide spreading,impulse water waves and modelling of the vajont rockslide[J].Rock Mechanics and Rock Engineering,2016,49(6):2413–2436. doi: 10.1007/s00603-015-0769-z [6] 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 [7] Li Ji.Barrier lake formation due to landslide impacting a river:A numerical study using a double layer-averaged two-phase flow model[J].Applied Mathematical Modelling,2020,80:574–601. doi: 10.1016/j.apm.2019.11.031 [8] Liu Wei,He Siming.A two-layer model for simulating landslide dam over mobile river beds[J].Landslides,2016,13(3):565–576. doi: 10.1007/s10346-015-0585-2 [9] Liu Wei,He Siming.Dynamic simulation of a mountain disaster chain:Landslides,barrier lakes,and outburst floods[J].Natural Hazards,2018,90(2):757–775. doi: 10.1007/s11069-017-3073-2 [10] 徐文杰.滑坡涌浪流–固耦合分析方法与应用[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 [11] 年廷凯,沈月强,郑德凤,等.海底滑坡链式灾害研究进展[J].工程地质学报,2021,29(6):1657–1675. doi: 10.13544/j.cnki.jeg.2021-0815 Nian Tingkai,Shen Yueqiang,Zheng Defeng,et al.Research advances on the chain disasters of submarine landslides[J].Journal of Engineering Geology,2021,29(6):1657–1675 doi: 10.13544/j.cnki.jeg.2021-0815 [12] Zhao Tao,Dai Feng,Xu Nuwen.Coupled DEM–CFD investigation on the formation of landslide dams in narrow rivers[J].Landslides,2017,14(1):189–201. doi: 10.1007/s10346-015-0675-1 [13] 年廷凯,张放,郑德凤,等.基于计算流体力学–离散元法耦合方法的海底黏性滑坡体运动行为模拟[J].岩土力学,2022,43(11):3174–3184. doi: 10.16285/j.rsm.2021.1856 Nian Tingkai,Zhang Fang,Zheng Defeng,et al.Numerical simulation on the movement behavior of viscous submarine landslide based on coupled computational fluid dynamics-discrete element method[J].Rock and Soil Mechanics,2022,43(11):3174–3184 doi: 10.16285/j.rsm.2021.1856 [14] Cundall P A,Strack O D L.A discrete numerical model for granular assemblies[J].Géotechnique,1979,29(1):47–65. doi: 10.1680/geot.1979.29.1.47 [15] Anderson T B,Jackson R.Fluid mechanical description of fluidized beds.equations of motion[J].Industrial & Engineering Chemistry Fundamentals,1967,6(4):527–539. doi: 10.1021/i160024a007 [16] Hirt C W,Nichols B D.Volume of fluid (VOF) method for the dynamics of free boundaries[J].Journal of Computational Physics,1981,39(1):201–225. doi: 10.1016/0021-9991(81)90145-5 [17] Nian Tingkai,Li Dongyang,Liang Qiuhua,et al.Multi-phase flow simulation of landslide dam formation process based on extended coupled DEM–CFD method[J].Computers and Geotechnics,2021,140:104438. doi: 10.1016/j.compgeo.2021.104438 [18] Di Felice R.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 [19] 彭恺然,刘红帅,平新雨,等.CFD–DEM耦合模拟中拖曳力模型精度[J].吉林大学学报(地球科学版),2021,51(5):1400–1407. doi: 10.13278/j.cnki.jjuese.20210017 Peng Kairan,Liu Hongshuai,Ping Xinyu,et al.Accuracy of drag force models in the CFD–DEM method[J].Journal of Jilin University(Earth Science Edition),2021,51(5):1400–1407 doi: 10.13278/j.cnki.jjuese.20210017 [20] Zhao T,Houlsby G T,Utili S.Investigation of granular batch sedimentation via DEM–CFD coupling[J].Granular Matter,2014,16(6):921–932. doi: 10.1007/s10035-014-0534-0 [21] Li Dongyang,Zheng Defeng,Wu Hao,et al.Numerical simulation on the longitudinal breach process of landslide dams using an improved coupled DEM–CFD method[J].Frontiers in Earth Science,2021,9:673249. doi: 10.3389/feart.2021.673249 [22] Kloss C,Goniva C,Hager A,et al.Models,algorithms and validation for opensource DEM and CFD–DEM[J].Progress in Computational Fluid Dynamics,an International Journal,2012,12(2/3):140. doi: 10.1504/pcfd.2012.047457 [23] Liao Haimei,Yang Xingguo,Lu Gongda,et al.Experimental study on the formation of landslide dams by fragmentary materials from successive rock slides[J].Bulletin of Engineering Geology and the Environment,2020,79(3):1591–1604. doi: 10.1007/s10064-019-01651-7 [24] Li Dongyang,Nian Tingkai,Wu Hao,et al.A predictive model for the geometry of landslide dams in V-shaped valleys[J].Bulletin of Engineering Geology and the Environment,2020,79(9):4595–4608. doi: 10.1007/s10064-020-01828-5 [25] 周礼,范宣梅,许强,等.金沙江白格滑坡运动过程特征数值模拟与危险性预测研究[J].工程地质学报,2019,27(6):1395–1404. doi: 10.13544/j.cnki.jeg.2019-037 Zhou Li,Fan Xuanmei,Xu Qiang,et al.Numerical simulation and hazard prediction on movement process characteristics of Baige landslide in Jinsha River[J].Journal of Engineering Geology,2019,27(6):1395–1404 doi: 10.13544/j.cnki.jeg.2019-037 [26] Ouyang Chaojun,An Huicong,Zhou Shu,et al.Insights from the failure and dynamic characteristics of two sequential landslides at Baige Village along the Jinsha River,China[J].Landslides,2019,16(7):1397–1414. doi: 10.1007/s10346-019-01177-9 [27] Zhang Shilin.Dynamics and emplacement mechanisms of the successive Baige landslides on the upper reaches of the Jinsha River,China[J].Engineering Geology,2020,278:105819. doi: 10.1016/j.enggeo.2020.105819 [28] Zhang Limin,Xiao Te,He Jian,et al.Erosion-based analysis of breaching of baige landslide dams on the Jinsha River,China,in 2018[J].Landslides,2019,16(10):1965–1979. doi: 10.1007/s10346-019-01247-y [29] Wang Wenpei,Yin Yueping,Zhu Sainan,et al.Investigation and numerical modeling of the overloading-induced catastrophic rockslide avalanche in Baige,Tibet,China[J].Bulletin of Engineering Geology and the Environment,2020,79(4):1765–1779. doi: 10.1007/s10064-019-01664-2 [30] An Huicong,Ouyang Chaojun,Zhou Shu.Dynamic process analysis of the Baige landslide by the combination of DEM and long-period seismic waves[J].Landslides,2021,18(5):1625–1639. doi: 10.1007/s10346-020-01595-0
