平底突擴三元水躍流場的數(shù)值模擬

寧利中，張 珂，寧碧波，渠亞偉，田偉利

(1.西安理工大學(xué)省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，西安 710048；2．嘉興學(xué)院 建筑工程學(xué)院，浙江 嘉興 314001；3．上海大學(xué) 建筑系，上海 200444)

0 引 言

水利樞紐的建設(shè)改變了原河流中水流的流態(tài)，當(dāng)水流經(jīng)過水閘時，下泄水流攜帶著較大的機械能，常伴有空化、空蝕、沖刷及霧化等不良的水力學(xué)問題[1]。消能工能夠起到促使流體通過發(fā)生紊亂而以熱能的形式來減少水流機械能和保護下游河床的作用[2-3]。所以選擇適宜的消能方式和消能工的合理布置是水利樞紐能夠安全運行的保障[4-5]。

底流消能作為傳統(tǒng)消能形式之一，其優(yōu)勢明顯。底流消能具有對地質(zhì)條件的適應(yīng)性較高，對河床沖刷小，泄洪霧化很小，消能率較高等優(yōu)點[6]。目前國內(nèi)的一些水利樞紐采用了底流消能形式，如亭子口水利樞紐[7]的5個底孔、恰甫其海水電站[8]的中孔泄洪洞和導(dǎo)流洞、還有葛洲壩二江泄水閘[9]，國外具有代表性的工程有特里水電站大壩[10]、薩揚舒申斯克水電站[11]等。

在排水行洪期，多孔閘門全部啟用，可能會產(chǎn)生經(jīng)典二元水躍，但在多數(shù)情況下，為部分閘門對稱開啟，閘后水流在突擴處開始向兩側(cè)擴散，水流會互相混摻碰撞，形成一種水面波動大，水面線難以確定的復(fù)雜三元突擴水躍。閘后水流在經(jīng)過突擴斷面時流速減小，會對上游水流的下泄產(chǎn)生一定的阻礙和頂托作用，從而使得水面線壅高，當(dāng)上下游水深達到共軛水深關(guān)系時，就會產(chǎn)生水躍現(xiàn)象，水利工程上常利用水躍進行消能，即底流消能。國內(nèi)外對棱柱體渠道內(nèi)的二元水躍和漸擴式渠道內(nèi)的三元水躍已做了大量研究，羅永欽[12]結(jié)合物理模型試驗和三維數(shù)值模擬，得到了適宜觀音巖水電站突跌漸擴式溢流壩消力池體型。一些學(xué)者指出，渠道的突然收縮或擴散能夠使水流在一定范圍內(nèi)劇烈紊動，可以作為消能的主要手段[13]。但是對于平底突擴三元水躍的研究相對較少。韓偉濛等[14]以某一水利工程為原型設(shè)計并進行了物理模型試驗，對水躍躍尾水深的計算公式進行了修正。文獻[15-19]對平底突擴式渠道內(nèi)的回流區(qū)水深提出不同假設(shè)，通過動量方程推導(dǎo)出更高精度的共軛水深計算公式。對突擴水躍的分類有著不同的標準，Rajaratnam等[20]以下游水位與產(chǎn)生水躍的形式將水躍分為R水躍和S水躍。Bremen等[21]以水躍發(fā)生的位置將平底突擴渠道中的水躍分成了R型水躍、S型水躍、T型水躍和經(jīng)典水躍。文獻[22-24]對水躍的紊流特性，摻氣特性進行了測量與直接數(shù)值模擬。本文基于前人所做試驗成果，對平底突擴式渠道中發(fā)生的S型淹沒式水躍，S型臨界式水躍和R型水躍進行了三維數(shù)值模擬計算，旨在總結(jié)平底突擴三元水躍的流速、流線、紊動能和最大流速等水力特性。

1 數(shù)學(xué)模型

1.1 控制方程

連續(xù)方程為

(1)

動量方程為

(2)

雷諾應(yīng)力方程為

(3)

較標準k-ε模型而言，RNGk-ε模型能夠更好地模擬這種對流線曲度較大的流動，本文采用RNGk-ε來封閉流體力學(xué)方程。

湍動動能k方程為

(4)

湍動動能耗散率ε方程為

(5)

聯(lián)合以上方程組，針對實際問題給定相對應(yīng)的邊界條件，便可定量求解。

1.2 VOF法

VOF法的基本思想是:定義一個表示計算區(qū)域內(nèi)流體體積與計算區(qū)域總體積的相對比例的體積率函數(shù)Fw(xi,t)=1。對于某一具體計算單元而言，當(dāng)函數(shù)值為1時，表示單元完全被液體充滿；當(dāng)函數(shù)值為0時，表示一個空單元；若函數(shù)值介于0和1之間時，表示液體占單元的一部分。

