燕尾坎水舌冲击特性试验及数值仿真研究

马 斌， 叶星宇， 缑文娟， 王孝群，2，3， 王宇航

(1.天津大学 水利工程仿真与安全国家重点实验室， 天津 300350； 2.河北工程大学 水利水电学院，河北 邯郸 056038； 3.河北省智慧水利重点实验室， 河北 邯郸 056038；)

1 研究背景

挑流消能具有施工方便，造价及维护成本低的优点，在我国高坝工程建设中被广泛使用，然而下泄水体能量集中而带来的巨大冲击力与强烈的水流脉动作用也是此类消能方式的特点[1]，多数工程通过挑坎的体型变化分散水舌、引导水舌、控制水舌的挑射距离等来优化防冲消能效果。近年来双曲面窄缝坎[2]、差动坎[3]、贴角扭面坎[4-5]、舌瓣坎[6]、翻转坎[7]等已经在国内重大工程上得以运用，有的则作为试验成果中的推荐体型用以解决工程上的局部动水荷载过大、岸坡冲刷和泄洪雾化等突出问题。对于高水头大单宽流量的高坝工程，水舌的分散有助于解决下游荷载集中与消能水体利用不充分的问题。根据Deng等[8]的研究，燕尾坎对于水舌的纵向拉伸效果明显，其团队与国内学者对这一新型挑坎的水力特性进行了相关的试验和研究[9-12]，在国内外一些重大工程的消能防冲设计中燕尾坎已经取得初步的实践效果[13-14]。总体来讲燕尾坎的研究和运用历史尚短，结合大型重力坝坝身泄洪应用的研究则更少，因此对其水力特性的研究仍具有学术价值。经过长期的探索和发展数值模拟技术已经成为泄洪消能研究的重要手段，燕尾坎出挑水舌流态复杂，可以借助数值模拟的可视化优势来详细了解其水力特性。从Xu等[15]的研究结果来看，雷诺时均法的k-ε模型对高坝泄洪中典型淹没射流的模拟是适用的；郭文思等[16]将k-ε模型模拟的射流场结果与PIV(particle image velocimetry)测量结果进行了比对，进一步验证了该湍流模型模拟此类淹没射流的可行性。Yakhot等[17]将重整化群RNG (renormalization group)理论引入湍流领域，在Standardk-ε双方程模型的基础上建立了RNGk-ε湍流模型，此后该模型便广泛应用于求解大弯曲流动问题。从文献 [18]～[22]的研究成果来看，RNGk-ε湍流模型对于扭面坎、窄缝坎、斜切坎、单孔双窄缝坎和多股水舌等复杂水舌的泄洪过程均能取得理想的模拟效果。

以上研究成果多以求解的水舌参数、流场、时均压强等作为检验数值模拟精度的指标或作为求解目标，但运用该湍流模型求解水垫塘边壁脉动压强的研究很少。雷诺时均法在时间域上对湍流场中的物理量进行了平均，从而失去了部分脉动信息。结合梁在潮等[23]、崔广涛等[24]、刘昉[25]的研究来看，大涡旋的紊动作用是水垫塘冲击区边壁脉动的主要成因，压力脉动呈现出幅值大、频率低的特点，低频脉动的能量贡献具有突出优势。雷诺时均法将每一步迭代计算的结果进行时均化处理，但在模拟时间序列中每一步计算结果时仍然具有瞬态性，若迭代计算的时间频率远大于模拟流场中水流的脉动频率，那么来描述该频率范围内低频的水流脉动理论上是可行的。

本文对连续坎和燕尾坎的泄洪过程进行了模型试验和数值模拟，对比了两种挑坎的水力特性和水垫塘底板动水荷载的分布特性；同时通过对两种挑坎泄洪过程的数值模拟验证RNGk-ε湍流模型求解水垫塘冲击区低频脉动的可行性。

2 研究方法

2.1 物理模型试验

模型试验依托黄河古贤水利工程，在天津大学高速水流实验室进行，模型按重力相似准则设计，几何比尺为1∶60，泄水建筑物与防冲建筑物制作材料为亚克力板，糙率为0.007，如图1所示。

