带闸墩溢流坝控泄工况二维与三维数值模拟

2011-09-06姜伯乐

水利水电科技进展 2011年6期

李静,李丹,姜伯乐

(1.长江水利委员会长江科学院,湖北武汉 430010;2.武汉大学水资源与水电工程国家重点实验室,湖北武汉 430072)

目前对带闸墩溢流坝过坝水流的三维数值模拟比较少见,这主要是由于三维模拟的计算网格单元数急剧增加,闸墩的存在加大了网格划分的难度并且使整个流场变得异常复杂。因此对过坝水流的研究工作通常较少考虑闸墩对流场的影响,而直接以二维流场为研究对象[1-3]。文献[4]通过分别计算过坝水流和闸墩绕流2个二维模型来近似模拟三维带闸墩的溢流坝流场,为解决三维过坝水流问题提供了一条途径,但是将二维数值模型应用于实际带闸墩溢流坝过坝水流计算还有一定的差距。主要表现在闸墩的存在对水面线位置、坝面压力、溢流流量等均产生显著的影响。文献[5]利用有限体积法对溢流坝三维湍流场进行数值模拟,但没有考虑闸墩的影响。文献[6-14]对溢流坝及湍流场进行了数值模拟研究。笔者以银盘水电站带闸墩溢流坝为例,以雷诺时均方程为基础,利用k-ε紊流模型模拟湍流流动,采用流体体积(VOF)法确定自由水面线的位置,采用有限体积法进行离散求解,并与物理模型试验结果进行对比,比较准确地给出了过坝水流的水面线位置及坝面压力分布,可为泄水工程的水力设计提供比较可靠的设计依据。

1 数学模型

1.1 模型控制方程

VOF模型适用于跟踪捕捉两相或多相流中互不穿透流体间的界面变化。VOF模型通过引入体积分数、求解体积分数方程来捕捉界面[15]。

设某一控制单元内第q相体积分数为αq(0≤αq≤1),则当 αq=0时,该控制单元内无第q相流体;αq=1时,该控制单元内充满第q相流体;0＜αq＜1时,该控制单元内存在相间界面。当控制单元中存在n相流体时,各控制单元内各相体积分数之和为1,即

本文所研究问题属于不可压缩流体的湍流流动问题,因此模型中采用标准k-ε模型,其基本控制方程如下:

1.2 计算区域和网格划分

溢流坝堰顶高程为195m,进口坝面头部采用方程为x2/72+y2/4.0752=1的1/4椭圆曲线,后接WES曲线,曲线方程为y=0.038x1.85,其后由1∶1.8斜坡与半径为69.80m的反弧段衔接。消力池底板高程为175m。闸孔宽度为15.5m,闸墩宽度为4.5m。本文研究工况:上游水位为215m,下游水位为199.66m,进口水深为40m,闸门开度为0.4,闸门开启孔数为8个,下泄流量为10850m3/s。物理模型取2孔,按照重力相似准则设计,长度比尺为1∶50。为了与物理模型结果进行比较,二维和三维模型与物理模型尺度相同。根据WES曲线的特性,x方向坐标原点取在堰顶处。二维模型计算区域见图1。

图1 二维模型计算区域

由于闸墩附近的几何形状十分复杂,为了保证计算精度,对包含闸墩的三维区域进行网格划分时对整个区域进行分区,在各分区上采用不同的网格系统。三维模型计算区域和网格划分分别见图2和图3。

图2 三维模型计算区域

图3 三维模型闸门局部网格划分

1.3 边界条件

气体入口边界为压力边界条件,入口边界处压强为大气压强。水的入口边界为速度入口边界条件。二维计算和三维计算中流速uin均为1.844m/s。入口紊动能k和耗散率ε的边界条件分别由式(6)、式(7)得出:

式中:Hin为进口水深。

1.3.1 出口边界

由于出口压力分布已知,因此出口边界采用压力边界条件,出口边界上压力分布为p=ρg(199.66-z)。过堰水流经过消力池消能后流态仍然非常紊乱。水流过水断面上流速分布需要经过一段足够长的距离才能逐渐调整至渐变流流速分布。因此,在数学模型建立以后,通过大量的重复测试,并与物理模型试验对比,确定了合适的出口边界位置。

1.3.2 壁面边界条件

壁面边界采用无滑移边界,即在壁面上u=v=w=0。近壁区的黏性底层采用标准壁函数处理。

1.4 数值求解方法

本文采用基于非结构网格的有限体积法对计算区域进行离散,并采用PISO算法进行求解。

2 物理模型

本文用于对比的试验数据采用物理模型试验获得。物理模型采用重力相似原则,比尺为1∶50。物理试验主要采集了水面线和坝面压力数据。物理模型试验流态见图4。

图4 物理模型试验流态

3 计算结果比较与分析

3.1 水跃区流线

图5 二维与三维模拟消力池内流线分布