描述的控制方程為

(6)

引入VOF后，ρ，μ為Fw(xi,t)的函數(shù)，其表達式為

ρ=Fwρw(1-Fw)ρa

(7)

μ=Fwμw(1-Fw)μa

(8)

式中：ρw和ρa分別為水和空氣的密度；μw和μa分別為水和空氣的分子粘性系數(shù)。

1.3 計算模型

基于盧士強[25]試驗，建立三維平底突擴水躍計算模型(圖1)，模型中突擴口至上游始端的距離和下游河道末端的長度分別為0.50 m和3.00 m，上下游側(cè)墻高分別為0.50 m和0.40 m，上、下游河寬分別為0.15 m和0.30 m，即擴散比為2.0。模型主要分為閘室段和下游擴散段。突擴口在計算模型的位置為x=0.00 m(yz平面)，上游x為負值，下游x為正值。

圖1 三維平底突擴水躍計算模型(單位：m)Fig.1 Calculation model of 3D hydraulic jump with flat bottom and abrupt lateral expansion

1.4 計算工況

對擴散比為2.0下的體型做了數(shù)值模擬計算。已知條件：閘門開啟高度為0.03 m，上游水位為0.44 m，躍前水深為0.025 m，躍前流速為2.53 m/s，躍長為0.90 m。再通過改變下游水位(h2=0.09，0.13，0.17 m)，使得水躍逐漸由R型水躍過渡到S型臨界式水躍，再變化到S型淹沒式水躍，從而實現(xiàn)水躍過程的變化。

1.5 邊界條件

邊界條件如下：進口為流速進口，在計算時給定上游水位0.44 m，并給定一個行進流速；出口為壓力出口，給定一個下游水位(h2=0.09，0.13，0.17 m)，并且下游出口處的水流為緩流，服從壓力靜態(tài)分布；模型上部空氣為壓力進口，進口斷面最頂部工作壓強設(shè)為101 325 Pa；側(cè)面和底面為無滑移固體壁面。

1.6 計算方法

采用VOF法進行自由表面的追蹤計算，流體力學(xué)氣液兩相流時均方程采用RNG模型進行封閉，控制方程采用有限體積法進行離散，速度與壓力耦合方程采用PIOS算法進行求解。時間步長為0.01 s，收斂標準為進出口質(zhì)量流通率的差值與進口質(zhì)量流通率之比≤1%。

利用Fluent軟件，采用非定常流動PISO算法，對平底突擴三元水躍進行了數(shù)值模擬計算。

2 計算結(jié)果及分析

2.1 計算模型的驗證

高遠[26]以盧士強[25]所做試驗數(shù)據(jù)為基礎(chǔ)，利用FLOW-3D軟件對平底突擴三元水躍進行了數(shù)值模擬計算。本文對閘門開啟高度為0.03 m，上游水位為0.44 m的工況時S型臨界式水躍，利用Fluent軟件進行了數(shù)值模擬，提取計算了x=0.10 m，x=0.65 m，x=1.10 m位置橫截面中垂線的流速平均值的計算結(jié)果并與高遠所做數(shù)值模擬計算結(jié)果進行對比，對比結(jié)果見表1。

表1 不同橫截面的平均流速對比Table 1 Comparison of average velocity for different cross sections

由表1可見，3個位置橫截面中垂線的平均流速值相近。另外，本文計算的S型臨界式水躍對應(yīng)的下游水位為0.13 m，躍尾位置在x=0.90 m附近，并經(jīng)計算得到躍尾斷面的平均流速為0.43 m/s，與盧士強所做試驗結(jié)果也相吻合?？沈炞C本文建立的計算模型是正確的。

2.2 計算結(jié)果及分析

2.2.1 流線分布

繪制3個工況下中軸立面對應(yīng)的流線圖，見圖2。由圖2可見，水躍發(fā)生的躍首位置。在沿程方向，3個工況下水流均發(fā)生了不同程度的水平旋滾，但旋滾的位置不同。流線的疏密情況反映出流速的大小，流線密的區(qū)域?qū)?yīng)的流速大。水流經(jīng)過突擴斷面后，由急流向緩流轉(zhuǎn)變的過程中發(fā)生了水躍，在水躍發(fā)生區(qū)域，水深沿程逐漸增大，在水流的上部區(qū)域形成逆時針方向的翻滾漩渦，空氣進入漩渦中，使得主流與旋滾的水流發(fā)生碰撞，從而實現(xiàn)能量交換與損失，該區(qū)域為水躍的主要消能區(qū)。