图1 试验模型及燕尾坎消能工体型装置(单位：m)

泄水建筑物拟采用挑流设计结合重力坝身布置，表孔进口高程为612 m，挑流鼻坎末端高程为480 m。两种挑坎挑角及反弧半径相同。下游设消能防冲水垫塘，水垫塘底板高程为450 m，水垫塘长度为300 m，末端二道坝净高为22 m；以水垫塘底板为基准，校核水位为178.75 m，本文模型试验与数值模拟均在该水位下进行。

泄水建筑物与下游防冲水垫塘沿表孔中轴线和水垫塘中轴线对称，物理模型与数值模型建立相同的坐标系，坐标原点设置在挑坎末端，即水垫塘起始点x=0，水垫塘末端x=300 m，数值模型和物理模型关于xoz面对称。

时均压强采用测压排测量，脉动压强采用YPS300微型数字压力传感器和SDA2000数据采集分析系统测量。根据奈奎斯特采样定理[26]，采样频率必须大于脉动频率的2倍。根据前人研究成果[27-29]可知，水垫塘中的水流脉动频率主要集中于5 Hz以内，传感器采样频率和数值模拟压力采样信号频率均取100 Hz，采样时长为180 s。

动水压强以1 m高度水柱的压强为计量单位，即9.8 kPa=1 m。时均压强测点沿水垫塘底板中轴线布置，从67.0～297.5 m处依次布置32个测点，水舌冲击区测点间距为6 m，冲击区上、下游间距为12 m。脉动压力传感器根据不同挑坎体型其位置做出相应调整。数值模型中监测点布置于水垫塘底板上表面0.01 m处，中轴线测点与试验测点保持一致；在水舌冲击区y向布置8列脉动压强测线，每列21个测点，共计168个，如图2所示。

图2 数值模型脉动压强测点布置

2.2 数值模型仿真

2.2.1 数学方程 RNGk-ε湍流模型紊流控制系列方程如下。

连续方程：

(1)

动量方程：

(2)

k方程：

(3)

ε方程：

(4)

(5)

(6)

式中：ρ为体积分数平均密度，kg/m3；μ为分子黏性系数；g为重力加速度，m/s2；p为修正压强，Pa；k为湍动能，J；ε为湍动能耗散率；μt为紊流黏性系数；Gk为平均流速梯度变化引起的紊动能项；σk和σε分别为k方程和ε方程的紊流普朗特数；C1ε和C2ε为ε方程常数，ε方程常数值如表1所示。

表1 RNG k-ε控制方程各常数值

2.2.2 计算域及网格划分处理 模拟范围涵盖上游库区，其范围为x向100 m、y向80 m，水库底部距进口高程100 m，中表孔泄槽、下游水垫塘、下游20 m河床及自由液面以上为水气流动区域。整体网格均采用六面体结构化划分，网格正交性与连续性良好。网格全局尺寸为1.0 m；下游水垫塘水流行进方向网格尺寸设置为1.0 m，垂直于水流向网格尺寸为0.8 m；流道弯曲段和水舌冲击区域以0.5～0.8 m进行局部加密；水舌运动区域水体与空气发生强烈的卷气和能量交换，其模拟精度关系到水舌的冲击效果，该部分网格与周边空气运动区独立划分，尺寸为0.5 m；水舌运动区与周边空气运动区采用interface连接，整体网格数量为236×104个，流体域网格划分如图3所示。

图3 数值仿真流体域网格划分

2.2.3 求解器及边界条件设置 数值模拟采用压力基求解器，计算时间为瞬态，压力速度耦合求解器选择PISO (pressure-implicit with splitting of operators)，动量方程离散取二阶迎风格式，湍动能及耗散率均为一阶迎风格式，收敛条件保持默认设置。

上游库区进口设置为压力进口给出等水位的压力水头，流速水头设置为0；下游出口设置为压力出口；自由液面、水舌运动区域和燕尾坎镂空段设置为压力进口，设置水的体积分数为0；壁面均设置为无滑移固壁边界，粗糙度设置为0.5 m，粗糙高度设置为0.001 m。整个模拟过程在标准大气压101.32 kPa环境下进行，空气密度设置为1.225 kg/m3。

