趙東陽，尹進(jìn)步，王國杰，楊 彥，楊 釗，張曙光

(西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，陜西 楊凌712100)

0 引 言

高水頭、大流量的高壩泄洪消能問題，往往是大型水利水電樞紐工程設(shè)計中重大而復(fù)雜的科學(xué)技術(shù)問題。需要綜合考慮擋水建筑物、泄水建筑物和興利建筑物三者之間的關(guān)系和合理布局，并根據(jù)壩址的地形、地質(zhì)和水文條件，合理選取泄洪方式、泄水建筑物的型式和消能方式[1]。同樣，本次研究借助的高水頭岸邊溢洪道之所以采用底流，是由于鼻坎位于岸邊，距離主河床較遠(yuǎn)，且水舌挑距有限，無法落到河床，故只能放棄挑流采用底流。原方案消力池邊墻不擴(kuò)散，通過模型試驗發(fā)現(xiàn)由于躍后水深不足，消力池內(nèi)未能形成淹沒水躍，為增加消能水體，采用邊墻擴(kuò)散來增加消力池尺寸。由于該消力池右側(cè)邊墻為山體且坡度較陡，開挖工程量和邊坡支護(hù)量較大，而左側(cè)邊墻在河道內(nèi)，開挖量不大，因此采用單側(cè)漸擴(kuò)。

目前使用數(shù)值模擬分析消力池內(nèi)水力特性的研究較多，Iwao，Katakam研究發(fā)現(xiàn)對稱擴(kuò)散消力池，消能率較高，流態(tài)較好[2,3]。王月華利用數(shù)值模擬提供的消力池Fr值、紊動能分布圖和流速矢量圖，比較全面地反映了三級消能池水流和消能情況[4]。鄭梅玲通過數(shù)值模擬分析了消力池內(nèi)部流場，發(fā)現(xiàn)漸擴(kuò)式能使水流橫向擴(kuò)散的空間增大，增強(qiáng)消力池首部水流的立軸旋滾[5]。上述只是關(guān)于對稱擴(kuò)散消力池的研究，而針對單側(cè)漸擴(kuò)消力池的研究較少，本文參考相關(guān)工程資料[6-11]，采用Flow-3d軟件中RNGk-ε雙方程紊流模型基于VOF方法，對某工程單側(cè)漸擴(kuò)消力池進(jìn)行數(shù)值模擬，得到水面線、流速、壓強(qiáng)和紊動能等水力參數(shù)，并與模型試驗結(jié)果進(jìn)行對比分析，為工程設(shè)計提供參考。

1 數(shù)學(xué)方法及模型建立

1.1 控制方程

采用RNGk-ε雙方程湍流模型對消力池進(jìn)行數(shù)值模擬，連續(xù)性方程，動量方程和k-ε方程表示如下：

連續(xù)性方程：

(1)

動量方程：

(2)

紊動能k方程：

耗散率ε方程：

(4)

1.2 VOF方法

采用VOF方法[12]追蹤水流自由表面，該方法定義流體體積函數(shù)F=F(x,y,z,t)表示計算區(qū)域內(nèi)流體的體積占據(jù)計算區(qū)域的相對比例。對于某個單元，F(xiàn)=1表示該單元被流體完全充滿，F(xiàn)=0表示該單元是個空單元，沒有流體，F(xiàn)=0～1表示該單元被流體部分充滿，VOF中自由表面的復(fù)雜變化可用函數(shù)描述為：

(5)

式中：Ax、Ay、Az分別為x、y、z三個方向可流動的面積分?jǐn)?shù)；u、v、w為流速；VF為可流動的體積分?jǐn)?shù)。

1.3 模型建立

計算區(qū)域包括進(jìn)水渠、控制段、泄槽和消力池及部分海漫，總長323.00 m(見圖1)。溢洪道為岸邊開敞式有閘控制正槽溢洪道，控制段設(shè)置兩孔閘門，泄槽陡坡段設(shè)有兩道摻氣坎。溢洪道出口采用底流消能，消力池采用單側(cè)漸擴(kuò)式，底板高程為834.00 m，進(jìn)口寬14.00 m，出口寬20.37 m，長84.53 m，消力池出口采用直立尾坎，為增加躍后水深，在直立尾坎上設(shè)置高2.50 m、寬1.50 m不連續(xù)的梯形墩，邊墻擴(kuò)散角度為5.6°。查閱資料得知，實際工程擴(kuò)散角度一般較小，通常小于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、溢洪道右側(cè)邊墻外立面為Y=0 m(右側(cè)墻內(nèi)壁為Y=2.50 m位置)、消力池底板下部6.10 m高程的水平面為Z=0 m，即消力池底板高度為Z=6.10 m。計算區(qū)域全部用結(jié)構(gòu)化網(wǎng)格來劃分，單元尺寸為 ，網(wǎng)格總數(shù)大約790萬。

