基于流固耦合方法的滑坡涌浪数值模拟

司鹏飞，田利勇

(1.上海市水利工程设计研究院有限公司，上海 200061; 2.上海市水务局防汛减灾工程技术研究中心，上海 200061)

库区滑坡是水库蓄水运行期间普遍存在的问题。滑坡体入水产生的涌浪破坏力巨大，由于波浪爬高及漫顶作用，会对库区挡水建筑物、水上航行、周边居民的生命财产造成严重威胁[1]。例如，1963年意大利的瓦依昂水库库岸发生滑坡，近3亿m3的滑坡体以30 m/s的速度滑入水库，掀起高达175 m的涌浪，摧毁了坝区的所有建筑物，66名水库技术人员全部遇难，涌浪以超出坝顶约100 m的高度冲入下游河谷，冲毁5个村镇，并夺去近3 000人的生命[2]。1961年我国发生了首例水库蓄水初期诱发的大型滑坡，湖南柘溪水库的大坝上游右岸1.5 km处的塘岩光发生了165万m3的大型滑坡，该滑坡以近20 m/s的速度入水后，产生的涌浪传播到对岸时浪高仍有21 m，并以3.6 m高的涌浪越过正在施工的坝顶，造成40余人死亡[3]。

为了研究滑坡涌浪的特性并对此类灾害进行预警及风险评估，学者们建立了多种不同的数学模型[4-5]。针对崩塌、巨石滚落以及历时较短、发展迅速的滑坡涌浪问题，很多学者将滑坡体简化为刚性体处理，即不考虑滑坡体在短时间运动过程中自身的变形特性。Wu等[6]利用大涡模拟模型以及VOF水面追踪方法，将物理模型试验中测得的滑块实际运动轨迹输入到数值模型中，研究了三维块体下滑造成的涌浪爬坡和回落现象；Heinrich[7]利用二维水动力NASA-VOF2D模型，计算了三角形滑块入水产生的波面，在整个模拟中利用已知的试验数据对滑块的运动参数进行了设定；Watts等[8]设定滑坡体位移-时间曲线为对数双曲余弦函数形式，利用Boussinesq方程模拟了滑坡体入水产生的涌浪传播过程。但是以上数学模型都是将滑坡体的运动作为已知参数直接输入到建立的水波模型中进行计算，没有考虑滑坡体入水后与水体的相互作用。也有许多学者直接假定滑坡体入水后作匀速运动或匀加速运动[9-11]，或假设滑坡体的位移-时间曲线为正弦或余弦形式[12-13]，或者假设水体对滑坡体的作用力与滑坡体速度满足一定的关系[14]。但是在实际情况中滑坡体的运动特性难以事先确定，特别是滑坡体入水后的速度变化规律更是难以预测，这样的假设往往会使模拟结果与滑坡涌浪的实际特性存在一定的偏差。

本文通过建立考虑滑坡涌浪过程中滑坡体与水体相互作用的流固耦合模型，利用FLOW-3D软件进行数值求解，研究刚性滑坡体入水生成涌浪的规模及其传播过程，并针对滑坡体在不同运动状态下所产生的涌浪特性进行对比，分析流场、涌浪爬高及其传播过程的异同。

1 数值模型

流体的控制方程采用考虑流固耦合作用的三维N-S方程：

(1)

(2)

式中：xi、xj(i，j=1,2,3)为坐标分量；ui为xi方向的流体速度分量；p为压力；νt为紊动黏性系数；gi为重力加速度在xi方向的分量；VF为计算域所占网格单元的体积分数，Ai为计算单元在xi方向的截面面积比例；ρ为流体的密度；t为时间。FLOW-3D采用FAVOR(fractional area/volume obstacle representation)方法定义计算单元并求得相应单元的VF与Ai的值[15]，随着滑坡体运动过程中位置和坐标的更新，流体计算单元的面积和体积分数也会被重新计算。

控制方程的空间离散采用三维矩形交错网格，数值离散采用有限差分法，同时采用VOF方法[16-17]对自由水面进行追踪处理。刻画自由水面的VOF函数满足方程：

(3)

式中：α为流体占计算单元空间的体积分数，取值范围为0～1。α=0说明计算单元内全为空气，α=1说明计算单元内全为水体，当0<α<1时，水气混合处于自由水面。式(3)中的右端项反映了滑坡体运动对流体空间的挤占作用。

采用k-ε双方程紊流模型模拟紊流效应，流体的紊动黏性系数νt=Cμk2/ε，其中经验系数Cμ取为0.09，紊动能k和紊动能耗散率ε的守恒方程分别表示为

(4)

(5)

式中：ν为水体的运动黏性系数；经验常数Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.0。

利用FLOW-3D中的通用刚体运动模型(general moving object, GMO)模拟滑坡体运动，可以选择指定滑坡体的运动轨迹或是通过流固耦合计算，当采用流固耦合模型时，滑坡体的控制方程可表示为

(6)

式中：Fi为滑坡体沿xi方向受到的合力，包括重力、流体的阻力和壁面摩擦阻力；m为滑坡体的质量；vi为滑坡体沿xi方向的运动速度。

2 模型验证

2.1 模型设置

