局部断面收缩的溃坝水流三维数值模拟

戎贵文，袁 岳，肖柏青，向英奇

(1．安徽理工大学地球与环境学院，安徽淮南 232001;2．水文水资源与水利工程科学国家重点实验室，江苏南京 210098)

病险水坝会对人民生命和财产安全构成潜在的威胁，溃坝洪水多年来一直是水利工程师和学者们重点关注的问题。溃坝水流具有突发性和灾害性，当溃坝发生时，准确预测溃坝洪水的演进过程对保护下游人民生命和财产安全具有重要意义。溃坝水流向下游演进过程中如果遇到变化的断面(如局部断面收缩)，水流流态会发生较大的变化。河道中的堤坝、树木、碎石以及其他一些建筑物都可以造成河道的局部断面收缩。在这种情况下，大坝下游各处随时间变化的水位、负波的形成与传播、正波和负波的速度等都是影响准确预测溃坝洪水演进过程的重要参数。

对于溃坝研究来说，由于难以获得随时间变化的速度、水位、流量以及波前、波峰的移动速度等水动力要素的现场测量数据，因此，人们通常在实验室中进行模型试验研究。Dressler［1］在水平矩形河道中对溃坝水流进行了试验研究。Bell等［2］研究了直线型河道和弯曲河道内的二维溃坝水流。Bellos等［3］在干-湿河床上试验研究了分汊河道内的溃坝水流。Lauber等［4］在矩形光滑的干水槽内进行了溃坝水流的试验研究，重点研究了溃坝水流正波和负波的形成和传播。虽然这些试验研究得到了一些关于溃坝水流的数据，但是由于试验环境的假定、试验模型的大小、测量装置的影响等使得试验结果的应用有一定的局限性。随着测量技术的提高以及计算机的高度发展，溃坝洪水研究结果变得更加准确和直观。Aureli等［5］采用试验手段和MacCormack计算格式数值模拟了局部收缩河道内的溃坝水流。戎贵文等［6］采用MacCormack预测-校正技术的隐式数值格式模拟了溃坝水流。符传君等［7］建立了溃坝水流运动的三维数学模型，使用k-ε模型对水流方程进行封闭，并对急弯河道进行了数值模拟计算。Liang［8］提出了一种求解不规则地形下的二维浅水水流数值模拟技术，模拟了局部收缩河道内的溃坝水流。王嘉松等［9］用TVD格式对瞬间全溃溃坝波的传播、反射和绕射进行了数值模拟。Chang等［10］采用光滑粒子流体动力学(SPH)方法求解浅水方程，研究了溃坝水流经过有三角形障碍物河道时的流动情况。Kocaman等［11］利用新的摄像技术试验研究了局部收缩对溃坝水流的影响。Lindsey等［12］利用试验和大涡模型(LES)数值模拟的方法研究了二维溃坝水流，得到了溃坝水流的一些水位变化剖面图。Jaan等［13］为研究溃坝水流遇到障碍物的运动过程，采用表面梯度逆风方法(SGUM)求解浅水方程(SWEs)，并与采用不可压缩光滑粒子流体动力学方法(ISPH)的模拟结果进行了比较。王晓玲等［14］建立耦合流体体积函数(VOF)法的三维标准k-ε紊流模型，模拟了三维溃坝洪水在复杂淹没区域的演进过程，并绘制了溃坝洪水风险图。虽然国内外对溃坝水流的研究已有很多，但是对三维局部断面收缩地形内溃坝水流流动的研究却很少，而且对该种地形下负波的形成与传播规律的研究更为罕见。

本文采用能够更好地处理高应变率及流线弯曲程度较大流动的RNG k-ε紊流模型［15-16］和具有倾斜网格校正能力、允许使用大时间步长、适合过渡流动计算的压力隐式算子分割(PISO)算法［17］，研究了三维局部断面收缩情况下溃坝水流的演进过程，分析了负波的形成和传播过程，得到了负波形成和传播的规律以及溃坝水流水位变化特征。