S型臨界式水躍(h2=0.13 m)z=0.024 m平面的流線圖見圖3。由圖3可見，水流經(jīng)過突擴斷面后向兩側(cè)發(fā)生了幾乎對稱的擴散運動，在靠近突擴口處，左右兩側(cè)水流形成了回流，主要發(fā)生在x=0 m斷面至x=0.20 m斷面之間，回流區(qū)的范圍不大。擴散后的水流經(jīng)過渠道側(cè)壁的阻礙作用，折沖擠壓主流，主流的流線加密，主流的流速有所增加，在x=1.00 m斷面至x=1.30 m斷面形成了漩渦，漩渦在左右兩岸基本對稱分布。從x=1.50 m斷面開始，水流平穩(wěn)。這里在x=1.00 m斷面至x=1.30 m斷面之間漩渦的不對稱可能是主流擺動引起的結(jié)果。

2.2.2 紊動能分布

3個工況下中軸立面紊動能分布見圖4。由圖4可見，在水流表面均形成了一個有空氣摻入的翻騰回旋滾動，水流內(nèi)部產(chǎn)生強烈的摩擦與摻混作用。紊動能表示流體紊動的程度，流體之間的剪切、摩擦混摻、碰撞作用越劇烈，紊動能值越大。紊動能的最大值均在水躍中部，且紊動能由旋滾區(qū)的內(nèi)部向外部逐漸減小。R型水躍發(fā)生紊動的范圍主要在從突擴口向下游x=1.00 m，z=0.10 m內(nèi)，最大值為0.35 m2/s2；S型臨界式水躍發(fā)生紊動的范圍較大，主要是在從突擴口向下游x=1.00 m，z=0.16 m內(nèi)，最大值為0.20 m2/s2；而S型淹沒式水躍發(fā)生紊動的范圍最小，主要是在突擴口向下游x=0.80 m，z=0.10 m內(nèi)，最大值為0.20 m2/s2。

R型水躍和S型淹沒式水躍相比，前者發(fā)生的紊動能較大，但水躍發(fā)生的范圍較大；S型臨界式水躍和S型淹沒式水躍相比，兩者的紊動能相近，但后者水躍發(fā)生的范圍較小。

2.2.3 流速分布

3個工況下中軸立面流速分布沿程變化圖見圖5。由圖5(a)可見，流速在水躍區(qū)各斷面從底板開始逐漸增大到最大值Umax，然后逐漸減小；由圖5(b)可見，在x=0.60 m斷面之前，流速在水躍區(qū)各斷面從底板開始逐漸增大到最大值Umax，然后逐漸減小到零，之后向上游反方向增加；在x=0.90 m斷面，流速在水躍區(qū)各斷面從底板開始逐漸增大到最大值Umax，然后逐漸減小；由圖5(c)可見，在x=0.90 m斷面之前，流速在水躍區(qū)各斷面從底板開始逐漸增大到最大值Umax，然后逐漸減小到零，之后向上反方向增加；在x=1.20 m斷面，流速在水躍區(qū)各斷面從底板開始逐漸增大到最大值Umax，然后逐漸減小。

圖5 3個工況下中軸立面流速沿x軸變化Fig.5 Velocity distribution of vertical section of central axis along x-axis under three working conditions

按照邊界層理論，邊界層厚度δ為各斷面上流速最大值至底板的間距。圖5(a)中，在x=1.20 m斷面之前，邊界層厚度δ沿程變化不大，大致在0.01 m高度，在x=1.20 m斷面之后，邊界層厚度δ逐漸增大到水面附近。圖5(b)中，在x=1.20 m斷面之前，邊界層厚度δ沿程略有增加，但增加趨勢平緩，在x=0.90 m斷面之后，界層厚度δ逐漸增大到水面附近。圖5(c)中，在x=1.20 m斷面之前，邊界層厚度δ沿程逐漸增加，增加趨勢明顯，在x=1.20 m斷面之后，邊界層厚度δ逐漸增大到水面附近。

橫斷面上流速分布規(guī)律為：R型水躍，水躍區(qū)內(nèi)，邊界層厚度的變化不大，之后在水躍區(qū)向明渠均勻流過渡的過程中變化明顯，且逐漸增大到水面附近；S型臨界式水躍區(qū)內(nèi)，邊界層厚度的逐漸增大，但增加趨勢不明顯，之后水躍區(qū)向明渠均勻流過渡的過程中變化明顯，且逐漸增大到水面附近；S型淹沒式水躍區(qū)內(nèi)，邊界層厚度的逐漸增大，且增加趨勢明顯，之后水躍區(qū)向明渠均勻流過渡的過程中逐漸增大到水面附近。

