高水头单侧渐扩消力池水力特性数值模拟研究

赵东阳，尹进步，王国杰，杨 彦，杨 钊，张曙光

(西北农林科技大学水利与建筑工程学院，陕西 杨凌712100)

0 引 言

高水头、大流量的高坝泄洪消能问题，往往是大型水利水电枢纽工程设计中重大而复杂的科学技术问题。需要综合考虑挡水建筑物、泄水建筑物和兴利建筑物三者之间的关系和合理布局，并根据坝址的地形、地质和水文条件，合理选取泄洪方式、泄水建筑物的型式和消能方式[1]。同样，本次研究借助的高水头岸边溢洪道之所以采用底流，是由于鼻坎位于岸边，距离主河床较远，且水舌挑距有限，无法落到河床，故只能放弃挑流采用底流。原方案消力池边墙不扩散，通过模型试验发现由于跃后水深不足，消力池内未能形成淹没水跃，为增加消能水体，采用边墙扩散来增加消力池尺寸。由于该消力池右侧边墙为山体且坡度较陡，开挖工程量和边坡支护量较大，而左侧边墙在河道内，开挖量不大，因此采用单侧渐扩。

目前使用数值模拟分析消力池内水力特性的研究较多，Iwao，Katakam研究发现对称扩散消力池，消能率较高，流态较好[2,3]。王月华利用数值模拟提供的消力池Fr值、紊动能分布图和流速矢量图，比较全面地反映了三级消能池水流和消能情况[4]。郑梅玲通过数值模拟分析了消力池内部流场，发现渐扩式能使水流横向扩散的空间增大，增强消力池首部水流的立轴旋滚[5]。上述只是关于对称扩散消力池的研究，而针对单侧渐扩消力池的研究较少，本文参考相关工程资料[6-11]，采用Flow-3d软件中RNGk-ε双方程紊流模型基于VOF方法，对某工程单侧渐扩消力池进行数值模拟，得到水面线、流速、压强和紊动能等水力参数，并与模型试验结果进行对比分析，为工程设计提供参考。

1 数学方法及模型建立

1.1 控制方程

采用RNGk-ε双方程湍流模型对消力池进行数值模拟，连续性方程，动量方程和k-ε方程表示如下：

连续性方程：

(1)

动量方程：

(2)

紊动能k方程：

耗散率ε方程：

(4)

1.2 VOF方法

采用VOF方法[12]追踪水流自由表面，该方法定义流体体积函数F=F(x,y,z,t)表示计算区域内流体的体积占据计算区域的相对比例。对于某个单元，F=1表示该单元被流体完全充满，F=0表示该单元是个空单元，没有流体，F=0～1表示该单元被流体部分充满，VOF中自由表面的复杂变化可用函数描述为：

(5)

式中：Ax、Ay、Az分别为x、y、z三个方向可流动的面积分数；u、v、w为流速；VF为可流动的体积分数。

1.3 模型建立

计算区域包括进水渠、控制段、泄槽和消力池及部分海漫，总长323.00 m(见图1)。溢洪道为岸边开敞式有闸控制正槽溢洪道，控制段设置两孔闸门，泄槽陡坡段设有两道掺气坎。溢洪道出口采用底流消能，消力池采用单侧渐扩式，底板高程为834.00 m，进口宽14.00 m，出口宽20.37 m，长84.53 m，消力池出口采用直立尾坎，为增加跃后水深，在直立尾坎上设置高2.50 m、宽1.50 m不连续的梯形墩，边墙扩散角度为5.6°。查阅资料得知，实际工程扩散角度一般较小，通常小于5°[13]，与本工程差别不大。此处校核工况下消力池底板以上水头为89.08 m，水头相对较高。工况情况为：校核工况闸门全开，上游水位为923.08 m，下游水位为846.81 m，泄量为513.30 m3/s，消能防冲工况时闸门对称开启，开度为3.68 m，上游水位为921.50 m，下游水位为846.08 m，泄量为280.00 m3/s。定义溢流堰上游立面前32.50 m断面为X=0+00.00 m、溢洪道右侧边墙外立面为Y=0 m(右侧墙内壁为Y=2.50 m位置)、消力池底板下部6.10 m高程的水平面为Z=0 m，即消力池底板高度为Z=6.10 m。计算区域全部用结构化网格来划分，单元尺寸为 ，网格总数大约790万。

图1 溢洪道结构(单位：m)Fig.1 Layout of the spillway

边界条件：计算区域上下游为压力边界，并给定水位；底部及左右两侧均为壁面边界；顶部为压力边界。压力求解器选用广义极小残差算法(GMRES)，基于压力隐式求解法(Implicit)；计算结果数据输出间隔2 s，输出数据有流速、压力、水的体积分数、自由液面高程等；初始时间步长定为0.000 1 s。

2 计算与试验对比