1 数学模型及数值求解方法

1．1 控制方程

式中:ρ为流体的密度，kg/m3;t为时间，s;xi、xj为坐标分量(i，j=1，2，3)，m;ui、uj为速度分量(i，j=1，2，3)，m/s;p为修正压力，Pa;μ为分子动力黏性系数，N·m/s;k为紊动能，m2/s2;ε为紊动耗散率，m2/s2;σk、σε分别为紊动能k和紊动耗散率ε对应的普朗特数;μt为紊动黏度;G 为紊动能生成项;cμ、C1ε、C2ε为经验常数，C1ε的表达式如下:

式中:Sij为平均应力张量。根据Yakhot等［15］的研究结果，σk、σε、cμ、C2ε的取值分别为 0. 7179、0. 7179、0.0845和1.68。

1．2 数值求解方法

采用有限体积法离散控制方程，对控制方程中的时间项使用隐式欧拉格式离散，对Navier-Stokes方程中的对流项采用TVD格式离散［18］，与此同时使用一阶迎风格式离散紊流模型控制方程中的对流项，其余各项采用二阶中心差分进行离散。物理变量采用交错网格布置，在每个控制单元中对微分方程进行积分，再将积分方程线性化，把控制方程离散为可以数值求解的代数方程，采用PISO算法对方程组进行求解［17］。该算法由 Issa［19］于 1986 年提出，适用于非稳态问题。PISO算法包括1个预测步骤和2个校正步骤，可以看作是在SIMPLE算法的基础上增加了1个校正步骤，通过一系列的预测-校正来求解每一时刻的流动物理量［20］。

1．3 自由表面的处理

溃坝水流自由表面采用VOF方法处理［21］。与其他处理自由表面的方法，如刚盖假定法、MAC法和LINC法相比，VOF方法用体积分数代替了标记点，只用一个函数就可以描述自由表面的各种复杂变化，大幅度节约了计算机内存和计算时间，而且计算精度也较高［16］。VOF方法的基本原理是通过研究网格单元中第q相流体体积和网格体积比函数αq来确定自由面，追踪流体的变化，而非追踪自由液面上质点的运动。若αq=1，说明该单元全部为第q相流体所占据;若αq=0，则该单元中没有第q相流体;若0＜αq＜1，表明该单元包含了第q相流体和一相或者其他多相流体的界面。ρ和μ由体积分数加权平均得出，计算公式如下:

式中:下标w和a分别代表水和气体;αw为水的体积分数。

2 数值模拟

2．1 物理模型

本文所采用物理模型是Kocaman等［11］的水槽模型(图1)，假定模型表面光滑。该模型整体上是一个六面体，长8.9 m，宽0.30 m，高0.34 m。用一块挡板将模型分为上下游两个部分，上游为有水区域，代表水库库区，下游为无水区域。初始时刻水槽上游水位为0.25 m，下游无水。在下游1.52 m处河道两侧有一个对称的三棱柱体，形成了一个局部断面收缩的区域，最窄处只有0.1 m。

图1 水槽模型示意图 (单位:cm)

2．2 网格划分及边界条件

2．2．1 网格划分

在笛卡儿坐标系下对上述物理模型进行网格剖分，x、y、z分别代表物理模型的长、高、宽方向，由于坐标原点的选择对文中的研究内容没有影响，理论上可以选择在任意位置，因此将坐标原点(0，0，0)选在挡板底部中心处。由于该模型是对称的，对其进行数值模拟计算时只需取计算模型的一半，所以在笛卡儿坐标系下只对模型的一半进行网格划分，划分的网格为均匀结构性网格，每个网格的大小为Δx=Δy=Δz=8 mm。

2．2．2 边界条件

