滑坡–堵江–涌浪灾害链模拟的DEM–CFD耦合分析方法及其应用

李东阳，年廷凯*，吴 昊，张彦君

(1.大连理工大学 土木工程学院，辽宁 大连 116024；
2.中国科学院 水利部成都山地灾害与环境研究所，四川 成都 610299；
3.中国地质调查局武汉地质调查中心，湖北 武汉 430205)

中国西南地区，尤其是青藏高原地区，沟壑林立且河流发育，一直是滑坡类灾害的高发区[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.1 控制方程

本研究所提出的DEM–CFD耦合数值方法中，滑坡体通过DEM建模，并通过求解牛顿第二定律来模拟滑坡体的运移过程[14]：

式（1）～（2）中，mi、vi和ωi分别为颗粒质量、速度及角速度，g为重力加速度，Fij和Mij分别为颗粒的接触作用力及扭矩，Ff为流体对颗粒的耦合作用力。

对于流体相，通过CFD求解颗粒相体积分数修正后的局部平均Navier–Stokes方程[15]：

式（3）～（4）中， ε为局部孔隙率，u、 ρf、 µf分别为流体速度、密度及黏度，p为流体压力，fpf为表征颗粒对流体阻力的动量源项，式（4）右侧其余3项则依次为压力梯度项、应力项及重力项。

1.2 VOF模型

在本研究中，将VOF法引入传统的局部平均DEM–CFD耦合方法中，从而允许对多相流体介质建模，涌浪的形成及演化过程通过对多相流体交界面的追踪来实现。本研究的扩展DEM–CFD耦合方法中，每个流体单元被分为水相和空气相，并用 αw和 αa分别表示相应相的体积分数。其中，水相的体积分数满足0 ≤αw≤1，则空气相的体积分数可表示为[16]：

因此， 0 <αw<1被定义为自由液面的存在。式（4）中的流体密度 ρf及黏度 µf，参考两相流体性质，通过线性插值计算。通过求解基于流体速度的输运方程，自由液面运动可计算如下[17]：

式中，ur为两相介质相对速度，ur=u1−u2。值得注意的是，方程中已经考虑颗粒相的存在。左侧第1项中ε的时间导数考虑到了体积分数变化对液相体积分数的影响。式（6）的前两项与式（4）具有相同的平流形式，第3项为人工项，用于减少数值耗散，以达到使自由液面更为清晰的目的。

1.3 耦合作用力计算

本研究中，流体与颗粒之间的耦合作用力Ff包括浮力Fb和拖曳力Fd，其中，浮力Fb基于压力梯度求解，可表示为：

式中，Vp为颗粒体积。拖曳力是由流体域颗粒间的相对运动引起的，作用于颗粒质心，其方向与相对运动方向相反。通过Di Felice模型来计算拖曳力[18]：

1.4 虚拟球模型

局部孔隙率 ε是颗粒相在流体控制方程中的重要表征，对于计算流体对颗粒的拖曳力至关重要。然而，在局部平均DEM–CFD耦合数值方法中，颗粒尺寸接近或大于最小流体网格尺寸时，会导致颗粒相占据的体积分数接近甚至大于网格体积，从而导致流体域孔隙度场不连续并使数值收敛困难[19]。因此，对于最小网格尺寸与最大颗粒粒径的比值，建议保持在4倍以上[20]。但对于真实案例模拟而言，复杂地形建模显然要求高分辨率的地形网格，因此限制了传统DEM–CFD耦合方法在滑坡–堵江–涌浪灾害链演化模拟中的应用。为此，本研究提出了一种基于虚拟球模型的局部孔隙率计算的改进方法，如图1所示。

图1 虚拟球模型示意图Fig. 1 Schematic of the virtual sphere model

在扩展DEM–CFD耦合系统中，孔隙率计算的网格范围通过虚拟球判定，其质心与真实颗粒相同，直径为4倍真实颗粒直径。真实颗粒对孔隙率计算的贡献体积被均匀划分至虚拟球确定的相关网格。此外，同一网格被允许关联至不同虚拟球。因此，某一网格单元的局部孔隙率可通过式（1）计算得出[21]：

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.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 单颗粒沉降模拟参数Tab. 1 Parameters for single particle sedimentation simulation

根据式（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.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 滑坡坝形成过程模拟参数Tab. 2 Simulation parameters for landslide dam formation process

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

2018年10月11日凌晨，西藏江达县与四川省白玉县交界处的白格村金沙江右岸边坡（东经98°42′17.98″，北纬31°4′56.41″）突发大规模滑动[25]。超过2.3×107m3的失稳高位滑坡体高速滑动并侵入金沙江，激发了最大高度约160 m的涌浪沿对岸快速攀爬侵蚀，在对岸留下了最大标高达3 040 m的涌浪侵蚀区域，明显扩大了灾害影响范围[26–27]。随后，滑坡体阻塞金沙江并形成堰塞湖，形成了顺河长度达1 500～2 000 m，横河宽度约510 m的堰塞坝，并在溃决前阻水达2.9×108m3[26]。

本研究以四川省测绘地理信息局提供的2017年实测数据的10 m分辨率灾前数字高程模型和2018年10月12日高性能无人机航空摄影测量的1 m分辨率灾后数字高程模型为基础建立模拟模型。“10·11”白格滑坡堵江模拟模型如图9所示，滑坡体被离散为2～5 m的颗粒，并通过DEM建模。CFD计算域内最小网格尺寸为2.5 m，考虑到灾害发生在金沙江汛期，参考调查信息设置金沙江初始水位高程为2 890 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 白格滑坡案例模拟参数Tab. 3 Simulation parameters for Baige landslide event

图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

本文提出了一种扩展DEM–CFD耦合数值方法，并将其应用至滑坡堵江及涌浪过程模拟和演化机制研究，进一步模拟了“10·11”金沙江白格滑坡堵江事件。具体结论如下：

1）通过引入VOF模型及虚拟球模型，所提出的扩展DEM–CFD耦合方法实现了对河流自由液面演化的追踪，并克服了网格颗粒临界尺寸的限制，实现了耦合作用力的精确求解，能够满足滑坡体粗颗粒材料的建模及高分辨率地形网格的要求。

2）模拟结果表明，滑坡坝形成过程涉及强烈的滑坡–河流相互作用现象。滑坡体驱动涌浪在对岸爬升及河道内传播；
而河流显著增加了滑坡体的动能耗散，且在上下游方向上对滑坡体的沉积产生相反的影响。

3）所提出的扩展DEM–CFD耦合数值方法准确再现了“10·11”金沙江白格滑坡堵江及涌浪演化过程，模拟得到的滑坡运移路径、沉积形态及涌浪侵蚀面积与实地调查结果保持了较高的一致性，从而进一步证明了方法的有效性。

猜你喜欢演化过程滑坡体滑坡模因论视角下韩语“먹다”表“喝”动作演化过程研究韩国语教学与研究(2021年2期)2021-11-24时间非齐次二态量子游荡的演化过程分析数学物理学报(2021年4期)2021-08-30滑坡推力隐式解与显式解对比分析——以河北某膨胀土滑坡为例河北地质(2021年1期)2021-07-21重庆万盛石林的形成时代及发育演化过程四川地质学报(2020年2期)2020-05-31秦巴山区牟牛沟滑坡体治理施工技术城市建筑空间(2018年12期)2018-08-26浅谈公路滑坡治理北方交通(2016年12期)2017-01-15强震下紫坪铺坝前大型古滑坡体变形破坏效应西南交通大学学报(2016年6期)2016-05-04堆载诱发型滑坡变形演化机理的模型试验研究石家庄铁道大学学报(自然科学版)(2016年1期)2016-04-20“监管滑坡”比“渣土山”滑坡更可怕山东青年(2016年3期)2016-02-28基于灰色系统理论的滑坡体变形规律研究地矿测绘(2015年3期)2015-12-22 相关热词搜索：涌浪，耦合，灾害，