3 水舌形态及挑坎水流运动特性

3.1 水舌形态

取数值模拟结果与试验结果对比两种挑坎的水舌形态及射流流态，进而研究其不同的水舌运动特性和冲击特性。两种挑坎的水舌形态及射流流态如图4和5所示，水舌形态参数试验值与模拟值汇总见表2。

图4 两种挑坎的水舌形态对比

图5 数值模拟两种挑坎射流流态纵剖面图(t=60s)

表2 两种挑坎的水舌形态参数汇总 m

由图4和5可看出，连续坎水舌出挑后有部分横向扩散，经过最高点后舌心区周边水体开始分散，至末端部分水体已经破碎。水垫塘上游0～170 m段水面无明显波动，强烈的紊动发生在冲击区及下游。燕尾坎泄洪过程中下泄水体行进至燕尾开口处时，水体下边界失去约束水体先从中间缺口处起挑，另一部分水体从两侧底板起挑。相比于连续坎，其水舌后缘大幅度向上游移动，前缘挑距也有明显的减小，水舌纵向获得大幅度拉伸，侧面呈扇形。水舌上缘与下缘之间水体掺气效果明显，从水垫塘上游x=70 m处水体开始受水舌后缘水体冲击，冲击区范围大幅度增加，水垫塘水体的紊动程度减弱且紊动区范围扩大。

两种挑坎水舌横断面的沿程变化如图6所示。由图6可见，连续坎水舌出挑后向两侧扩散，中间水体沿z向上凸，沿x向运动趋势增强，水体主要分布于断面中上部；燕尾坎水舌出挑后断面则呈下凹状，沿x向运动，两侧水体向射流中轴线集中，至x=80 m处断面开始向z向上凸，相比于连续坎水舌其横向扩散幅度减小，水体集中于射流上部，射流下部掺气明显。

图6 两种挑坎水舌横断面沿程变化

3.2 挑坎水流运动特性

为了进一步探究两种挑坎的水舌扩散特性，对两种挑坎的水流运动特性进行分析。如图7所示，截取t=60 s时刻挑坎x=-17 m、x=-11 m、x=-5 m处3个横断面，提取各横断面上的水气二相分布、压强和不同方向流速分量u、v、w的时均分布云图，u、v、w分别为x、y、z方向的流速分量。

图7 两种挑坎水流运动特性研究断面选取

3.2.1 挑坎断面水气二相及压强分布 图8、9分别为两种挑坎各选取断面的水气二相分布和压强分布情况。

图8 两种挑坎各断面水气二相分布

图9 两种挑坎各断面压强分布

如图8、9所示，连续坎断面上水体的自由液面靠近挡墙的两侧，略高于流道中间部分，整体呈水平分布，沿x方向没有明显的起伏变化；燕尾坎上水体行进至开口处时液面已经有下凹迹象，沿x方向过水断面呈“漏斗状”向底板缺口处拉伸。连续坎上压强在z方向近似线性分布；燕尾坎压强两侧大中间小，在y方向形成明显的梯度变化；两种挑坎上的断面压强分布均沿x向减小。

3.2.2 挑坎断面流速分布 两种挑坎各选取断面流速分布如图10所示。分析图10可知，连续坎上水流流速分布较为均匀，除固壁边界的粘滞作用影响整体分布均匀，没有明显的梯度变化，流速u沿x向运动过程中受挑坎阻力影响流速减小；流速v沿x向运动过程中没有明显的梯度变化；由于挑坎上挑方向与流速w保持一致，因此流速w沿x向则有明显的增大。燕尾坎流速u的核心区集中于中间缺口处，沿x向运动过程中向下移动；流速v在沿x向运动过程中两侧流速增势明显，在y向与中间部分形成较大流速差。

图10 两种挑坎各断面x、y、z向的流速分量分布

