邓成进 ，党发宁，陈兴周，袁秋霜，陈莉丽

(1. 西安理工大学 岩土工程研究所，陕西 西安 710048；2. 中国电建集团西北勘测设计研究院有限公司，陕西 西安 710065；3. 西安科技大学，陕西 西安 710065)

水库库区存在的滑坡、变形体、堆积体等不良地质体，在蓄水后的不利条件下可能失稳高速入水产生巨大的涌浪，尤其近坝库区滑坡涌浪会对大坝安全及下游构成巨大威胁。因此，坝前涌浪危害分析和研究一直是国内外学者关注的重点问题。

库区滑坡涌浪分析方法主要有经验公式、模型试验和数值分析等。常用的经验公式[1]有潘家铮法、水科院法、美国土木工程学会法等，经验公式法计算简便、应用较为广泛；而物理模型试验的成本较大、试验周期长，应用相对较少。随着计算机技术的发展，基于流体动力学理论的数值分析方法发展成熟，数值分析将逐渐成为解决此类问题的重要手段。周桂云等[2]基于浅水控制方程通过自编程序GEO-FLOW实现了滑坡涌浪传播动态过程的模拟；邓成进等[3]基于FLUENT软件采用有限体积法模拟了二维滑坡涌浪过程。除了传统的有限差分法和有限元法，还有可模拟涌浪产生过程的SPH方法[4-5]和基于分步有限元法的LEVEL SET数值模拟方法等[6]。Viroulet等[7]对比分析了SPH法和有限体积法模拟二维楔形滑坡涌浪规律，并采用模型试验验证。任坤杰等[8]采用非规则网格的流动标点法建立二维模型计算广家崖松散堆积体滑坡涌浪。上述涌浪数值分析主要集中在二维模型或小范围的三维模型，并与物理模型试验取得相对一致的结果，验证了基于Navier-Stokes方程并引入湍流模型和自由表面处理技术的可行性。但大尺寸复杂的三维数值分析往往需要足够多的网格和时间步，需付出巨大的计算代价及成本，因此考虑复杂的库区地形条件高坝库区涌浪的研究和应用较少，大多仅作为辅助研究手段[9]。

本文对某水库3#变形体整体滑动激起的涌浪进行三维数值分析，分析网格尺寸对计算精度的影响规律，采用水工物理模型试验进行验证，并与常规潘家铮法计算的涌浪高度进行对比分析，研究涌浪对大坝安全的影响，为大坝安全评价提供技术支撑。

1 某库区滑坡概况

某库区3#变形体在大坝上游约1 300 m，岩性为花岗岩体。水库库岸冲沟、山梁相间。3#变形体边坡位于两冲沟之间，冲沟底部有基岩出露，边坡表部分布有松动破碎岩体和崩坡堆积体。岩体以倾倒崩塌破坏为主，根据边坡形成拉张裂缝及结构面特征，总方量约360万m3。

3#变形体内存在的陡倾岸内结构面（图1）是其发生倾倒变形的决定性因素，并与发育的陡倾岸外的结构面相互切割。水库蓄水后原有岩体的平衡被改变，导致水下岩体发生倾倒变形，下部较小的变形传递至中上部岩体使之发生大的倾倒变形，最终引起顶部产生“倾倒-拉开-塌陷”变形，累计最大综合位移达40.2 m。目前分析认为，3#变形体以倾倒崩塌破坏为主，不存在沿深部特定结构面产生整体滑动的可能。

3#变形体岩体倾倒变形规模和量级大，且变形条件及成因机制复杂，一旦失稳下滑，产生的涌浪将严重危及大坝安全。因此，需对其在极端情况下发生整体滑动的涌浪问题进行专门研究。

图1 3#变形体岩体结构剖面Fig. 1 Rock structure of 3# deformation body

2 库区滑坡涌浪三维数值模拟

2.1 理论基础

库区滑坡涌浪模拟基于流体动力学的研究，将水流运动视为不可压缩的黏性流体运动，采用离散完整的Navier-Stokes方程，网格中的障碍物采用体积分数和面积分数来定义[10-11]，矩形网格能够精确地模拟复杂的几何形状，表达式如下：