提取3個工況下的z=0.012 m平面x=0.10 m和x=0.64 m 2個位置的流速分布見圖6。

由圖6(a)可見，R型水躍的流速整體上最大，S型臨界式水躍次之，S型淹沒式水躍最小，但流速分布規(guī)律相同，因為側(cè)壁為無滑移壁面，在壁面處流速為零。在靠近兩邊側(cè)墻附近的一定距離內(nèi)為副流區(qū)，其流速比主流區(qū)要小，呈現(xiàn)出向上游凹陷分布；由圖6(b)可見，R型水躍的流速整體上依然最大，S型臨界式水躍次之，S型淹沒式水躍最小，但S型淹沒式水躍的流速分布趨于均勻。3個工況下副流區(qū)范圍都有所增大，主流區(qū)減小。

圖6 3個工況下不同x時在z=0.012 m平面上的流速分布Fig.6 Velocity distribution of z=0.012 m plane at different x under three working conditions

由此得出，R型水躍流速最大，攜帶能量大，對河床的沖刷最為嚴重，躍尾位置距突擴口的距離最遠；S型臨界式水躍流速較小，攜帶能量較小，流態(tài)不穩(wěn)定，易發(fā)生折沖水流，躍尾位置距突擴口的距離減??；S型淹沒式水躍流速最小，攜帶能量最小，流態(tài)較為穩(wěn)定，對河床的沖刷最小，躍尾位置距突擴口的距離最小。

3個工況下橫斷面最大流速沿程變化見圖7。由圖7可見，3個工況下，橫斷面最大流速沿程逐漸減小，不同之處在于：R型水躍，在x=1.50 m處最大流速值變化開始減小，減小到0.8 m/s左右時最大流速趨于平穩(wěn)；S型臨界式水躍，在x=1.20 m處最大流速值變化開始減小，減小到0.6 m/s左右時最大流速趨于平穩(wěn)；S型淹沒式水躍，在x=1.00 m處最大流速值變化開始減小，減小到0.5 m/s左右時最大流速趨于平穩(wěn)。由此得出：S型淹沒式水躍在較短的區(qū)域內(nèi)，最大流速值以較快速率減小，并且區(qū)域平穩(wěn)時的最大流速值相比另外兩種工況下要小。

圖7 3個工況下橫斷面最大流速沿x軸變化Fig.7 Variation of maximum velocity of the cross section along x-axis under three working conditions

3 結(jié) 論

1)通過對比前人所做數(shù)值模擬結(jié)果，驗證了本文建立的平底突擴三元水躍數(shù)值模擬計算模型的正確性，VOF法能夠較為準確地跟蹤計算自由水面，RNGk-ε湍流模型能夠得到在消力池發(fā)生水躍的瞬態(tài)過程，能夠得到一定精度的流速、紊動能等水力特性，表明采用Fluent軟件對平底突擴水躍進行三維數(shù)值模擬計算具有較強的適用性。

2)得到3種水躍紊動能的分布特點，發(fā)現(xiàn)水流紊動能最大值在水躍中部，R型水躍發(fā)生的紊動能較大，發(fā)生紊動的范圍主要在從突擴口向下游x=1.00 m，z=0.10 m內(nèi)，最大值為0.35 m2/s2；S型臨界式水躍和S型淹沒式水躍兩者的紊動能相同，最大值均為0.20 m2/s2，但是S型臨界式水躍發(fā)生紊動的范圍主要是在從突擴口向下游x=1.00 m，z=0.16 m內(nèi)，后者主要是在突擴口向下游x=0.80 m，z=0.10 m內(nèi)。

3)得到3個工況下中軸立面流速分布及邊界層變化的特點，發(fā)現(xiàn)在相同位置的橫斷面上，R型水躍流速均較大，R型水躍邊界層變化不大，S型臨界式水躍邊界層增長緩慢，S型淹沒式水躍邊界層增加趨勢明顯，最終均增大到水面附近。

4)得到3個工況下z=0.012 m平面沿程流速分布，R型水躍流速最大，躍尾位置距突擴口的距離最遠；S型臨界式水躍流速較小，流態(tài)不穩(wěn)定，易發(fā)生折沖水流，躍尾位置距突擴口的距離減小；S型淹沒式水躍流速最小，流態(tài)較為穩(wěn)定，躍尾位置距突擴口的距離最小。

5)S型淹沒式水躍最大流速值以較快速率減小，并且最大流速值相比另外兩種工況下要小。