圖1 溢洪道結(jié)構(gòu)(單位：m)Fig.1 Layout of the spillway

邊界條件：計算區(qū)域上下游為壓力邊界，并給定水位；底部及左右兩側(cè)均為壁面邊界；頂部為壓力邊界。壓力求解器選用廣義極小殘差算法(GMRES)，基于壓力隱式求解法(Implicit)；計算結(jié)果數(shù)據(jù)輸出間隔2 s，輸出數(shù)據(jù)有流速、壓力、水的體積分?jǐn)?shù)、自由液面高程等；初始時間步長定為0.000 1 s。

2 計算與試驗對比

為確保計算結(jié)果的可靠性，此處分別將消能防沖工況和校核工況下計算得到的水面線、流速、壓強(qiáng)等水力參數(shù)與相應(yīng)工況的模型試驗數(shù)據(jù)進(jìn)行對比。用于對比的試驗?zāi)Ｐ桶粗亓ο嗨茰?zhǔn)則設(shè)計，模型比尺為 1∶40。

2.1 水面線

兩工況下計算水面線與試驗水面線對比見圖2。從圖2(a)可以看出，消能防沖工況試驗水面線與計算水面線吻合較好，在樁號0+224.26 m兩者重合，相對誤差為0%。在樁號0+280.00 m處，兩者相對誤差最大，為12.30%，但考慮消力池內(nèi)水流翻滾劇烈，水面波動較大，仍認(rèn)為具有一定的可信度。

圖2 水面線對比Fig 2.Comparison of surface line

校核工況試驗水面線與計算水面線對比由圖2(b)可見，兩者相差甚微，最大相對誤差出現(xiàn)在樁號0+290.00 m處，僅為7%。

對比兩個工況計算結(jié)果發(fā)現(xiàn)，隨著流量的加大，水面爬升梯度越大。由于此消力池邊墻擴(kuò)散角度較小，僅5.6°，計算擴(kuò)散側(cè)水面線(見圖2)與計算中心線位置水面線對比發(fā)現(xiàn)兩者基本重疊，故而無必要對漸擴(kuò)側(cè)水面線進(jìn)行研究。

2.2 流 速

此處首先對計算得出的中心線處底、表流速與相應(yīng)的試驗值進(jìn)行對比如圖3所示，可見兩者相差無幾，在此基礎(chǔ)上，對計算中心線位置沿程流速分析之后，又詳細(xì)分析了水平面流場。

圖3 消力池流速對比Fig.3 Comparison of stilling basin velocity

消能防沖工況下沿程流速分布見圖3(a)所示，可見在樁號0+253.00 m以前，底部正向流速較大，超過10 m/s，在樁號0+268.50 m以前，表面流速為負(fù)值，說明該區(qū)域表面水流流向為回流。在樁號0+259.00 m以前，表面反向流速較大，超過3 m/s。

分析校核工況下沿程流速分布[見圖3(b)]發(fā)現(xiàn)，表面流速在樁號0+279.50 m以前為負(fù)值，水流流向為回流。反向流速在樁號0+270.00 m以前流速較大，超過3 m/s。底部正向流速在樁號0+270.00 m以前流速較大，超過10 m/s，對兩個工況計算與測試結(jié)果對比發(fā)現(xiàn)，在樁號0+224.50～270.00 m區(qū)域，兩工況底部正向流速超過10 m/s，表面反向流速超過3 m/s，對實際工程而言，該正向流速與反向流速可能偏大，容易形成底板磨蝕破壞。

試驗發(fā)現(xiàn)擴(kuò)散側(cè)表面時而存在立軸旋渦，消能防沖工況時在樁號0+278.00 m位置[見圖4(a)]，校核工況時在樁號0+289.00 m[見圖4(b)]。然而消力池內(nèi)流態(tài)復(fù)雜，水體來回波動頻繁，試驗測量池內(nèi)水面及以下流場不便，為了對池內(nèi)流場分布規(guī)律進(jìn)行深入探索，下面采用數(shù)值模擬對水平面流場進(jìn)行分析。