式中：Ax、Ay、Az分别为x、y、z这3个方向可流动的面积分数；u、v、w为3个方向的速度分量；ρ为流体密度，本文中即为水的密度；VF为可流动的体积分数，自由液面的体积分数即为0.5；Gx、Gy、Gz为3个方向的重力加速度；fx、fy、fz为3个方向的黏滞力加速度；p为水的压强。

FLOW3D软件采用有限差分法进行数值离散，将模型离散成空间矩形网格，速度和面积分数参数位于矩形网格边界面的中心点上，求解方法采用广义的极小残差算法（GMRES算法）。

常用的紊流模型有Reynolds应力模型、标准k-ε模型和RNG k-ε紊流模型等，RNG k-ε湍流模型可以很好地模拟高应变率和流线弯曲程度较大的流动[12]。滑体入库及库区内涌浪相互作用时会导致水面出现巨大的变形、破碎，计算模型选用RNG k-ε湍流模型，方程表达式如下：

式中：kt为紊动能量；Gt为由浮力引起的紊动动能产生项；Pt为由速度梯度引起的紊动动能产生项；εt为紊动耗散率；Dt、Dε为紊动扩散项；C1、C2、C3为无量纲的经验系数。

计算采用VOF法追踪自由水面[13]，引入体积函数αq的概念实现对自由面的追踪：空气体积函数αq=0，流体内部体积函数αq=1.0，自由液面的体积函数为αq=0.5。

2.2 计算模型及初始化

根据现场地形资料将滑床简化为59°斜坡，水库库岸三维模型由多个河谷剖面相接而成，挡水建筑物为双曲拱坝，大坝坝顶高程2 460 m，最大坝高250 m，水库正常水位为2 452 m，模型不考虑坝身孔洞、防浪墙等细部结构。根据3#变形体一次落水方量，将滑体简化为矩形刚性滑体。滑体顺河向宽220 m、厚度为26 m，落水点距对岸的水面宽度约为530 m。库区地形三维模型、挡水建筑物及滑体模型均采用犀牛软件完成，并导入FLOW3D计算。三维模型上游、左右岸及下游大坝采用Wall边界条件，水在边界上速度为零；上部自由水面边界满足自由面法向速度为零，还满足水体与气体交界面应力平衡条件。模型初始化时自由水面高程为2 452 m，模拟计算时间步△t=1.0 s，总计算时间步为160 s。水库库岸涌浪三维计算模型见图2。

3#变形体倾倒变形破坏机理较为复杂，且破坏运动及堆积形态受河谷地形、滑体与水阻力的影响，数值模拟及模型试验均很难精确实现其运动及堆积过程。本文采用常用的潘家铮法估算3#变形体整体滑动速度，3#变形体整体滑动最大滑速为24.6 m/s。为了考虑滑体解体后滑体堆积对水体的抬升作用，滑体滑至库底后，采用一梯形刚体块体由库底向上运动，以模拟不同时刻水库中岩土体的堆积形态(图3)。

图2 水库库岸涌浪三维计算模型Fig. 2 Three-dimensional model of reservoir surge

图3 水库滑体及滑体堆积过程模拟Fig. 3 Simulation of reservoir sliding body and sliding body accumulation

3 库区滑坡涌浪数值分析计算精度

3.1 网格尺寸对计算精度影响分析

首先通过二维模型分析网格尺寸对计算精度的影响，计算模型网格尺寸分别为5.0、2.5、1.0和0.5 m，其他计算条件一致。计算结果表明，当涌浪形成时，水体受滑体的挤压，水面形成巨大的浪花，不同模型计算的涌浪形成过程基本一致，表明网格尺寸对模拟库区涌浪变化过程影响较小。在t=5 s时水面达到最高，模拟的涌浪形态及涌浪高度见图4，其中涌浪高度（ηmax为最大涌浪高度）为水面超过初始水面高程2 452 m以上的高度。

图4 不同网格尺寸下二维模型计算的涌浪形态及涌浪高度Fig. 4 Calculation of initial surge shape and surge height by two-dimensional model with different grid sizes

