多孔溢洪道水力特性的三维数值模拟

宁利中, 宁碧波, 宁景昊, 田伟利

(1.西安理工大学 省部共建西北旱区生态水利国家重点实验室，西安 710048;2.嘉兴学院 建筑工程学院, 浙江 嘉兴 314001；3.上海大学 建筑系, 上海 200444)

0 引 言

溢洪道是水利工程中的重要泄水建筑物。其作用是渲泄超过水库调蓄能力的洪水，确保水利工程的安全[1]。受地形，水利枢纽布置的限制，溢洪道可分为正槽式溢洪道、侧槽式溢洪道、井式溢洪道等。正槽式溢洪道是经常采用的一种泄流设施。关于正槽式溢洪道或者溢流坝水力特性的研究已经有了一些研究成果[2-6]。溢洪道水力特性的研究一般分为数值模拟方法和模型试验方法。数值模拟具有成本较低，模拟运算周期短等优点。现在，数值模拟方法广泛应用于水利工程泄水建筑物水力特性的研究中。

溢洪道分为单孔和多孔布置形式。多孔溢洪道的水力特性较单孔溢洪道复杂，计算工作量较大。文献[7-8]采用标准k-ε紊流模型对溢洪道水力特性进行了三维数值模拟，但未考虑门槽、闸墩的影响及多孔泄流的特性。文献[9-10]模拟了门槽和闸墩对水力特性的影响,发现溢洪道不同孔泄流时，水力特性不同。但未就多孔泄流时闸墩后水流的汇合问题进行研究。在模拟自由水面时，文献[11]在溢洪道水流计算中,文献[12]在对张峰水库溢洪道水力特性的模拟中，文献[13]对“S”型溢洪道流动特性的数值模拟中，文献[14]在溢洪道两相流水流模拟中，都采用了VOF法。文献[15-16]对泄水建筑物水力计算进行了探讨。文献[17-19]对溢洪道流动特性及掺气槽等进行了研究。因此，泄水建筑物的水力计算与数值模拟是一个值得深入研究的问题。本文以某工程3孔溢洪道为例，使用Gambit对模型进行网格划分，紊流时均方程使用RNGk-ε紊流模型封闭，自由表面采用VOF算法模拟，采用PISO算法进行压力速度耦合方程求解。通过Fluent软件对多孔溢洪道水流特性及水流在闸墩后泄槽内汇合与扩散情况等主要水力特性进行了数值模拟。模拟得出3种工况下溢洪道水面线、底流速的变化，速度分布，压强等水力参数。其结果与模型试验数据吻合一致。通过对溢洪道流场分析，进一步对溢洪道闸墩后多孔水流汇合与扩散边界问题进行了研究。

1 工程概况

某水利枢纽的泄水建筑物由布置在河床中部的3孔坝身表孔溢洪道组成，溢流坝段宽79 m，堰顶高程561.00 m，单个孔口尺寸为15 m×18 m，闸墩宽4.5 m。溢流堰是堰面曲线为y=0.053 71x1.85的WES实用堰，泄槽宽度为54 m，坡度为1∶0.75。泄槽后接挑流鼻坎，反弧半径为30 m，挑角为18.19°，挑流鼻坎顶高程481.00 m。溢洪道布置与体形见图1。试验模型按重力相似设计，几何比尺Lr=80。

图1 多孔溢洪道布置与体形Fig.1 Layout and shape of spillway with multi-opening

2 控制方程与数值模拟方法

2.1 控制方程

考虑水流是不可压缩的，采用RNGk-ε湍流模型，方程为[20]

连续方程：

(1)

动量方程：

(2)

雷诺应力方程：

(3)

式中：t为时间；xi为沿i方向空间上的坐标(i=1,2,3)；ui为沿i方向上时均速度分量；μt为湍动黏度；μ和ρ分别为流体分子黏性系数和密度；p为压强；gi为i方向上的重力加速度；k为湍动动能；δij为Krone-cker函数：δij=1(i=j)，δij=0(i≠j)。

k方程 ：

(4)

ε方程：

(5)

2.2 数值方法

运用GAMBIT对模型计算区域采用了有限体积法离散，采用PISO算法进行压力速度耦合求解，PISO执行了相邻校正和偏斜校正，校正后的速度比SIMPLE更满足连续方程和动量方程。PISO算法在每次迭代上花费稍多的时间，但最大程度地减少了迭代次数，从而提高了计算效率。基于 VOF(The Volume of Fluid)模型捕捉自由表面。时间步长为0.001 s。以进出口断面流量差低于总流量的5%为判断收敛的依据。