由于模型是对称的，将模型对称面的边界条件设为镜面对称;模型的下游出口边界条件设为自由出流;由于上游没有来水条件，将上游边界设为壁面边界;该试验是在常规大气压下进行的，故将模型的顶部边界条件设为镜面对称;其他壁面的边界都设为壁面边界。

2．3 模型验证

在初始时刻，上游长4.65m、宽0.15m、高0.25m的区域内充满水，其他区域无水。在模型轴线处设置了6个测点，测点位置见图2。试验模拟了这6个测点处水位随时间的变化规律，通过与Kocaman等［11］的试验数据的比较，验证文中所采用模型的可行性。

图2 测量水位的位置点 (单位:cm)

图3显示了6个测点处水位随时间变化的模拟值与测量值，图3表明该模型可以准确地模拟局部断面收缩地形内的溃坝水流，验证了该模型的可行性。

2．4 模拟结果

2．4．1 负波的形成与传播

将挡板突然移去后，模拟的溃坝水流迅速流向下游。当溃坝水流到达局部收缩的断面时，一部分水流经过该收缩断面后继续向下游流动;另一部分水流拥堵在局部收缩段，抬高了这部分的水位，形成了一个波峰。当波峰达到一定高度时就会向上游移动形成负波，从图4可以清楚地看出形成负波的过程。随着挡板的突然移去，水流在t=0.9 s时到达收缩段，部分水流通过收缩段继续流向下游。在t=0～1.2s时，水流渐渐拥堵在收缩段中，水位不断升高。从t=1.5s开始收缩段上游形成了强烈的水跃，从收缩段到突扩段水位迅速降低。在t=2.1 s时，水流到达收缩段的最窄断面，此时水位迅速上升，部分水流反向上游流动。在t=3.0 s之前，收缩段的水位一直在升高，自由液面的形状为凸形。与此同时，流向突扩段的水流自由液面呈现出凹形。

图3 各测点处水位随时间变化的测量值与模拟值

图4 负波的形成与传播过程

溃坝水流由于通过局部收缩的断面，形成了负波。按时间顺序将负波的变化过程分为两个阶段，即负波的形成阶段和负波的传播阶段。从图4可以看出，t=0.9～3 s为负波的形成阶段，t=3～5 s为负波的传播阶段。笔者从能量的转化角度来解释负波的形成和传播过程。由于收缩断面的存在，当水流的动能不足以使水流通过断面时，水流就会拥堵在收缩口处将动能转化为势能。这个能量转化的过程导致了收缩口上游的水位逐渐升高，水流的流态也相应地发生变化，从上游的急流变为缓流。根据水力学的知识，在自然条件下当水流流态由急流转变为缓流时必然发生水跃现象，这就很好地解释了从t=1.5 s起开始形成水跃的原因。

图5 测点处水流速度随时间的变化过程

图5为测点处水流速度、负波速度随时间的变化过程，图5(a)反映了溃坝水流向下游传播时，上述6个测点处流速的变化情况，图5(b)反映了局部收缩断面处负波速度随时间的变化情况。由于取溃坝水流正常流动方向为正方向，所以向上游传播的负波速度为负值，图5(b)中以－ν波表示负波速度。从图5(a)可以看出，在t=1～3 s时，收缩断面处的测点P4、P5、P6处向下游流动的速度大幅度减小，这主要是由于上游大量水流到达最窄断面后水流不能及时流出，拥堵在收缩断面附近。图5(b)中，在t=2～3 s时，随着负波的移动，P4、P5、P6 处的速度增大，反映了负波的波前速度随着负波的移动逐渐增加。在t=3～5 s时，向下游流动的溃坝水流速度与向上游移动的负波速度逐渐保持不变，此时负波正向上游传播，并与向下游传播的水流处于能量平衡状态，这与图4中3～5 s时负波的自由液面渐渐呈现水平状态保持一致。

2．4．2 随时间变化的溃坝水流水位

