溢洪道斜切型挑坎挑射水舌三维数值模拟

溢洪道斜切型挑坎挑射水舌三维数值模拟

1在水舌和工程中的应用集水区能源消耗是当前节水工程重要的节水手段之一。它结构简单，投资省，在项目中得到广泛应用。理论分析和模型试验是研究挑流水舌水力特性的传统方法,目前已经有不少相关研究成果,如刘宣烈等给出了空中水舌的运动机理、掺气水舌运动轨迹方程、水平射距公式、入水速度和曲线长度表达式等;刘士和等通过试验得到了水舌断面形态的总体变化特征,并给出了断面面积和湿周的沿程变化。由于水舌运动具有强紊动性,模型测量和原型观测都具有一定的难度,尤其当模型比尺较大时,试验数据精度受测量设备、手段影响较大,如挑射水舌的空中形态、空中水舌流速、水舌入水角度等,常规仪器尚难以精确观测。因此,近年来一些学者引入数值模拟作为水舌理论和试验研究的补充,如刁明军等模拟了二维情况下射流在空中的轨迹线和水垫塘的自由水面及底板冲击动压;沙海飞等建立了垂向二维水气两相流紊流数值模型,计算得到了坝后挑流水舌风场;叶茂等对导流洞挑射水舌进行了模型试验和数值模拟研究。但是水舌运动属于非线性的水气二相流问题,水舌在空中本身掺气量很大,且上下缘都与空气接触,存在自由面,厚度也较薄,对计算区域网格划分要求较高,因而水舌在空中的形态模拟一直是计算流体力学的难点之一,尤其是进行三维数值模拟。本文结合某狭窄河道水电工程溢洪道挑流消能进行数值模拟研究。由于下游河道较窄,溢洪道轴线与下游河道交角较大,水舌的扩散受到限制,枢纽的消能防冲成为该工程关键性控制技术难题。为了避免溢洪道水舌冲刷下游河道左岸,使水舌尽可能沿河道纵向充分扩散且落入下游河道中心,从而有利于消能和保护岸坡稳定,经过试验优化采用了斜切型挑坎。本文利用优化后的挑坎体型建立三维模型,选用RNGk-ε紊流模型,以流体体积方法(volumeoffluid,VOF)跟踪自由水面,较为精确地模拟了空中水舌形态及水舌入水区域的流态,为试验提供了可靠的依据。2数学模型2.1水的体积分数w动量方程:k方程:ε方程:式中:ρ和μ分别为体积分数平均的密度和分子黏性系数;p为修正的压力;μt为紊流黏性系数,它的表达式为:,其中Cμ为经验常数,Cμ=0.0845;常数αk=αε=1.39;,其中,η0=4.377,β=0.012,C1ε=1.42;常数C2ε=1.68;Gk为由平均速度梯度引起的紊动能产生项,它可由下式定义:。自由水面的跟踪采用VOF方法,该方法假设二相或者多相流体之间不发生质量交换,在一个计算单元内,定义了每相流体的体积分数,单元内所有流体相的体积分数之和为1。如果在单元中第q相流体的体积分数记为αq,那么就会出现以下3种可能:αq=0表示第q相流体在单元中是空的;αq=1表示第q相流体在单元中是充满的;0<αq<1表示单元中包含了第q相流体和一相或者其他多相流体的界面。这里ρ和μ由体积分数加权平均得出,用αw表示水的容积分数,用下标w和a分别代表水和气体,则密度ρ和μ可由下式表示:水的体积分数αw的控制微分方程为:式中:t为时间,ui和xi分别为速度分量和坐标分量。水气界面的跟踪即通过求解该连续方程来完成。2.2广义源项s式中:ϕ为广义变量,可以为速度、温度或浓度等一些待求的物理量;Γ为相应于ϕ的广义扩散系数;S为广义源项。对于速度压力耦合采用PISO算法,与SIMPLE算法相比,PISO算法增加了一个修正步,能更好地同时满足动量方程和连续性方程,在每一时间步内的计算精度高于SIMPLE算法,而且减少了收敛所需要的计算时间,对于瞬态问题有明显的优势。3计算区域网格和边界条件的定义3.1计算网格划分为保证数值模拟来流条件与实际情况相吻合,建模时考虑了上游库区、溢流堰、泄槽段、掺气坎,其网格间距均采用1m(见图1)。由于水舌区和下游河道是本次模拟的重点研究区域,故对该区域的网格进行了加密,同时也兼顾计算机的计算能力,水舌区和水舌与河道水面碰撞区域网格间距设为0.5m。计算区域网格以结构化网格为主,挑坎段少量为非结构化网格(如图2所示),总的计算单元约为163万。模型上游库区正常蓄水位高程为260m,为该工程的正常水位工况,下泄流量为1000m3/s,溢流堰堰顶高程为248m,斜切型挑坎左侧坎出口高程为207.3m,右侧坎出口高程为212.5m,下游河道高程为200m。考虑实际情况下水舌入水区有冲坑存在,根据模型试验数据确定冲坑最低高程为189m,为使计算空间有一定富余,模型最低高程设定为187m。3.2边境条件4计算4.1自由水面的确定本工程下游河道较窄,属于狭长形河道,且溢洪道轴线与河道交角较大,因此优化时主要考虑将水舌横向拉开,使水舌沿河道纵向落入河床。通过模型试验可以看出,斜切型挑坎达到了使水舌纵向扩散的效果,且掺气较明显,解决了水舌集中及入水点靠近左岸的问题(图3)。数值模拟结果采用水的体积分数为0.4确定自由水面。由文献的研究成果可知,水和空气之间的粘滞作用,导致水舌表面在过渡段从外观上看是不透明的连续体,且有很多水颗粒逸出并产生大波纹,分裂段卷入大量的空气使水舌分裂为水束甚至大水片。由于水舌表面水体的逸出、破裂及空气的掺入很难用常规的多相流模型进行模拟,VOF法也是在假设水体不破碎,水气间无混杂存在的基础上进行模拟的,所以模拟水舌的自由水面与实际情况略有差异,但水舌扩散运动的总体形态和实际是一致的。通过观察可以看出水舌在空中扩散较充分,水舌呈扇形散开,厚度较薄,尤其是水舌左侧落水点靠近河道中心,数值模拟水舌形态与模型试验水舌形态吻合良好,参见图4。4.2工程模型试验的数值结果分析水舌挑距是工程中较关注的问题。本文采用该工程1∶40比尺的整体水工模型试验数据作为验证依据。由于水舌形状呈扇贝型,可认为水舌右侧为最大挑距,水舌左侧为最小挑距,通过图5可得数值模拟水舌挑距。表1比较了水舌挑距,数值计算挑距与模型实测挑距基本吻合。误差的出现可能由以下原因造成:(1)数模计算采用的是非恒定流模式,在每一时刻,水面总会有一定的波动,水面高程的波动会导致挑距发生变化;(2)对于计算掺气量大、散裂度高的水体,VOF法具有一定的局限性;(3)受计算机计算能力的制约,数值模拟中水舌入水区域的网格划分还不够细密;(4)模型试验中由于下游水体紊动比较剧烈,很难精确确定水舌入水点位置,试验测量数据也存在一定的误差。在工程中,水舌入水角度及入水速度是研究下游冲刷的重要参数,但是这些数据在模型试验当中很难精确获得,尤其是斜切型挑坎导致挑流水舌成为非对称形态,入水角度的最大最小值往往出现在水舌的左右两侧,因此只能借助于理论分析结果同数值模拟结果对比。文献利用抛射体理论计算公式推导出了水舌的入水角度公式与水舌入水速度公式。本工程中左侧挑坎水舌出挑角度α=11.17°,右侧挑坎水舌出挑角度为α=30°,挑坎位置水流平均初始速度v0=28m/s,左侧挑距L=40m,右侧挑距L=72m,左侧砍顶到水舌入水处水面的垂直距离H=3.3m,右侧坎顶到水舌入水处水面的垂直距离H=11.5m。将以上试验数据代入公式得到理论计算结果,并与图6数值模拟结果对比。分析表2结果:由于文献中的公式未考虑空气阻力因素的影响,所以数值模拟的水舌入水速度略小于理论计算值;因水舌厚度较薄,内缘和外缘入水角度相差不大,数值模拟结果采用水舌内缘入水角较为精确,而公式在此处只适合计算外缘入水角,因此导致模拟结果略大于理论计算值。不过计算结果相对误差均小于3%,说明数值模拟结果与理论计算结果吻合较好,从而验证了数值模拟计算的可靠性。4.3水舌水舌反应时水舌的紊动能分布下泄水流的能量通过水舌空中扩散和在下游河道水垫内剧烈紊动而消耗掉,然而工程中所关心的这些水力参数在模型试验和原型观测中都很难得到确切的数据,但是借助数值计算却能够较容易获得结果。通过图7(a)可以看出,水舌进入下游水垫后,水舌的后缘区形成了一个顺时针漩涡,靠近水舌的前缘区域形成了一个较小的逆时针漩涡,并在水垫内靠左岸方向形成了另一个较大的逆时针漩涡;同时也可以看出,射流进入水垫后,主流流速沿水深迅速衰减,在到达冲坑高程189m时流速降到4~7m/s。在图7(b)中由于左侧水舌出挑角度较小,导致水舌几乎成水平方向射出且入水角度很小,所以在水舌的后缘区只形成了一个较小的顺时针漩涡;当射流进入水垫后主流被分成两股:一部分射流在水面掺气翻滚,流速逐渐衰弱;另一部分射流潜入水垫中,在未到达左岸时流速已衰减为4~5m/s。从水舌右侧紊动能k剖面图8(a)上看出:水舌在空中阶段紊动剧烈,在快落入水垫时紊动能逐渐衰弱;水舌落入下游水垫后,因强烈的掺气和翻滚,紊动能出现了一次小的峰值。结合流场图7(a)可以看出,该峰值位置出现在水舌前缘区域的逆时针漩涡处和后缘区的顺时针漩涡处。结合流场图7(b)和水舌左侧紊动能k剖面图8(b)可以看出:紊动能随着空中水舌的运动轨迹逐渐减弱,但是在水舌后缘顺时针漩涡位置和水舌主流分层处紊动能依然达到了一次小的峰值,说明此处水体紊动比较剧烈,有利于下游河道水垫内水流的进一步消能。图9分别为水舌右侧和左侧剖面上紊动能耗散率ε的分布,结合紊动能κ的分布规律可以看出:水舌在分裂段和破碎段由于空气的阻滞作用和大量空气的卷入,水舌运动变得极不稳定,因此水舌在空中阶段已经消耗掉一部分能量;当水舌进入下游河道水垫时,由于空气的卷入和水流的强烈紊动,大部分能量又在河道水垫中消耗掉,这与实际情况是吻合的。总体上,本文所建立的模型计算结果可定性地反映出水舌及下游水垫中能量耗损的分布规律,为模型试验和原型观测提供了一定的补充。4.4流速分布对比分析挑流消能将泄流的绝大部分能量消散在下游冲坑形成的水垫中,由于潜入河底的高速主流会对河床和岸坡基岩造成很大的破坏,因此水舌入水位置及下游岸边流速也是工程中关注的重点问题。本文分别选取了下游河道高程为203m、195m和189m剖面的流速分布进行了对比分析。从图10(a)中的流速分布情况可以大致看出水舌的入水位置,因为图中高程接近河道水面高程,所以水舌前缘逆时针漩涡位置出现了一个低流速区,但是紧接着高速主流又重新出现。由于右侧水舌入水角度较大,从图10(b)依然可以分辨出水舌主流入水后的位置,但是主流的位置已经逐渐偏向左岸并且流速已经有所降低,而靠近水舌左侧的流速已经降到了4~8m/s。图10(c)表明:水舌入水后的主流已在水垫中消散,在原先主流的位置流速降到4~7m/s,与试验冲坑数据对比发现最低冲坑位置就在该流速范围所在区域。通过数值模拟结果分析与模型试验观察到的下游河道流态基本吻合,说明数值模拟结果是可靠的。5优化网格结构研究结果表明,本文的数值模拟方法能够较好地模拟溢洪道或者泄洪洞出口挑流水舌的形状及获得水舌挑距、水舌入水角度及入水速度等相关水力参数。通过获得的水舌入水角度及入水速度,结合下游河道水深、河床地质条件及抗冲流速的分析,可以计算出冲坑深度和最低冲坑所在位置,为泄水建筑物体型优化设计提供依据,精度能够满足设计要求。该方法可有效用于方案的比选,在一定程度上可以替代物理模型试验,在需要通过模型试验进行验证的情况下,也可以先通过数模进行多方案筛选,使试验方案更有针对性,从而大幅度提高研究效率。值得说明的是,为了提高计算精度,以往通常的做法是先把重点模拟区域锁定在一个规则几何体内,然后再把该几何体划分成更加细密的网格,但是这种方法往往同样对非重点关注区域进行了网格加密,浪费了计算资源,不利于计算效率的提高。本文研究表明,在计算机计算能力允许的情况下,可在