猴子岩水库色玉滑坡涌浪灾害链CFD–DEM耦合数值模拟

2023-02-19肖华波王泽皓石伟明王东坡欧阳朝军

工程科学与技术 2023年1期

肖华波，王泽皓，石伟明，王东坡，欧阳朝军

(1.中国电建集团成都勘测设计研究院有限公司，四川成都 610072；2.成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川成都 610059；3.中国科学院水利部成都山地灾害与环境研究所，四川成都 610041)

位于河道、水库周边的大型滑坡、崩塌等地质灾害一旦发生，可能造成堵江并产生巨大的涌浪，对人民的生命财产安全、船舶安全及水电站的运行等构成威胁[1–3]。例如：2003年7月，位于三峡库区的千将坪地区突然发生了一起方量为1.5×107m3的滑坡，该滑坡形成的堰塞坝将青干河堵塞，造成24名人员死亡，直接经济损失约700万美元[4]。2013年7月，云南省永善县黄华镇黄坪村发生山体滑坡，滑坡体冲入金沙江激起15～20 m高的涌浪，造成12人失踪，多部车辆损毁[5]。2015年6月，重庆市巫山县红岩子发生滑坡，其方量约为2.3×105m3，坡体入水后产生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模型中，流体界面的演化通过求解以下方程来描述：

式中，i为网格数，ρi、µi分别为第i个网格的流体密度和黏度。

在CFD–DEM耦合模型中，流体的控制方程为Navier–Stokes方程，如式（5）～（6）所示：

1.2 离散单元法控制方程

在CFD–DEM耦合模型中，颗粒的运动状态由牛顿第二定律进行描述，包括平动及转动两部分。除颗粒之间的相互作用力和力矩外，平动及转动两部分控制方程中还分别添加了流体作用于颗粒上的力和力矩，方程如下：

式（9）～（10）中，mi为颗粒i的质量，ui为颗粒i的平动速度，nj为与颗粒j接触的颗粒总数，Fij为颗粒i与颗粒j的相互作用力，Ffi为流体作用在颗粒i上的合力，Ii为转动惯量， θi为颗粒i的旋转角速度，Mij为颗粒j作用于颗粒i的力矩，Mif为流体作用于颗粒i的力矩，g为重力加速度。

颗粒间相互作用力的计算采用基于Hertz接触理论的Hertz–Mindlin模型[21]。当两个颗粒的接触距离小于两者半径之和时，将产生相互作用力，法向包括弹簧力和阻尼力两项，切向包括剪切力和阻尼力两项。该模型计算公式如式（11）所示：

式中，kn、kt分别为法向、切向弹性系数， δn为法向重叠量， δt为切向位移， γn、 γt分别为法向、切向阻尼系数，vn为法向相对速度，vt为切向相对速度。

1.3 流体–颗粒相互作用

在CFD–DEM耦合模型中，流体–颗粒间的相互作用通过在控制方程中交换相互作用力来考虑。本文主要考虑拖拽力和浮力的影响[21]，即：

式中：Cd为拖拽力系数；up为颗粒速度； χ为孔隙率的一个经验修正系数，此系数使得该关系式能够适用于更大范围的雷诺数的情况。

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 颗粒堆积体坍塌–涌浪试验计算参数Tab. 1 Calculation parameters of the particle accumulation collapse–surge experiment

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×104m3，失稳高度最大达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。

表2 色玉滑坡计算参数Tab. 2 Calculation parameters of Seyu landslide

图5 色玉滑坡3维模型Fig. 5 3D model of Seyu landslide

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可以看出，涌浪运动过程主要分为3个阶段，具体如下：

图11 涌浪演进过程剖面Fig. 11 Profile of surge evolution process

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个波峰，后续波峰峰值衰减较大。沿程涌浪是一个逐渐衰减过程，最大涌浪高度随着传播距离的增加而减小。