采用建立的模型模拟滑坡体沿斜坡滑入水中产生的涌浪特性及其传播情况，为验证数值模型的有效性，将数值计算结果与Wu等[6]在港池中进行的滑坡涌浪试验观测数据进行对比。数值模型设置如图1所示，试验中的涌浪爬高测点和传播过程中的涌浪波高测点布置如图2所示，其中1号、2号用来测量涌浪沿斜坡的爬高值，3号、4号、5号、6号用来测量涌浪传播过程中不同位置处波高的变化情况。模型区域设置为6.6 m(长)×3.7 m(宽)×3.3 m(高)，3个方向的网格数量分别为60个、100个和60个，共计36万个网格；斜坡位于模型区域一端，长高比为2∶1。滑块剖面形状为三角形，相应尺寸为0.91 m(长)×0.61 m(宽)×0.455 m(高)，滑块密度设为2 640 kg/m3。初始水深设为3 m，初始时刻滑块下表面刚好与静止水面齐平，斜坡壁面摩擦系数设为0.2，模拟时间为4 s。

图1 模型设置示意图 (单位：m)

图2 涌浪爬高、波高测点布置(单位：m)

2.2 边界条件

2.3 结果分析

利用所建立的流固耦合数学模型进行数值求解，分别将波高测点和涌浪爬高测点的数值计算结果与物理试验结果进行对比，见图3和图4。

从图3中可以看出，数值计算结果与试验数据整体吻合较好，流固耦合模型能够较准确地模拟出滑坡入水后产生的首浪高度和涌浪传播速度；4号和6号测点处第二个波峰的计算高度略大于试验值，原因可能是由于模型中侧壁设置为对称反射边界，而试验中的港池侧壁不一定能够完全实现水波的反射，因此导致了波高计算值偏大。不过整体数值模拟结果还是能够较准确地反映涌浪的产生和传播过程。

从图4可以看出，模型计算所得涌浪爬高与试验观测值吻合较好，表明所建模型能够准确模拟滑坡入水后涌浪的传播与爬坡过程。

图3 不同测点涌浪波高计算值与试验值对比

图4 不同测点涌浪爬高计算值与试验值对比

3 不同滑坡体运动状态下的涌浪特性

3.1 匀加速运动

假设忽略水体对滑坡体的作用力和滑坡体所受摩擦力，即滑坡体在自身重力的作用下匀加速下滑，则滑坡体沿xi方向的加速度ai和速度vi可通过下式计算:

mgisinβ=mai

(7)

vi=ait

(8)

式中：β为斜坡倾角。

若不考虑滑坡体与水体的相互作用，将式(7)(8)得到的滑坡体运动轨迹直接作为输入参数，利用所建立的数学模型进行数值求解，3号、4号波高测点的数值计算结果与2.3小节中考虑流固耦合模型时的计算结果对比如图5所示。

图5 滑坡体匀加速运动与流固耦合作用下涌浪波高对比

从图5可以看出，若指定滑坡体在重力作用下作匀加速运动，则滑坡体入水产生的涌浪规模及其传播速度与利用流固耦合模型计算得到的结果存在较大的差异。采用指定滑坡体作匀加速运动的计算模型所得的首浪高度要比实际情况大30%左右，而且涌浪传播速度较流固耦合模型算出的结果也有明显的增加，但两种方法计算所得的整体波形基本一致，说明涌浪的产生和发展主要是由于滑坡体排开其运动轨迹前方一定体积的水而形成波峰，水体与滑坡体的相互作用只是改变了波峰的高度以及涌浪的传播速度，对整体波形影响较小。通过上述研究可以得知，在滑坡涌浪模拟过程中若直接指定滑坡体的运动轨迹进行计算所得到的结果是不合理的。

3.2 流固耦合+匀(加)速运动

为了研究流固耦合作用在滑坡体不同运动阶段对涌浪特性所起的作用，将滑坡涌浪过程分为两个阶段：滑坡体从开始下滑到完全浸入水中(入水阶段)和滑坡体完全浸入水体后(水下阶段)。入水阶段采用流固耦合模型，水下阶段忽略水体对滑坡体的阻力作用，采用给定滑坡体运动轨迹的方式进行计算。

首先利用流固耦合模型模拟计算滑坡体入水阶段(0≤t≤1.5 s)的涌浪特性，在模拟结束后，利用所得计算结果作为忽略水体阻力的水下阶段(1.5 s

图6 滑坡体分段运动状态下的涌浪波高对比

由图6可知，将滑坡体完全浸入水体时作为分界点分段计算所得的结果与全过程均用流固耦合模型计算所得的结果高度吻合，并且两种不同运动轨迹(匀速、匀加速)的计算结果也相差不大。表明当滑坡体完全浸入水体之后，其运动状态对涌浪的特性及传播过程几乎没有影响，说明滑坡体与水体的相互作用对涌浪特性的影响主要发生在滑坡体接触水面到完全浸没的过程中，因此在此阶段采用流固耦合模型来考虑滑坡体与水体相互作用对涌浪产生的影响是非常必要的。

4 结 论

a. 建立的滑坡涌浪三维流固耦合模型具有较高的精度，能够准确模拟出涌浪的波高、传播速度、沿岸爬高等特征参数。

b. 若不考虑滑坡体与水体的相互作用，直接给定滑坡体运动轨迹进行数值模拟所得的涌浪波高及传播速度要远远大于实际情况。

c. 滑坡体与水体相互作用对涌浪特性的影响主要体现在滑坡体入水阶段，在此阶段运用流固耦合模型进行数值求解，而在水下阶段，根据前一阶段流固耦合模型计算所得结果指定滑坡体入水后的运动轨迹所计算得到的结果与实际观测结果非常吻合。表明本文提出的滑坡体入水阶段采用流固耦合模型，水下阶段通过直接给定滑坡体运动轨迹的方法，既能较好地模拟出涌浪的产生与传播过程，又可以简化计算，对实际大型滑坡涌浪问题的模拟研究具有借鉴意义。