图6为各测点处溃坝水流的水位变化过程。从图6(a)可以看出，随着挡板的突然移去，P1处的水位迅速降低，而随着水流的运动，P2处水位随后迅速升高，直到t=1.6s左右。在t=1.6～5.0s之间，P2处的水位基本保持恒定，水位约为0.116 m。Ritter［22］研究得出此处的水位 h=4/9h0(h0为来流水位)，P2处模拟得到的水位与Ritter的研究结果吻合。当t=5.0 s时，收缩段产生的负波影响了P2处的水位，P2处的水位随之迅速升高。同样的情形也可以从P1、P3处的水位变化看出。P1、P2、P3处的水位随时间的变化，间接地反映了负波的传播过程。从图6(b)中可以看出，t=0.9 s左右时水流到达收缩段的最窄部分，P4、P5、P6处的水位随后迅速上升。在t=3 s时P5处的水位达到最大值，约为0.20m，随后逐渐降低。图6(b)恰好说明了t=0～3 s为负波形成的时间，t=3 s后负波开始传播，P4、P5、P6处的水位保持在一定高度，直至水流越过收缩段最窄处，随后这3个测点处的水位逐渐降低。

图6 各测点处溃坝水流的水位变化过程

3 结论

本文利用三维RNG k-ε紊流模型对局部断面收缩情况下的三维溃坝水流进行了数值模拟，给出了水槽内负波的形成与传播过程、溃坝水流在测点处的水位和速度变化，分析了局部收缩对溃坝水流水力特性的影响。通过与已有试验数据的比较，验证了该模型的可行性，表明该模型可以准确地模拟局部断面收缩地形内的溃坝水流。采用该紊流模型得到了三维局部断面收缩水槽内溃坝水流的水位变化过程以及负波的形成和传播规律，可为研究实际情况下的溃坝水流提供参考。

但本文在研究溃坝水流时，假设大坝的溃坝形式为瞬间全溃，下游河道壁面光滑，同时溃坝水流中只含有气液两相。与全溃坝相比，部分溃坝液面变化复杂，并且上游形成的波形也不同;在相同的时间内溃口的大小对溃坝波推进的距离有影响，特别是溃坝初期溃口的大小对推进距离影响较大。实际的溃坝过程中，大坝并非瞬间全溃，而是部分逐渐溃决，而且由于水流的冲刷作用，水体中含有泥沙以及其他杂物，所以该模型还存在一定的局限性。为了更准确地模拟真实的溃坝水流，可以在模拟时加入泥沙，假设壁面有一定的粗糙度等，也可以将计算区域内的网格剖分得更精细，模拟的准确度将会更高，但计算时间和计算成本也会相应增加。

［1］DRESSLER R F．Comparison of theories and experiments for the hydraulic dam-break wave［J］．Int Ass Sci Hydrol Publi，1954，38(3):319-328．

［2］BELL S W，ELLIOT R C，CHAUDHRY M H．Experimental results of two dimensional dam-break flows［J］．Journal of Hydraulic Research，1992，30(2):225-252．

［3］BELLOS C V，SOULIS J V，SAKKAS J G．Experimental investigation of two dimensional dam-break induced flows［J］．Journal of Hydraulic Research，1992，30(1):47-63．

［4］LAUBER G，HAGER W H． Experiments to dam-break wave: horizontal channel ［J］． Journal of Hydraulic Research， 1998， 36 ( 3) : 291-308．

［5］AURELI F，MIGNOSA P，TOMIROTTI M．Numerical simulation and experimental verification of dam-break flows with shock［J］．Journal of Hydraulic Research，2000，38(3):197-206．

［6］戎贵文，魏文礼，严建军．溃坝洪水演进数值模拟研究［J］．黑龙江水专学报，2008，35(3):12-14．(RONG Guiwen，WEI Wenli，YAN Jianjun．Numerical simulation of 2D dam-break flood waves［J］．Journal of Heilongjiang HydraulicEngineering，2008，35(3):12-14．(in Chinese))