为确保计算结果的可靠性，此处分别将消能防冲工况和校核工况下计算得到的水面线、流速、压强等水力参数与相应工况的模型试验数据进行对比。用于对比的试验模型按重力相似准则设计，模型比尺为 1∶40。

2.1 水面线

两工况下计算水面线与试验水面线对比见图2。从图2(a)可以看出，消能防冲工况试验水面线与计算水面线吻合较好，在桩号0+224.26 m两者重合，相对误差为0%。在桩号0+280.00 m处，两者相对误差最大，为12.30%，但考虑消力池内水流翻滚剧烈，水面波动较大，仍认为具有一定的可信度。

图2 水面线对比Fig 2.Comparison of surface line

校核工况试验水面线与计算水面线对比由图2(b)可见，两者相差甚微，最大相对误差出现在桩号0+290.00 m处，仅为7%。

对比两个工况计算结果发现，随着流量的加大，水面爬升梯度越大。由于此消力池边墙扩散角度较小，仅5.6°，计算扩散侧水面线(见图2)与计算中心线位置水面线对比发现两者基本重叠，故而无必要对渐扩侧水面线进行研究。

2.2 流 速

此处首先对计算得出的中心线处底、表流速与相应的试验值进行对比如图3所示，可见两者相差无几，在此基础上，对计算中心线位置沿程流速分析之后，又详细分析了水平面流场。

图3 消力池流速对比Fig.3 Comparison of stilling basin velocity

消能防冲工况下沿程流速分布见图3(a)所示，可见在桩号0+253.00 m以前，底部正向流速较大，超过10 m/s，在桩号0+268.50 m以前，表面流速为负值，说明该区域表面水流流向为回流。在桩号0+259.00 m以前，表面反向流速较大，超过3 m/s。

分析校核工况下沿程流速分布[见图3(b)]发现，表面流速在桩号0+279.50 m以前为负值，水流流向为回流。反向流速在桩号0+270.00 m以前流速较大，超过3 m/s。底部正向流速在桩号0+270.00 m以前流速较大，超过10 m/s，对两个工况计算与测试结果对比发现，在桩号0+224.50～270.00 m区域，两工况底部正向流速超过10 m/s，表面反向流速超过3 m/s，对实际工程而言，该正向流速与反向流速可能偏大，容易形成底板磨蚀破坏。

试验发现扩散侧表面时而存在立轴旋涡，消能防冲工况时在桩号0+278.00 m位置[见图4(a)]，校核工况时在桩号0+289.00 m[见图4(b)]。然而消力池内流态复杂，水体来回波动频繁，试验测量池内水面及以下流场不便，为了对池内流场分布规律进行深入探索，下面采用数值模拟对水平面流场进行分析。

2.3 水平面流场

平面流场能清楚反映水流的运动轨迹，为准确捕捉流场的运动过程，每个工况分别均截取三个高程的水平面流场，如图5所示，消能防冲工况截取的平面分别为Z=6.20 m，Z=10.20 m和Z=18.73 m，校核工况截取的平面分别为Z=6.20 m，Z=8.24 m和Z=18.73 m。两种工况Z=18.73 m截面均无上游流速矢量，是由于消力池上游水面较低，未被剖切到。

图4 消力池流态图Fig.4 Flow pattern of still basin

图5 消力池内不同高程流场分布图Fig.5 Distribution of flow in the still basin

消能防冲工况下水平面流场见图5(a)，Z=6.20 m截面接近消力池底板，渐扩侧尾部的立轴旋涡涡心在桩号0+280.00 m(与试验测得桩号0+278 m接近)，范围大约在X=281.4～290.3 m，Y=19.26～20 m，漩涡范围较小。由Z=10.20 m截面可见左边渐扩侧在桩号0+221.75 m处产生翻滚漩涡，即跃首位置，同时渐扩侧尾部的立轴旋涡基本消失。由Z=18.73 m水平面流场可见，旋滚区发展到0+275.00 m停止，即跃尾位置。

校核工况见图5(b)，由Z=6.20 m流场分布图可见，在消力池渐扩侧尾部的立轴旋涡涡心桩号为0+291.00 m(与实测的桩号0+289.00 m接近)。旋涡范围较大，大约为X=281～298 m，Y=13.52～20 m。观察Z=8.24 m平面流场可见，跃首在桩号0+229.01 m处，位于左边渐扩侧的漩涡范围缩小。由Z=18.73 m得出，跃尾在桩号0+296.00 m位置。跃尾横轴在渐扩侧向下游倾斜的角度较消能防冲工况更大。

对比两个工况水平面流场发现，渐扩侧尾部存在立轴旋涡，靠近消力池底板旋涡范围最大，随高程升高逐渐减小，且随入池流量加大，旋涡范围越大。跃首横轴在渐扩侧向上游移动，跃尾横轴在渐扩侧向下游移动。分析其原因，水流受到右侧边墙约束，只能向左边约束较小的渐扩侧扩散，由此在底板附近产生的立轴旋涡将跃尾横轴推向下游，将跃首横轴推向上游。