2.3 水平面流場

平面流場能清楚反映水流的運動軌跡，為準(zhǔn)確捕捉流場的運動過程，每個工況分別均截取三個高程的水平面流場，如圖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 消力池流態(tài)圖Fig.4 Flow pattern of still basin

圖5 消力池內(nèi)不同高程流場分布圖Fig.5 Distribution of flow in the still basin

消能防沖工況下水平面流場見圖5(a)，Z=6.20 m截面接近消力池底板，漸擴(kuò)側(cè)尾部的立軸旋渦渦心在樁號0+280.00 m(與試驗測得樁號0+278 m接近)，范圍大約在X=281.4～290.3 m，Y=19.26～20 m，漩渦范圍較小。由Z=10.20 m截面可見左邊漸擴(kuò)側(cè)在樁號0+221.75 m處產(chǎn)生翻滾漩渦，即躍首位置，同時漸擴(kuò)側(cè)尾部的立軸旋渦基本消失。由Z=18.73 m水平面流場可見，旋滾區(qū)發(fā)展到0+275.00 m停止，即躍尾位置。

校核工況見圖5(b)，由Z=6.20 m流場分布圖可見，在消力池漸擴(kuò)側(cè)尾部的立軸旋渦渦心樁號為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處，位于左邊漸擴(kuò)側(cè)的漩渦范圍縮小。由Z=18.73 m得出，躍尾在樁號0+296.00 m位置。躍尾橫軸在漸擴(kuò)側(cè)向下游傾斜的角度較消能防沖工況更大。

對比兩個工況水平面流場發(fā)現(xiàn)，漸擴(kuò)側(cè)尾部存在立軸旋渦，靠近消力池底板旋渦范圍最大，隨高程升高逐漸減小，且隨入池流量加大，旋渦范圍越大。躍首橫軸在漸擴(kuò)側(cè)向上游移動，躍尾橫軸在漸擴(kuò)側(cè)向下游移動。分析其原因，水流受到右側(cè)邊墻約束，只能向左邊約束較小的漸擴(kuò)側(cè)擴(kuò)散，由此在底板附近產(chǎn)生的立軸旋渦將躍尾橫軸推向下游，將躍首橫軸推向上游。

2.4 時均壓強(qiáng)

試驗中消力池中心線處時均壓強(qiáng)采用玻璃管測試，由試驗結(jié)果與計算結(jié)果對比圖圖6可見，兩工況在樁號0+260.00 m以前計算與試驗結(jié)果相差甚微，最大相對誤差不足5%，在樁號0+260.00 m以后誤差稍大，消能防沖工況相對誤差小于10%，校核工況相對誤差最大值僅為5%，計算與試驗結(jié)果基本吻合。

各工況下時均壓強(qiáng)基本均呈現(xiàn)先增后減再增之后逐步趨于穩(wěn)定，消力池內(nèi)均未出現(xiàn)負(fù)壓。樁號0+221.50～232.30 m區(qū)域為反弧段，受離心力作用，水流劇烈碰撞消力池底板，造成局部壓強(qiáng)迅速增大。之后離心力減小，壓強(qiáng)逐漸下降。在樁號0+242.50 m位置，校核工況比消能防沖工況的上游水頭大，弗汝德數(shù)大，躍前水深小，即靜水壓強(qiáng)小，所以該位置校核工況時均壓強(qiáng)值比消能防沖工況時小。該部位流速較大，壓強(qiáng)較低，在實際工程中，摻氣量不足的情況下，較容易破壞，而本工程中設(shè)有兩道摻氣坎，且水面較薄，摻氣量足夠，因而發(fā)生破壞的可能性較小。樁號0+242.50～270.00 m區(qū)域上部為水躍，水面迅速爬升，且水躍有向下游的推力，因此壓強(qiáng)回升較快。樁號0+270.00 m之后水面逐漸平穩(wěn)，壓強(qiáng)也隨之平穩(wěn)。

擴(kuò)散側(cè)壓強(qiáng)如圖6所示，計算中心線位置與計算擴(kuò)散側(cè)水面線基本重疊，在消能防沖工況下，兩者相對誤差最大僅為8‰，校核工況下最大僅差2%，差別甚微，因此文中不再對擴(kuò)散側(cè)壓強(qiáng)進(jìn)行分析。