2.3 模型计算区域的网格划分及计算工况

数值模拟以某水利枢纽溢洪道尺寸为计算区域，建立了溢洪道三维模型，以堰顶为坐标原点O。桩号从堰头算起，坐标原点的桩号为3.88 m。确定溢洪道上游水库长度是5倍的堰上水头，长度为90 m。上游水库和溢洪道泄槽宽度为54 m。鼻坎后为下游河道，长度为300 m，宽度为74 m。计算区域全长504.8 m，宽54～74 m，高146 m，为了满足计算精度并节省计算时间，上游库区和下游河道采用六面体网格划分，闸室段、泄槽段和挑坎段选用混合网格，并进行加密处理。全部网格数量约为68万。

数值模拟采用了3种泄流工况：①设计水位：水位578.95 m，泄流量7 228 m3·s-1，3孔全开；②校核水位：水位579.75 m，泄流量5 200 m3·s-1，左孔，中孔全开；(3)500年一遇洪水：水位571.80 m，泄流量3 230 m3·s-1，3孔全开。

2.4 边界条件

①进口边界，溢洪道计算模型分为上部空气进口和下部水流进口。计算模型的水体自由表面与大气联通。水流进口采用速度进口，取来流方向为X正向。空气进口选用压力进口，压力值为一个大气压。②出口边界，下游出口边界为压力出口。③固壁边界，溢洪坝底板及侧壁均选用无滑移边界条件。

计算初始条件为上游水库以计算水位充满水，计算水流从堰顶向下游流动。

3 模拟结果分析

3.1 水流流态

溢洪道溢流堰正对水流流动方向，溢流堰前库区水流平稳，来流顺畅。设计水位时，溢洪道3孔弧门全开，溢流堰前库区水流顺畅，由于闸墩墩头的影响，在闸墩墩头处形成绕流。挑流水舌落点较为集中。校核洪水时，中孔、左孔开启，形成不对称泄流，这时泄槽内形成不对称折冲水流，泄槽内形成局部水面抬高，不同断面的水面变化较大。水流经挑流鼻坎平顺挑向下游，水舌有大的横向扩散，挑距也较大。

3.2 水面线变化

数值模拟获得了3种运行水位情况下溢洪道横断面中心线处的纵向水面线变化。以水气交界面处含气浓度50%为标准确定水面线的位置。设计条件下模型试验获得的溢洪道中心线上水面线变化与数值模拟结果见图2。溢洪道溢流堰抛物线段桩号0-04.00时，计算获得的水位高程是573.8 m。桩号0+10.00时计算获得的水位高程是570.68 m。桩号0+30.00时计算获得的水位高程是549.38 m，溢洪道溢流堰抛物线段水位高程沿程明显降低。说明闸室中水流明显跌落。桩号0+30.00时试验获得水深12.5 m，数值模拟获得水深为12.3 m，模拟值小于试验值1.6%。水深在泄槽内沿程逐渐减小，数值模拟获得的平均水深为5.72～4.5 m。挑流鼻坎末端附近水深约3.0 m。桩号0+99.129时试验获得水深为3.5 m，数值模拟获得水深为3.9 m，数值模拟的水深大于试验值4%。表明溢洪道的沿程水深在闸室段跌落明显，在泄槽段水流平稳，水深缓慢降低。数值模拟的水深与试验的水深的误差小于4%，说明两者吻合一致。

3.3 临底流速变化

数值模拟获得了3种运行水位情况下沿溢洪道横断面中心线处的纵向临底流速变化。设计条件下模型试验获得的溢洪道中心线上临底流速的变化与数值模拟结果见图3。由图3可见，数值模拟获得的临底流速与试验观测获得的临底流速的沿程变化趋势一致，吻合良好。水流进入溢洪道溢流堰后水面明显跌落，流速逐渐增大。在溢洪道溢流堰堰顶附近，临底流速达到10.5 m·s-1，然后，临底流速逐渐增大，到反弧段入口附近临底流速为35.5 m·s-1左右，到反弧段最低点临底流速为37.9 m·s-1，随后到反弧段挑坎末端前，临底流速有所降低。