如图11所示，进一步对两种挑坎末端x=-5 m断面上的流速u、v、w分量按监测面Ⅰ、Ⅱ分块进行局部监测，连续坎和燕尾坎监测面Ⅰ、Ⅱ流速u、v、w均值汇总分别如表3、4所示，其中，u+w为流速u和流速w在平行于平面xoz的纵断面上的合速度；v+w为流速v和流速w在平行于平面yoz的横断面上的合速度；θ11和θ21分别为监测面Ⅰ、Ⅱ上平行于xoz的纵断面上合速度与x轴正方向的夹角；θ12和θ22分别为监测面Ⅰ、Ⅱ上平行于yoz的横断面上合速度与z轴正方向的夹角。

图11 两种挑坎x=-5 m断面流速监测面分块

表3 两种挑坎监测面Ⅰ流速均值汇总表

表4 两种挑坎监测面Ⅱ流速均值汇总表

由表3、4可知，连续坎监测面Ⅰ、Ⅱ合速度u+w分别为47.64、50.39 m/s，挑角分别为24.71°、26.00°，根据类平抛运动，监测面Ⅱ对应的中间部分水体起挑速度和起挑角度更大，成为水舌外缘；监测面Ⅰ对应的起挑水体成为水舌内缘，但两者相差不大，这与连续坎水舌入水宽度集中的特点相符合。燕尾坎监测面Ⅰ、Ⅱ合速度u+w分别为41.52、49.46 m/s，挑角分别为21.81°和9.90°，可见燕尾坎中部起挑水舌与两侧底板起挑水舌存在较大的挑角差，两侧水体获得更大的挑角成为水舌外缘，中间水体成为水舌内缘；两者合速度的巨大差异形成燕尾坎水舌纵向的大尺度拉伸。

对比监测断面Ⅰ连续坎和燕尾坎合速度u+w，两种挑坎分别为47.64、41.52 m/s，挑角分别为24.71°、21.81°，由此可以解释等反弧半径挑角的燕尾坎前缘挑距小于连续坎。

对比两挑坎监测面Ⅰ、Ⅱ速度v均值的变化，两者断面中部的速度v近乎为零，两侧水体的流速燕尾坎和连续坎分别为2.07和0.18 m/s，可知燕尾坎y向流速梯度远大于连续坎，在监测面Ⅰ中v+w与z轴正方向的夹角θ12为7.65°，这也解释了燕尾坎水舌在横向上的收缩现象。

4 水垫塘底板压强分布特性

4.1 时均压强分布特性

两种挑坎相应的水垫塘底板时均压强试验值与模拟值分布如图12所示。由图12可看出，连续坎泄洪时x=50～150 m区间时均压强分布较为平缓；x=210～250 m为其冲击区，该区间主流受水垫与底板的折冲作用流速急剧减小，动能转化为势能形成冲击峰，峰值点在x=234.5 m；壁射流区水体附着于底板向上下游扩散，该部分水体动能增加、势能减小、时均压强在峰值点上下游呈下凹状。燕尾坎时均压强分布整体呈上升趋势，x=90～200 m为水舌冲击区，时均压强在该区间形成两个连续的冲击峰，其局部峰值点在底板x=145.0 m与x=186.5 m，相比于连续坎时均荷载核心区向上游移动了70 m，冲击区范围大幅增加，形成的冲击压强减小。

连续坎峰值点试验值与模拟值分别在x=234.5 m和x=240.5 m，其冲击压强值分别为3.09和3.67 m；燕尾坎时均压强试验值与模拟值的分布规律相似，两者的数值误差均在5%以内。两种挑坎时均压强的沿程分布规律与试验结果吻合，其数值误差均在10%以内。

4.2 脉动压强均方根分布特性

泄洪过程中水垫塘的水流运动属于紊流，因此流场中压强的变化随时间呈现瞬态性特征。图13为连续坎泄洪过程中水垫塘底板冲击区时均压强峰值点去均值后的试验实测和模拟压强时程曲线，由图13可看出，模拟压强信号与实测压强信号分布规律相似。

图12 两种挑坎相应的水垫塘底板时均压强试验值与模拟值分布

图13 连续坎水垫塘底板冲击区峰值点去均值后压强时程曲线