2.4 时均压强

试验中消力池中心线处时均压强采用玻璃管测试，由试验结果与计算结果对比图图6可见，两工况在桩号0+260.00 m以前计算与试验结果相差甚微，最大相对误差不足5%，在桩号0+260.00 m以后误差稍大，消能防冲工况相对误差小于10%，校核工况相对误差最大值仅为5%，计算与试验结果基本吻合。

各工况下时均压强基本均呈现先增后减再增之后逐步趋于稳定，消力池内均未出现负压。桩号0+221.50～232.30 m区域为反弧段，受离心力作用，水流剧烈碰撞消力池底板，造成局部压强迅速增大。之后离心力减小，压强逐渐下降。在桩号0+242.50 m位置，校核工况比消能防冲工况的上游水头大，弗汝德数大，跃前水深小，即静水压强小，所以该位置校核工况时均压强值比消能防冲工况时小。该部位流速较大，压强较低，在实际工程中，掺气量不足的情况下，较容易破坏，而本工程中设有两道掺气坎，且水面较薄，掺气量足够，因而发生破坏的可能性较小。桩号0+242.50～270.00 m区域上部为水跃，水面迅速爬升，且水跃有向下游的推力，因此压强回升较快。桩号0+270.00 m之后水面逐渐平稳，压强也随之平稳。

扩散侧压强如图6所示，计算中心线位置与计算扩散侧水面线基本重叠，在消能防冲工况下，两者相对误差最大仅为8‰，校核工况下最大仅差2%，差别甚微，因此文中不再对扩散侧压强进行分析。

图6 消力池压强对比Fig.6 Comparison pressure stilling basin

2.5 探讨计算紊动能与试验脉动压强

紊动能是能够体现水流由紊动而产生的能量耗散程度的重要指标，其可以衡量流体的紊动状态，对考察水跃消能的防冲建筑物性能具有重要意义。本文采用的 RNGk-ε双方程紊流模型可以计算水流的紊动能：

(6)

式中：u′为脉动流速。

图7 消力池底板紊动能云图Fig.7 Calculated turbulent energy of stilling basin

由消力池底板紊动能云图图7可见，紊动能向渐扩侧扩展明显，随水位流量加大，范围逐渐向消力池中心线移动。因此消力池重点防护区域偏向渐扩侧。

图8 紊动能纵剖面图 Fig.8 Turbulent energy longitudinal section

图9 消力池底板脉动压力Fig.9 Stilling basin floor pulsation pressure

取紊动能最靠近上游位置纵剖面(见图8)，与试验中相应位置脉动压力(见图9)进行对比发现，消能防冲工况下，在桩号0+225.00 m处脉动压力最大，脉动压力均方根最大值为34.80 kPa，对应到紊动能云图上为红色区域的起始位置[见图8(a)]。校核工况下，脉动压力在桩号0+238.00～242.50 m范围最大，脉动压力均方根最大值为47.20 kPa，对应到紊动能云图同样为红色区域的起始位置[见图8(b)]。由此可见紊动能云图可能反映脉动压力分布情况，分析其原因，紊流由大小不同尺度的涡体(或漩涡)所组成，尺寸各异的漩涡在高速水流中旋转，转向与水流方向一致的一边流速增大，因而压强变小，另一边则流速减小压强增大；两边的压强差产生的推力，对于涡体而言就是升力。当升力足以克服水流的黏滞阻力时，涡体将上升做横向运动，进入流速不同的流层，即产生水流质点的混掺作用。由于紊流中有无数个转向与大小不同的涡体，所以水流的混掺作用一直杂乱无章地进行着，从而导致了水流流速的脉动；根据能量方程就导致了压强的脉动，两者直接关联。所以脉动压力必将受到脉动流速影响，紊动能 在一定程度上反映脉动压力分布情况。

3 结 论

本文采用Flow-3d软件对单侧渐扩消力池进行数值模拟，得到了水面线、流速、压强、湍动能等水力参数分布规律，并与水工模型试验结果进行了对比验证，得出结论如下：

(1)Flow-3d软件中的RNGk-ε双方程紊流模型能较好地模拟单侧渐扩消力池中的水力特性，数值模拟结果与模型试验结果吻合较好，为进一步研究奠定了基础。

(2)详细分析消力池内流场，发现跃首横轴在渐扩侧向上游移动，跃尾横轴在渐扩侧向下游移动，且随流量的加大，移动的距离会逐渐加大。

(3)通过水平面流场分析发现，消力池渐扩侧尾部产生立轴旋涡，接近消力池底板漩涡范围最大，向上逐渐减小，且随流量加大，范围将越来越大。

(4)对比计算紊动能和试验脉动压力得出：紊动能在一定程度上反映脉动压力分布情况；反弧段及消力池底板前部的脉动压力较大，设计上应对这些区域加强防护。