图2 设计水位情况下沿溢洪道中心剖面的水面线Fig.2 Water surface line along the center section of spillway under design water level

图3 设计水位下沿溢洪道中心剖面的临底流速Fig.3 Bottom velocity along the center section of spillway under design water level

图4 流速垂向分布Fig.4 Vertical distribution of velocity

3.4 流速垂向分布

在3种水位情况下，模型试验时沿溢洪道的中轴线选择典型断面观测了流速垂向分布，数值模型也可以获得相应断面的流速垂向分布。设计水位下数值模拟获得的典型断面的流速垂向分布与实测资料的比较见图4。图4(a)为桩号0+017.10 m断面，模型试验实测的临底流速为16.4 m·s-1，数值模拟值获得的临底流速为17.1 m·s-1。模型试验实测的临底流速与模拟结果误差小于4%。从流速垂向分布的对比可见，模型试验实测的表流速与模拟的表流速均略小于流速垂向分布最大值。图4(b)为桩号0+040.00 m断面，模型试验实测的临底流速为27.69 m·s-1，表面流速增至28.36 m·s-1。数值模拟值获得的临底流速为28.4 m·s-1，表面流速增大到29.6 m·s-1。模型试验实测的流速与模拟结果吻合一致。从流速垂向分布的变化可见，在壁面附近流速由零迅速最大，急剧变化，存在边界层。在边界层厚度以外，流速沿垂向分布基本为常数。由图4(a)与图4(b)的比较可见，边界层厚度以外的流速随着流程的增加而增大。

图5 设计水位下沿溢洪道中心断面的底板压强Fig.5 Floor pressure along the central section of spillway under design water level

3.5 压强变化

数值模拟获得了不同运行水位下沿溢洪道横断面中心线Y=-27 m处的底板压强变化。设计条件下模型试验获得的溢洪道中心线上底板压强变化与数值模拟结果见图5。由图5可见，在WES堰顶上游压强较大。随后在溢洪道溢流堰顶附近桩号0～0+013.0 m之间出现负压区，至抛物线顶端降至最低，但负压不大。之后负压有所减小。到达泄槽直线段后压强逐渐增加。在反弧段，由于离心力作用引起动水压强显著增加，在反弧最低点压强达到了最大值，其数值为180 kPa左右。在反弧段末端，由于挑流出口的影响，压强逐渐减小。数值模拟获得的溢洪道横断面中心线处的底板压强与模型试验数据的变化趋势基本一致，并且两者吻合良好。说明溢洪道的体形设计合理。

3.6 多孔溢洪道闸墩后水流的汇合问题

图6 设计水位下不同断面处表流速分布Fig.6 Velocity distribution near the water surface in different section under design water level

图7 水流汇合长度与佛氏数的关系Fig.7 Relation between flow convergence length and Froude number

数值模拟获得了不同运行水位不同断面处的表流速分布。设计水位下不同断面处的表流速分布见图6。图6(a)为桩号0+023.90 m断面的表流速分布，左孔、中孔及右孔的流速分布比较均匀，最大表流速19.5 m·s-1左右。由于该断面位于闸室中，横坐标从15～19.5 m及34.5～39 m是闸墩厚度，没有流动，相应的流速为0。闸室闸墩末端断面桩号0+035.00 m，图6(b)为桩号0+043.90 m断面的表流速分布，左孔、中孔及右孔主流的表流速分布比较均匀，最大表流速31.0 m·s-1左右。由于该断面已经离开闸室，水流开始向横坐标从15～19.5 m及34.5～39 m的闸墩原有位置汇合，汇合处表流速低于左孔、中孔及右孔主流的表流速，其表流速为5.5 m·s-1左右。总表流速分布形成一个“山”形结构，3个主流区流速较大，存在两个低流速区。图6(c)为桩号0+048.90 m断面的表流速分布，左孔、中孔及右孔主流的表流速分布比较均匀，最大表流速34.0 m·s-1左右。在横坐标从15～19.5 m及34.5～39 m的位置水流已经汇合，汇合处表流速较图6(b)有明显提高，其表流速达20.0 m·s-1左右，但仍然低于左孔、中孔及右孔主流的表流速。总表流速分布仍然呈“山”字形结构。图6(d)为桩号0+058.90 m断面的表流速分布，表流速分布仍然呈“山”字形结构且大为改观。左孔、中孔及右孔主流的最大表流速37.0 m·s-1左右。在横坐标从15～19.5 m及34.5～39 m的汇合处表流速接近主流的最大表流速，其表流速达33.0 m·s-1左右，但表流速横向分布仍然不均匀。到图6(e)桩号0+063.90 m断面时，表流速横向分布基本均匀，沿泄槽宽度方向基本形成矩形分布,最大表流速40.0 m·s-1左右。如果以泄槽内闸墩后水流的表流速变均匀作为多孔溢洪道闸墩后水流汇合完成的标准。设计水位下水流从闸墩末端到3孔汇合完成沿泄槽底面经过了约48.1 m。同理，对校核洪水位及500年一遇洪水作出表流速分布图，可以获得相应的水流从闸墩末端到汇合完成的长度。