取连续坎冲击区时均压强峰值点压强时程信号进行傅里叶变换，得到的水垫塘冲击区典型功率谱如图14所示。分析图14可知，试验值与模拟值的主频分别为0.016 7和0.013 8 Hz，优势频率均集中于0～1 Hz；从频谱的能量分布来看试验值与模拟值的功率谱能量主要集中于0～2 Hz，由于真实的流场包含更小的涡旋结构， 因此试验结果在更高频率段也有一定的能量占比，模拟功率谱大于3 Hz的能量分布几乎为0，这符合RNGk-ε湍流模型本身的低通滤波属性；从整体分布趋势看在0～2Hz能量均随频率的增加而递减，这与前人对水垫塘冲击区频谱特性的研究结论是吻合的。

两种挑坎相应的水垫塘底板脉动压强均方根试验值和模拟值的沿程分布如图15所示。由图15可看出，连续坎在x=50～180 m区间脉动压强均方根均在1 m以内；x=180～280 m区间为波峰段，峰值点在x=234.5 m，峰值为5.92 m；峰值段靠近水垫塘下游，以峰值点为对称点40 m范围内下游均方根值大于上游。相比于连续坎，燕尾坎脉动压强均方根从起始测点至末端测点整体分布呈波峰状，在测线上的分布更为均匀，整体降幅明显，其峰值点在x=207.5 m，相比连续坎向上游移动了27 m，峰值大小为2.18 m，同比降幅高达57.14%。

图14 水垫塘底板脉动压力功率谱密度

图15 两种挑坎相应的水垫塘底板脉动压强均方根值沿程分布

连续坎试验值与模拟值的脉动均方根峰值点分别在x=234.5 m和x=240.5 m，数值大小分别为5.92 m和6.84m；燕尾坎试验值与模拟值峰值点分别在x=207.5 m和x=201.5 m，数值大小分别为2.18和2.41 m。对比模拟结果与实测结果，冲击区脉动压强均方根偏差均在20%以内；相比模型试验中的脉动压强均方根峰值降幅为63.18%，数值模拟中相应的降幅为64.77%，可见数值模拟能较好地反映脉动压强均方根峰值随挑坎体型的变化。

4.3 水舌冲击特性

提取连续坎和燕尾坎两种挑坎水垫塘冲击区监测点的脉动压强均方根模拟值绘制其分布云图，如图16所示。由图16可看出，连续坎脉动压强均方根核心区峰值优势明显(图16(a))，从核心区沿x轴正方向脉动压强均方根呈衰减趋势，而在y轴方向扩散半径沿程增大；沿x轴负方向脉动压强均方根的传递减弱，在y轴方向扩散半径也沿程减小。相比于连续坎，燕尾坎脉动压强均方根整体大幅度减小(图16(b))，沿核心区向x轴正方向及x轴负方向衰减速度也更均匀。两者在y轴方向的衰减速度均大于x轴方向，这也意味着能量的消杀主要集中在射流轴线附近。

图16 两种挑坎相应的水垫塘冲击区脉动压强均方根值平面分布

5 结 论

(1)通过数值模拟结果与试验结果的对比和分析可知，RNGk-ε湍流模型对连续坎和燕尾坎的泄洪过程模拟效果良好，水舌形态、水垫塘底板时均压强分布的模拟结果与实测结果基本吻合，脉动均方根值的对比结果表明该湍流模型求解水垫塘大尺度水流脉动是可行的。

(2)燕尾坎上水流边界的突变使断面压力和流速分布产生了明显的横向梯度，梯度指向断面中轴线；挑坎中间缺口出挑水流流速大、挑角小，形成了水舌内缘，两侧出挑水体流速小、挑角大，沿程向射流轴线集中形成了水舌外缘，水舌横向收缩明显而纵向得到大幅度拉伸。

(3)与连续坎相比，燕尾坎泄洪时水舌冲击区范围及水垫塘水体紊动区范围增大，水舌冲击区紊动程度减弱；水垫塘底板的时均压强和脉动压强核心区分别向上游移动了70和27 m，冲击压强和脉动压强均方根峰值均大幅减小，荷载能量集中程度减弱。