物理模型试验结果表明,该控泄工况下消力池内发生了淹没水跃。图5给出了消力池内水跃的二维模拟与三维模拟结果。图5(a)为二维模拟计算结果,图5(b)～(d)分别为三维模拟第 1孔、闸墩和第2孔中轴线剖面上的流线分布。从图5可以看出,三维计算中由于模拟的对称性,第1孔与第2孔的流线分布基本相同,二维计算的流线分布中存在2个涡,而三维计算的流线分布中只存在1个涡。在三维模型中计算区域、网格划分是完全对称的。但从图5(b)、5(d)中可以看出第 1孔与第 2孔中轴线剖面上的流线分布并不完全一致。这一现象看似与计算区域对称性矛盾。实际上,在三维模拟中水流与空气的掺混机理十分复杂,而且两者的掺混过程属于非恒定过程,因此三维模拟无法得出完全对称的结果。

3.2 水面线位置

图6给出了水面线的二维模拟、三维模拟以及与物理模型结果的对比。由图6可以看出,二维模拟得出的水跃发生位置与物理模型结果有较大偏差,三维模拟得出的水面线与物理模型结果吻合较好,水跃发生位置与物理模型结果比较一致。在水跃发生点之前,二维与三维模拟得出的水面线与物理模型结果均十分接近。但是水跃发生位置2种模拟结果显著不同。二维模拟的水跃发生位置比三维模拟提前许多,出现了较大的偏差。

图6 二维模拟、三维模拟与物理模型水面线对比

3.3 坝面压力分布

图7 二维模拟、三维模拟与物理模型坝面压力分布对比

图7给出了坝面压力分布的二维模拟、三维模拟与物理模型结果的对比。坝面压力采集范围为桩号-11.70～37.81m。从图7中可以看出,由于闸墩的影响,进口压力二维模拟与三维模拟结果有较大差别,三维模拟结果与物理模型结果吻合得更好,但在水跃区二维模拟与三维模拟得出的压力相差不大。坝面压力的数学模拟结果与物理模型结果相比还存在一定的差距,主要原因是数学模型是在物理模型的基础上对进口条件进行简化后建立的。进口条件的简化对计算结果造成的影响在二维模拟中表现得非常突出。由于没有闸墩、进口侧收缩的影响,堰顶上游段水流流速比物理模型试验所测得的流速大,因此造成二维模拟结果中进口段坝面压力偏小。

4 结语

以银盘水电站带闸墩溢流坝为例,用紊流模型模拟了控泄工况下二维和三维数值模拟水面线位置、坝面压力分布以及流线分布情况,并与物理模型结果进行了对比。与二维模拟结果相比,三维模拟结果与物理模型试验结果吻合程度较好,特别是对进口压力与水跃发生点的模拟。但由于数学模型进出口边界条件与物理模型存在差异,数学模拟结果与物理模型试验结果存在一定的差异。相比二维模拟中进出口边界条件简化造成的影响,三维模拟受到的影响较少,因此推荐采用三维模型进行模拟。

[1]苏远波,马震岳,张宏战,等.太平哨水电站溢流坝闸墩裂缝成因数值分析[J].水利水电科技进展,2009,29(S1):103-106.

[2]陈群,戴光清,刘浩吾.阶梯式溢流坝研究综述[J].水利水电科技进展,2002,22(1):47-50.

[3]魏文礼,刘玉玲,金忠青.溢流坝面明流流动水力特性数值分析[J].河海大学学报:自然科学版,2000,28(3):107-110.

[4]谭仁忠,王世夏.有闸墩的溢流坝过坝水流数学模型[J].水利学报,1988,19(8):40-44.

[5]马福喜,王金瑞.三维紊流数值研究[J].水动力学研究与进展:A辑,1995,10(2):115-124.

[6]王志东,汪德.含闸墩溢流坝三维过坝水流数值模拟[J].水科学进展,2004,15(6):735-738.

[7]王志东,汪德.不同墩型溢流坝水动力特性研究[J].河海大学学报:自然科学版,2004,32(4):391-394.

[8]陈群,戴光清,刘浩吾.带有曲线自由水面的阶梯溢流坝面流场的数值模拟[J].水利学报,2002,33(9):20-26.

[9]陈永明,吉庆丰.阶梯溢流坝水流数值模拟及消能分析[J].灌溉排水学报,2006,25(2):68-71.

[10]刘晓川,杨纪伟,李书芳,等.阶梯溢流坝自掺气水流的数值模拟[J].水利科技与经济,2008,14(3):190-192.

[11]程香菊,罗麟,赵文谦,等.阶梯溢流坝自由表面掺气特性数值模拟[J].水动力学研究与进展:A辑,2004,19(2):152-157.

[12]FERZIGER J H,PERIC M.Computational methods for fluid dynamics[M].Berlin:Springer,1996.