由图4可知：不同模型模拟计算的浪花翻滚形态不同，当网格间距为2.5、1.0和0.5 m时，可以清晰模拟出滑坡滑动激起的水面翻滚、破碎；而网格间距为5.0 m时，无法模拟出涌浪翻滚的形态。不同模型计算得到的浪高差别较大，模型网格尺寸越小，越能清晰模拟涌浪浪花翻滚、破碎的微观形态，计算的涌浪高度值越大；当网格间距分别为2.5、1.0和0.5 m时，最大涌浪高度分别为35.97、48.28和54.50 m。上述研究表明，网格尺寸对数值分析的精度影响较大，网格尺寸越小，数值分析计算精度越高，但计算代价越大。

考虑大尺寸三维数值分析计算成本和计算收敛性，以及网格尺寸对数值分析精度的影响，三维模型采用网格间距分别为5.0和2.5 m，在t=5 s时初始涌浪形态见图5所示。

图5 不同网格尺寸下三维模型计算的初始涌浪形态Fig. 5 Calculation of initial surge shape by three-dimensional model with different grid sizes

由图5可知：三维模型网格间距为5.0 m时，无法模拟出涌浪翻滚的形态；而网格间距为2.5 m时，可模拟出水面翻滚的形态，初始涌浪高度值较大。将二维、三维计算结果进行对比可知，三维计算的初始涌浪高度比二维计算的略小。随着涌浪的传播，三维、二维模拟计算的涌浪传播过程有着明显的差异，二维模型不能直观展示涌浪传播过程。

3.2 涌浪传播过程计算精度影响规律

根据三维数值分析成果，选择落水点对岸、落水点沿岸距离330 m及坝前位置的涌浪高程变化过程，分析涌浪传播过程中网格尺寸对计算精度的影响程度。三维模型网格间距2.5和5.0 m计算涌浪高度变化曲线见图6。

从图6可知：不同模型计算的涌浪高度变化趋势基本一致，但各波峰值有一定差异，模型网格间距2.5 m时计算的波峰值大。在落水点对岸处，网格间距2.5 m计算的首浪波峰值比网格间距5.0 m的大10.5%；在距落水点沿岸330 m处，网格间距2.5 m的仅大5.8%；在坝前处，两者计算结果基本一致。落水点对岸和沿岸330 m处的第2次波峰、坝前第4、5次波峰时，两种模型计算的波峰相差较大，此时涌浪波动频率相对较快，水面变化相对剧烈。

图6 三维模型不同网格尺寸计算涌浪高度变化曲线Fig. 6 Surge height curves calculated by 3D model with different grid sizes

上述分析表明，数值分析的计算精度不仅与模型网格尺寸有关，还与涌浪传播距离、涌浪剧烈程度相关。随着涌浪传播，离落水点的距离越远，模型网格尺寸的影响越小，计算精度越高，当涌浪传播至坝前时，计算精度基本满足要求；当水面波动剧烈时，网格尺寸的影响增大，计算精度会有所减小。

4 库区滑坡涌浪水工物理模型试验验证

为分析和研究滑坡涌浪对大坝安全的影响，建立了3#变形体库区涌浪水工物理试验模型，试验模型几何比尺为1∶280。试验水位为水库正常蓄水位2 452 m，通过将滑块放至不同高度下滑以获得相应入水速度，涌浪高度采用DJ800型水工数据采集系统收集。为方便与数值分析计算成果进行验证，将落水点对岸及坝前位置的涌浪高度变化与三维模型数值分析计算结果进行对比(图7)。

图7 各监测点数值模拟和模型试验涌浪高度对比Fig. 7 Comparison of surge heights between numerical simulation and model test at each monitoring point

由图7可知：水工物理模型试验和三维数值分析得到的涌浪高度变化曲线的整体趋势基本一致，数值分析可展示滑坡涌浪传播整个宏观过程。两种方法得到水面整体变化趋势的一致性，验证和保证了数值分析计算的合理性，但两者计算涌浪形成过程仍存在一定误差。

水工物理模型试验得到波浪峰值比数值模拟结果大，数值上存在一定误差：在落水点附近处，两种方法计算的涌浪高度误差较大；传播至坝前时，两者计算的误差减小。数值上的误差主要是由网格划分精细程度引起的，涌浪微观形态模拟需更小的网格，三维数值分析的网格很难达到一定的精细程度，故数值分析计算结果小。但随着涌浪传播过程中能量的损耗，浪花的翻滚能量减小，网格尺寸对涌浪高度的影响逐渐减小，离落水点的距离越远，涌浪高度计算的误差越小。

