山区河流型水库碎土石滑坡涌浪三维数值模拟研究

赵 树 正， 田 宇 霖， 曹 山， 薛 宏 程

(重庆交通大学 河海学院，重庆 400016)

0 引 言

随着国家长江经济带黄金水道的建设，沿江两岸已探明的滑坡就有千余处[1]，若有巨型滑坡体高速滑入水中，其形成的涌浪不仅会影响过往船舶的航行安全，还会对沿岸居民的生命财产安全造成严重威胁[2]，因此有效预测滑坡涌浪这类次生灾害的影响范围具有重要的科学意义。目前研究滑坡涌浪的方法主要有物理模型试验和数值模拟[3]，模型试验能够有效地观测到涌浪生成、传播和爬坡等过程，如岳书波等[4]和Heller等[5]通过物理模型试验，总结了二维水槽中涌浪的生成形态和初始波幅经验公式；殷坤龙等[6]根据实际三维物理模型试验，发现涌浪在三维空间的波幅远小于二维情况。但是物理模型存在缩尺效应显著、数据测量较难和试验成本较高等问题。因此，部分研究人员利用数值模拟方法对滑坡涌浪形态变化、传播过程开展模拟研究，如孔增增[7]采用数值模拟等方法对V型河谷滑坡涌浪的最大波面高度以及涌浪的爬高进行模拟，优化了涌浪最大爬坡高度经验公式；Huang等[8]用非线性Boussinesq水波模型研究了三峡库区红岩子滑坡涌浪。研究结果表明涌浪传播过程受河谷地貌条件影响较大，滑坡对岸的城镇是受涌浪影响的主要区域。数值模拟方法作为一种新兴的研究手段，在解决地形复杂的水库和河道滑坡涌浪问题时具有明显优势。但是滑坡涌浪属于典型的流固耦合问题，滑坡体高速滑入水中还会卷入大量空气[9,10]，对计算区域网格划分要求较高，尤其是进行三维数值模拟难度较大。

本文利用CFD数值仿真技术，基于RNGκ-ε紊流模型与VOF方法，将碎土石滑坡体视为固体颗粒和液体的混合物，对已有模型试验研究成果进行数值模拟，发现数值模拟结果与模型试验吻合较好。在此基础上，模拟了某实际山区河流型水库中碎土石滑坡涌浪的生成和传播全过程，对涌浪的形态变化、波幅变化和传播速度进行了分析。研究结果可为实际工程中潜在的滑坡涌浪灾害预测提供科学依据。

1 数学模型

(1)

(2)

(3)

(4)

本文封闭连续方程(公式1)和运动方程(公式2)采用经典的RNGκ-ε紊流模型(公式3和公式4)，并利用有限差分法将控制方程离散为代数方程，从而进行数值求解。涌浪自由面的追踪采用VOF方法，该方法不考虑气体对于模型的影响，而只考虑纯液体单元的影响，即在一个计算网格内，所有流体相的体积分数之和为1。另外，滑坡体用颗粒体来模拟，是固体颗粒和液体的混合物，这种混合物被视为不可压缩流体，其边界可以是自由表面。

2 数值模拟方法可靠性分析

根据已有模型试验研究成果建立滑坡涌浪计算模型，对散粒体滑坡入水产生涌浪的过程进行数值模拟，以验证数值模拟方法的可靠性。

2.1 计算模型及网格划分

在长10 m，宽0.60 m，高1 m的矩形水槽内开展了散粒体滑坡涌浪试验研究[11]。该试验的滑坡体材料为沙石颗粒(72%的颗粒粒径为0.01～2 mm，18%的颗粒粒径为2～5 mm，10%的颗粒粒径为5～15 mm)，配比后的滑坡体密度为2 100 kg/m3，滑坡体的入水速度通过特制的滑动装置来控制，各项参数如图1所示，其中h为水深，vs为滑坡体入水速度，ls为滑坡体长度，s为滑坡体厚度，α为入水角度[11]。

图1 模型试验中各参数示意图

如图2所示，本算例的计算区域与上述模型试验一致，计算区域全部为四边形网格，网格间距为20 mm，网格总数约100万；边界条件包括速度进口、压力出口和无滑移固壁边界；初始相为水、空气和滑坡体。现以试验[11]中的一组试验参数值和数据作为数值模拟的依据，该组试验滑坡体下滑弗氏数F=1.65，滑坡体相对厚度S=0.23，水深h=0.3 m，滑坡体下滑速度vs=2.83 m/s，滑坡体长度ls=0.75 m，入水角度α=33°。滑坡体的平均颗粒粒径d取1.8 mm，滑坡体紧密堆积体积分数取0.63，自然休止角取34°。

图a 计算区域示意图 图b 局部网格示意图图2 数值模拟计算区域及网格示意图

2.2 计算结果对比与分析