圖6 消力池壓強(qiáng)對比Fig.6 Comparison pressure stilling basin

2.5 探討計算紊動能與試驗脈動壓強(qiáng)

紊動能是能夠體現(xiàn)水流由紊動而產(chǎn)生的能量耗散程度的重要指標(biāo)，其可以衡量流體的紊動狀態(tài)，對考察水躍消能的防沖建筑物性能具有重要意義。本文采用的 RNGk-ε雙方程紊流模型可以計算水流的紊動能：

(6)

式中：u′為脈動流速。

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

由消力池底板紊動能云圖圖7可見，紊動能向漸擴(kuò)側(cè)擴(kuò)展明顯，隨水位流量加大，范圍逐漸向消力池中心線移動。因此消力池重點防護(hù)區(qū)域偏向漸擴(kuò)側(cè)。

圖8 紊動能縱剖面圖 Fig.8 Turbulent energy longitudinal section

圖9 消力池底板脈動壓力Fig.9 Stilling basin floor pulsation pressure

取紊動能最靠近上游位置縱剖面(見圖8)，與試驗中相應(yīng)位置脈動壓力(見圖9)進(jìn)行對比發(fā)現(xiàn)，消能防沖工況下，在樁號0+225.00 m處脈動壓力最大，脈動壓力均方根最大值為34.80 kPa，對應(yīng)到紊動能云圖上為紅色區(qū)域的起始位置[見圖8(a)]。校核工況下，脈動壓力在樁號0+238.00～242.50 m范圍最大，脈動壓力均方根最大值為47.20 kPa，對應(yīng)到紊動能云圖同樣為紅色區(qū)域的起始位置[見圖8(b)]。由此可見紊動能云圖可能反映脈動壓力分布情況，分析其原因，紊流由大小不同尺度的渦體(或漩渦)所組成，尺寸各異的漩渦在高速水流中旋轉(zhuǎn)，轉(zhuǎn)向與水流方向一致的一邊流速增大，因而壓強(qiáng)變小，另一邊則流速減小壓強(qiáng)增大；兩邊的壓強(qiáng)差產(chǎn)生的推力，對于渦體而言就是升力。當(dāng)升力足以克服水流的黏滯阻力時，渦體將上升做橫向運動，進(jìn)入流速不同的流層，即產(chǎn)生水流質(zhì)點的混摻作用。由于紊流中有無數(shù)個轉(zhuǎn)向與大小不同的渦體，所以水流的混摻作用一直雜亂無章地進(jìn)行著，從而導(dǎo)致了水流流速的脈動；根據(jù)能量方程就導(dǎo)致了壓強(qiáng)的脈動，兩者直接關(guān)聯(lián)。所以脈動壓力必將受到脈動流速影響，紊動能 在一定程度上反映脈動壓力分布情況。

3 結(jié) 論

本文采用Flow-3d軟件對單側(cè)漸擴(kuò)消力池進(jìn)行數(shù)值模擬，得到了水面線、流速、壓強(qiáng)、湍動能等水力參數(shù)分布規(guī)律，并與水工模型試驗結(jié)果進(jìn)行了對比驗證，得出結(jié)論如下：

(1)Flow-3d軟件中的RNGk-ε雙方程紊流模型能較好地模擬單側(cè)漸擴(kuò)消力池中的水力特性，數(shù)值模擬結(jié)果與模型試驗結(jié)果吻合較好，為進(jìn)一步研究奠定了基礎(chǔ)。

(2)詳細(xì)分析消力池內(nèi)流場，發(fā)現(xiàn)躍首橫軸在漸擴(kuò)側(cè)向上游移動，躍尾橫軸在漸擴(kuò)側(cè)向下游移動，且隨流量的加大，移動的距離會逐漸加大。

(3)通過水平面流場分析發(fā)現(xiàn)，消力池漸擴(kuò)側(cè)尾部產(chǎn)生立軸旋渦，接近消力池底板漩渦范圍最大，向上逐漸減小，且隨流量加大，范圍將越來越大。

(4)對比計算紊動能和試驗脈動壓力得出：紊動能在一定程度上反映脈動壓力分布情況；反弧段及消力池底板前部的脈動壓力較大，設(shè)計上應(yīng)對這些區(qū)域加強(qiáng)防護(hù)。