上述分析表明，当涌浪传播至坝前时，涌浪波动剧烈程度已大幅减弱，网格尺寸对计算精度影响减小，可满足坝前涌浪高度的计算精度要求，可为库区滑坡涌浪灾害预测及其大坝安全研究提供完整的数据。

5 库区涌浪高度计算的讨论和分析

5.1 库区滑坡涌浪高度分析

下面采用常用的潘家铮经验公式计算库区涌浪高度，并与三维数值分析、水工物理模型试验得到的涌浪高度进行对比。潘家铮法计算时，取库水深h=145 m，滑体平均厚度λ=26 m，系数m=1.17，河谷宽取平均值B=530 m，滑坡平均宽度L=220 m。潘家铮经验法建立了一系列的假定，无法考虑挡水建筑物对涌浪的影响[14]。同时为了与潘家铮法计算结果进行对比，采用无挡水建筑物的三维模型进行数值分析，得到坝址处涌浪高度仅为4.29 m。采用多种方法计算得到的涌浪高度见表1。

表1 多种方法得到的库区涌浪高度Tab. 1 Surge heights calculated by various calculation methods in the reservoir area 单位：m

由表1可知：水工物理模型试验得到涌浪高度最大，数值分析计算的涌浪高度次之，潘家铮法估算结果最小。在落水点对岸，多种方法得到的涌浪高度差别较大。涌浪传播至坝前时，各方法得到涌浪高度的差别均有所减小。这是由于模型试验和数值模拟计算的最大浪高为浪花翻滚的最高点，落水点附近的涌浪高度计算结果相差较大。当涌浪传播至坝前时，浪花翻滚的高度减小，此时各方法计算误差减小，计算结果也较为接近。但近坝库区滑坡涌浪传播过程受挡水建筑物的影响较大，常规的潘家铮公式无法考虑挡水建筑物的壅水作用，因此得到的坝前涌浪高度与水工物理模型试验及数值分析的差别较大，不能作为大坝安全评价依据。

5.2 库区滑坡涌浪对大坝安全影响分析

滑坡涌浪作用在大坝面上力的大小与坝前涌浪形态密切相关[15]。由图7（b）可知，当涌浪传播至坝前时，在t=53、86、131和151 s时共有4次波峰超过坝顶，其中第4次波峰高度达到最大，为16.37 m，超过坝顶约8.37 m，此时数值分析获得的水面形态见图8所示。由图8可知，3#变形体库区滑坡涌浪在两岸坝肩处有小范围漫坝，每次漫坝持续时间较短，过流量有限，考虑拱坝超载能力较强，因此产生的涌浪不会对大坝及下游安全构成威胁。

图8 第4次波峰时涌浪漫坝的形态Fig. 8 Shape of surge over the dam at the 4th wave peak

6 结 语

（1）滑坡涌浪数值分析的模型网格尺寸越小，越能清晰模拟涌浪微观形态，数值分析计算精度越高。同时，计算精度还与涌浪传播距离、涌浪剧烈程度密切相关，当涌浪传播至坝前时，涌浪波动剧烈程度大幅减弱，网格尺寸对计算精度影响也减小。

（2）模型试验和三维数值模拟的涌浪高度变化曲线的整体趋势基本一致，数值分析可展示滑坡涌浪传播整个宏观过程。随着涌浪传播，离落水点的距离越远，两者得到涌浪高度差别也越小，坝前涌浪计算基本满足精度要求，可作为高坝大库下游涌浪灾害预测和预警的依据。由于潘家铮法无法考虑挡水建筑物壅水作用过程，计算的坝前涌浪高度与数值分析、模型试验相差较大，不能作为大坝安全评价依据。

（3）坝前最大涌浪发生在左岸坝肩位置，在53 s时坝前的首浪高度达到10.86 m，由于库区水面的反复震荡、叠加作用，在151 s时第4次波峰造成坝前涌浪高度达到16.37 m，而后逐渐衰减。坝前涌浪高度超过坝顶8.37 m，仅占最大坝高的3.3%。拱坝超载能力较强，所以产生的涌浪不会对大坝安全构成威胁。