图3(a)～3(c)左侧分别为模型试验在相对无量纲时间t*=0.46、t*=2.52和t*=3.43时的初始涌浪形态(以滑坡体与水面恰好接触的时刻为t*=0)。而图3(a)～3(c)右侧分别选取了与 试验[11]时刻对应的三个数值模拟涌浪形态，从图3(a)(t*=0.46)可以看出数值模拟中的滑坡体冲击水面后，水体的运动趋势与试验得到的涌浪形态基本一致，但由于滑坡体冲击水面形成的薄水舌和溅起的水花属于强非线性水气二相流，在模拟这种水气间强烈混掺的情况时VOF法存在局限性，故未能模拟出溅射的水花部分；从图3(b)(t*=2.52)中可以看出当数值模拟的滑坡体完全滑入水中后，滑坡体的扩散形态与文献中的滑坡运动形态较为接近；由于y轴方向上数值模拟采用了对称边界条件，而文献中水槽两侧的钢化玻璃对涌浪传播起到了限制作用，所以从图3(c)(t*=3.43)中可以看到滑坡体堆积形态和最大波幅高度与模型试验相比基本吻合，但最大波幅出现的位置比试验传播的稍远。总体上看，数值模拟得到的滑坡体初始堆积形态与试验结果吻合较好，波幅计算精度能够满足要求。

(a) t*=0.46

(b)t*=2.52

(c) t*=3.43图3 不同时刻下的涌浪形态对比[11]

3 实际工程应用

某水电工程位于澜沧江一级支流的下游河段，是以发电为主的混合式开发水电站，属于典型的山区河流型水库。工程主要建筑物包括拦河大坝、泄水建筑物及引水发电建筑物等，水库总库容3 403万m3，电站总装机容量100 MW，为中等规模Ⅲ等水电工程。工程勘测阶段发现坝址上游左岸约200 m处有一古滑坡体(滑坡后缘高出正常蓄水位约60 m)，由于水库蓄水后可能会引起两岸山体地下水位抬升，导致岸坡岩土体强度降低，从而诱发滑坡。故本文针对该水电工程库区蓄水后潜在的滑坡涌浪灾害开展了数值模拟研究。

图4分别给出了时刻为1.5 s、3 s、4.5 s、7.5 s、10.5 s和13.5 s时库区涌浪的波高变化过程(静止水面高程设置为0)。当滑坡体滑入水中，在近场处形成近30 m高的初始涌浪(t=4.5 s时)，随着涌浪向四周扩散传播，其沿滑坡体两侧河岸传播时呈现出明显的衰减趋势，但主波区波高仍可达十余米。当涌浪传播到对岸山体时，波幅衰减到十米内，但涌浪最大波峰离坝体较远。

图4 库区涌浪过程波高变化过程 (a)t=1.5 s; (b)t=3 s; (c)t=4.5 s; (d)t=7.5 s; (e)t=10.5 s; (f)t=13.5 s

图5给出了不同时刻下涌浪传播速度的变化过程。滑坡体在入水时具有较高的下滑速度，前缘入水时的速度可达18 m/s，后缘在7.5 s后也基本滑入水中，入水时最大速度可达24 m/s，具有较大的动能。初始涌浪形成后，其传播波速在库区可达10～16 m/s，若有船舶经过，海事风险极大。

4 结 语

采用RNGκ-ε紊流模型与VOF方法对散粒体滑坡涌浪进行数值模拟。通过与已有文献的试验结果对比分析，发现数值模拟得到的滑坡体水下初始堆积形态与试验结果吻合较好，且波幅计算精度较高。基于此，建立了某山区河流型水电工程库区的三维模型，模拟了碎土石滑坡涌浪在库区内的生成与传播过程，并分析了涌浪的波幅和波速特征。研究结果表明，本文的数值模拟方法能够较好地模拟山区河流型水库内的碎土石滑坡涌浪生成与传播过程。通过数值模拟获得的涌浪波幅和传播速度等参数，可以为大坝安全评估和涌浪灾害影响范围预测提供依据，还可为通航河流航行船舶的避险范围提供参考。值得说明的是，将碎土石滑坡体视为固体颗粒和液体的混合物，在滑坡下滑的初始阶段能够较为真实地反映滑坡体在岸坡底部的堆积形态，但依然无法改变其流体的本质，在库区水面波动趋于稳定后，滑坡体会水平堆积于河道底部，但这并不影响前期涌浪生成和传播数据的采集。另外，VOF方法在模拟强非线性水气二相流时具有局限性，滑坡体入水时产生的水舌飞溅较难模拟，在实际工程应用中可建立更加细致的三维模型和网格，尽可能确保滑坡体初始形态与实际情况吻合，从而进一步提高模拟计算精度。

图5 库区涌浪过程水面速度变化过程 (a)t=1.5 s; (b)t=3 s; (c)t=4.5 s; (d)t=7.5 s; (e)t=10.5 s; (f)t=13.5 s