由于水流佛氏数越大，水流急流越快，需要的水流汇合长度越长。因此，选择佛氏数作为变量，分析其与闸墩后水流汇合的相对长度的关系。3种工况下闸墩后水流汇合的相对长度与水流佛氏数的关系见图7，虽然只有3个点，但规律性良好，可见汇合相对长度随着佛氏数的变化。可给出水流从多孔溢洪道闸墩末端到汇合完成的长度的经验公式为

(6)

式中：Lx为从闸墩末端到水流汇合完成的长度;h为水流汇合处的水深;q为单宽流量;g为重力加速度。

3.7 多孔溢洪道闸墩后水流的扩散边界

图8 校核水位下不同断面的表流速分布和流动扩散Fig.8 Velocity distribution near the water surface and flow divergence in different sections under check water level

图9 流动单侧扩散宽度Fig.9 Unilateral divergence width of flow

在校核洪水位情况下左孔和中孔全开，水流从闸室到泄槽不同断面表流速分布见图8。来自左孔和中孔的水流汇合过程与图6的设计水位情况时水流的汇合过程类似。离开闸室末端的水流沿泄槽底面约56.6 m左右汇合完成。水流在泄槽中的横向扩散可由离开中孔的水流右流速边界的变化见图8(a)～(f)。若以流速为零的点作为水流扩散右边界的标准，就可以确定不同断面的水流扩散右边界的位置。图8(a)桩号0+023.9 m位于闸室中。图8(b)桩号0+043.90 m处水流刚离开闸墩末端不远，中孔右侧流速最大值到零的变化比较陡，水流向右的扩散距离较小。随着桩号的增加，流程变长，中孔右侧流速最大值到零的变化比较缓，水流向右的扩散距离变大。经分析，水流随着流程的扩散角大约为15.94°。与文献[1,21]中的结果类似。若以水流单侧横向扩散宽度Δb/H为纵坐标，以水流离开闸墩末端沿泄槽底板的距离L/H为横坐标，可以给出中孔水流单侧横向扩散的变化，见图9。经分析，中孔水流单侧横向扩散的经验方程为

(7)

式中：H为堰上水头；Δb为单侧横向扩散宽度；L为水流离开闸墩末端沿泄槽底板的距离。

4 结 论

使用Gambit对模型进行网格划分，紊流时均方程使用RNGk-ε紊流模型封闭，自由表面采用VOF算法追踪。通过数值模拟研究了某工程溢洪道的水力特性。可以得出如下结论。

1)数值模拟获得了溢洪道水面线、底流速的变化，速度分布，压强等水力参数，它们分布合理，并与模型试验数据吻合一致。说明本文数值模拟方法可以模拟多孔溢洪道的水力特性。多孔溢洪道的体形设计合理。数值模拟方法与模型试验相结合是水工水力学研究的一个有效途径。

2)通过对溢洪道表流速分布的分析，以泄槽内闸墩后水流的表流速变均匀作为多孔溢洪道闸墩后水流汇合完成的标准。给出了多孔溢洪道闸墩后水流汇合完成长度的经验公式。

3)研究了泄槽内多孔溢洪道闸墩后水流的单侧横向扩散问题。以流速为零的点作为水流扩散右边界的标准。给出了水流单侧横向扩散的扩散角及水流单侧横向扩散宽度的经验关系式。

4)本文给出的多孔溢洪道闸墩后水流汇合完成长度的经验公式和水流单侧横向扩散宽度的经验关系式，可应用于工程实际，适用于多孔溢洪道闸墩后溢流坝面上的水流汇合和横向扩散的水力计算。