［7］符传君，练继建． 溃坝水流在复杂河道中传播的三维数 值模拟［J］． 水利学报，2007，38 ( 10) : 1151-1157． ( FU Chuanjun，LIAN Jijian． Three-dimensional numerical simulation of dam-break flow in complicated river sections［J］． Journal of Hydraulic Engineering，2007，38 ( 10 ) :1151-1157． ( in Chinese) )

［8］LIANG Qiuhua．Simulation of shallow flows in non-uniform open channels［J］．Journal of Fluid Mechanics，ASME，2008，130(1):51-59．

［9］王嘉松，倪汉根，金生．瞬间全溃溃坝波的传播、反射和绕射的数值模拟［J］．水动力学研究与进展:A辑，2000，15(1):1-7．(WANG Jiasong，NI Hangen，JIN Sheng．Simulation of propagation，reflection and diffraction of bores caused by sudden and full destruction of a dam［J］．Journal of Hydrodynamics:Ser A，2000，15(1):1-7．(in Chinese))

［10］CHANG Tsangjung，KAO Hongming，CHANG Kaohua，et al．Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics［J］．Journal of Hydrology，2011，408(1):78-90．

［11］KOCAMAN S，OZMEN C H．The effect of lateral channel contraction on dam break flows:laboratory experiment［J］．Journal of Hydrology，2012(432/433):145-153．

［12］LINDSEY A L，JASIM I M，HANIF C．Experimental and numerical investigations of two-dimensional dam-break flows［J］．Journal of Hydraulic Engineering，2012，139(6):569-579．

［13］JAAN H P，SHAO Songdong，HUANG Yuefei，et al． Evaluations of SWEs and SPH numerical modelling techniques for dam break flows ［J ］． Engineering Applications of Computational Fluid Mechanics，2013，7( 4) : 544-563．

［14］王晓玲，张爱丽，陈华鸿，等．三维溃坝洪水在复杂淹没区域演进的数值模拟［J］．水利学报，2012，43(9):1025-1033．(WANG Xiaolin，ZHANG Aili，CHEN Huahong，et al．Three-dimensional numerical simulation of dam-break flood routing in the complex inundation areas［J］．Journal of Hydraulic Engineering，2012，43(9):1025-1033．(in Chinese))．

［15］YAKHOT V，ORZAG S A．Renormalization group analysis of turbulence:basic theory［J］．J Scient Comput，1986，1(1):3-51．

［16］戎贵文，魏文礼，刘玉玲，等．涌潮作用下丁坝附近水流运动特性的数值模拟研究［J］．水利学报，2012，43(3):296-301．(RONG Guiwen，WEI Wenli，LIU Yulin，et al．Study on flow characteristics near spur dikes under tidal bore［J］．Journal of Hydraulic Engineering，2012，43(3):296-301．(in Chinese))

［17］魏文礼，王德意．计算水力学理论及应用［M］．西安:陕西科学技术出版社，2001．

［18］刘玉玲，王玲玲，周孝德，等．二维溃坝洪水波传播的高精度数值模拟［J］．自然灾害学报，2010，19(5):164-169．(LIU Yuling，WANG Lingling，ZHOU Xiaode，et al．High-resolution numerical simulation of 2D dam-break flood waves［J］．Journal of Natural Disasters，2010，19(5):164-169．(in Chinese))

［19］ISSA R I．Solution of the implicitly discretized fluid flow equations by operator-splitting［J］． Journal of Computational Physics，1986，62:40-65．

［20］李人宪．有限体积法基础［M］．北京:国防工业出版社，2008．

［21］HIRT C W，NICHOLS B D．Volume of fluid(VOF)method for the dynamics of free boundary［J］．Journal of Computational Physics，1981，39:201-225．

［22］RITTER A．The propagation of water waves［J］．Journal of German Engineers Association，1892，36(33):947-954．