王海軍, 張 凡, 李會平
(1.天津大學 水利工程仿真與安全國家重點實驗室, 天津 300350;2.天津大學 建筑工程學院, 天津 300350)
隨著我國高壩建設的快速發展,泄水建筑物水頭越來越高。溢洪道是水利樞紐中重要的泄洪建筑物。隨著下泄流量增大,對溢洪道的過流能力及結構的安全性要求越來越高。溢洪道的過流能力及其壁面壓強、臨底流速的分布情況等因素都對溢洪道泄水、消能有重要的影響[1]。臨底流速是底流消能工的一項重要水力學指標,同時也是影響結構安全的一項重要因素。傳統的溢洪道消力池水力特性分析主要是采用物理模型試驗法[2]以及原型觀測法,尤其是前者,在已做的工程設計和水力靜動力特性研究中起到了極為重要的作用。楊敏等[1]利用物理模型試驗的方法對跌坎消力池的臨底流速進行了探究,分析比較了跌坎消力池與傳統消力池的臨底流速分布規律;王立杰等[3]采用物理模型試驗的方法針對新型射流消能的流態演變與極限臨底流速進行了研究,同時也發現,在采用物理模型試驗時,模型建立較為復雜、成本較高、測量精準性難以保證,且存在比較大的簡化,比尺效應影響較大。隨著近年來計算機技術的發展及數值計算方法的完善,采用數值模擬方法分析研究復雜水流也成為經濟有效的手段。基于計算流體動力學的大型商業軟件 Fluent、Flow-3D等的開發利用,水利領域的數值仿真模擬也得到快速發展[4]。李玲等[5]研究了三維VOF模型及其在溢洪道水流計算中的應用;薛宏程等[6]對溢洪道出口斜切型挑坎挑射水舌進行了三維數值模擬;王青等[7]基于Flow-3D對陡坡彎道水流進行了三維數值模擬;李樹寧等[8]對不同體型參數的跌坎型消力池進行了數值模擬計算,為工程設計提供依據等。
本文采用 Fluent軟件,對某水電站溢洪道泄洪進行精細數值模擬,結合該水電站原型觀測數據,驗證所采取數值模擬方法的合理性。并在此基礎上探尋溢洪道水流臨底流速的分布規律及臨底流速梯度對溢洪道結構安全的影響。
臨底流速的數值模擬主要包括基本的流體控制方程、自由液面VOF模型、模型求解方法。
將連續性方程和不可壓縮黏性流體運動的Navier-Stokes方程作為流體運動的控制方程,紊流模型采用的是RNGk-ε模型[9]。
(1) 連續性方程
(1)
式中:ui為i方向的速度分量,m/s;xi為坐標方向;下標i為變量的不同方向,取值1、2、3。
(2) 運動方程
(2)
式中:fi為作用于單位質量水體的質量力,m/s2;ρ為流體密度,kg/m3;p為壓強,N/m2;t為時間,s;uj為j方向的速度分量,m/s;xi、xj為坐標方向(下標i、j為變量的不同方向,取值1、2、3)。
(3)k方程
(3)
式中:Gk為剪切產生項,m2/s3;Gb為浮力產生項,m2/s3;k為單位質量紊動動能,m2/s2;ε為紊動動能耗散率,m2/s3;ν為運動黏性系數,m/s2;νt為紊流運動黏性系數,m/s2。
(4)ε方程
(4)
式中:系數C1ε、C2ε、C3ε為常量;σk、σε為k方程和ε方程的Prandtl數;Rε為附加項,m2/s4,以適應瞬變流和流線彎曲的影響[10]。
VOF模型(volume of fluid)基本原理是通過研究網格單元中流體與網格體積比函數來追蹤自由液面。相比于Mixture模型與歐拉模型更加適用于計算任意氣液分界面的分界面。VOF模型中不同的流體共用一套運動方程,通過引進相體積分數這一變量,實現對每一個計算單元相界面的追蹤。在每個控制容積內,所有相體積分數額總和為1,所有變量及其屬性在所控制容積內各相共享,且代表了容積的平均值。在任何給定控制容積內的變量和其屬性代表了一相或多相的混合。假設單元內第m相流體的體積分數為αm,則存在3種情況:
(1)αm=0,單元體積內不存在第m相流體;
(2)0<αm<1,單元體積內包含了第m相流體和其他一相或多相流體,且存在著不同相流體之間的交界面;
(3)αm=1,單元體積內充滿第m相的流體。單元體積內具有相同的速度和壓力場。
在VOF模型中,追蹤第m相流體自由水面的控制微分方程為:
(5)
通過對該方程的求解可以完成自由水面的追蹤[11-14]。
流場計算采用的是有限體積法,其是一種分塊近似的計算方法。具體的將計算區域劃分為一系列不重復的控制體積,并使每個網格點周圍有一個控制體積,這樣計算區域被離散成許多小的體積單元,然后將待解的微分方程對每一個小單元體進行積分,便得出一組離散方程解。其求解方法分為分離解法和耦合解法兩大類。具體的擴展分為SIMPLE算法以及在其基礎上改進的SIMPLEC、SIMPLER、PISO算法。由于PISO算法對瞬態問題有明顯的優勢[15],故本模擬采用PISO算法。
某工程采用了溢流表孔底流消能的泄洪消能方式。泄洪建筑物的泄水孔為開敞式表孔,總共分為5孔,單孔尺寸為13 m×20 m。邊墩、中墩厚度分別為4.0和5.0 m,溢流壩段總寬為93.0 m;孔口設有檢修平板門和弧形工作門,閘墩前緣呈尖圓流線形,溢流壩堰頂高程為1 398.0 m。堰面為WES曲線,曲線末端接1∶0.75直線段和半徑R=55 m的反弧段;溢流壩后為長約200.0 m的泄槽和170.0 m的消力池,入池跌坎高6.0 m,池底高程1 276.5 m,消力池后設海漫與河床平順連接。溢洪道縱剖面示意圖如圖1所示。

圖1 溢洪道剖面示意圖
模型計算求解步驟為: (1)在Auto CAD (前處理器) 中建立幾何模型。(2)將模型導入ICEM中,進行網格的劃分。(3)將劃分完成的網格模型導入Fluent中,并設置相關參數。參數設置包括計算模型選取、材料屬性和邊界條件設定以及求解方法的選取和運行環境的設置等。(4)對流場的初始化、求解。(5)進行計算結果的后處理及結果的分析。
網格的劃分是流場數值模擬的重要環節。網格質量的高低直接影響計算精度及收斂速度。溢洪道泄流的模擬范圍大、流道體型復雜,本文研究分析的
重點是臨底流速分布,因此在網格劃分時,近底處進行了細化加密,上部網格適當放大。計算模型約525×104個節點,約498×104個網格。
模型自由表面的處理采用VOF模型,邊界條件采用壓力進出口邊界條件。模型進流邊界為壓力進口,給定進水口水頭;出流邊界為壓力出口,給定出口水頭;固壁面采用無滑移壁面,法向和切向速度均為0。
利用建立的溢洪道三維數學模型,對已進行過原型觀測的工況進行模擬計算。采用3#、4#孔全開工況(表1中工況1)的模擬結果進行驗證,此工況堰上水頭為117.15 m,下游水頭13.50 m,溢洪道流量3 336 m3/s。比較數值模擬計算值和原型觀測值,對所建立的數學模型進行驗證。
溢洪道縱剖面水氣相分布如圖2所示。從圖2可以看出水流出閘室后,沿溢流面下泄,水流在溢流堰曲線段水流平順下泄,出閘室后的水流橫向擴散。水流沿溢流堰面急速下泄至消力池內,水流在消力池內劇烈摻混,充分消能,出消力池后水流趨于平順泄向下游,消力池消能效果較好。整個過程水流在溢流堰反弧段與泄槽陡坡段水深均較小。模擬水流流態與原型觀測流態分布基本一致,原型觀測如圖3所示。
提取自由液面水深和實測水深對比,圖4為3#孔中心線處溢洪道水深模擬計算值與原型實測值比較圖(橫坐標為分析位置距離壩趾的距離),二者吻合度較好。溢流堰陡坡處水深下降迅速,至消力池后水深迅速增加,出消力池后水深開始減小。
通過流態和水面線驗證,數值仿真模擬與原型觀測結果基本吻合,可見數值仿真模擬結果可作為原型觀測數據的補充。

圖2溢洪道縱剖面水氣相分布圖3原型觀測水流流態

圖4水深計算值與觀測值對比圖圖5縱向最大臨底流速變化趨勢圖
分別選取原型觀測中出現的最大流量、最小流量以及中間流量3種典型工況進行數值模擬,從而分析臨底流速分布規律。具體工況如表1所示,分析3種模擬工況,提取3#孔中軸線的縱向最大臨底流速分布,如圖5所示。從數值仿真結果可知,縱向不同斷面臨底流速最大值順水流方向呈現先增大后減小的趨勢。同時觀察到在1-1斷面(距離壩址102 m,見圖1)和2-2斷面(距離壩址200 m,見圖1)的臨底流速最大值在工況3時分別達到了42.24和44.95 m/s,如表2所示。
圖6~7為3種工況下,3#孔中軸線在1-1斷面與2-2斷面處距離底板不同位置臨底流速變化曲線圖,可以發現在這兩個斷面處近底側流速梯度均較大,最大值出現在工況3的斷面2處,達到了52.3 m/(s·m),具體值見表2。

表1 3種模擬工況表

表2 不同斷面最大臨底流速和流速梯度

圖6 1-1斷面臨底流速曲線圖圖7 2-2斷面臨底流速曲線圖
3種工況下溢洪道中軸線處近底板處的湍流動能如圖8所示,可見在1-1斷面和2-2斷面處湍流動能值并不是最大值,其湍流動能值相對于最大值處在一個較小的范圍內。
結合原型觀測情況,發現在1-1斷面和2-2斷面處前后均發現了底板以及墻體的沖刷破壞,如圖9、10所示。對比發現與數值模擬發現的臨底流速變化規律相吻合。即原型底板發生破壞處臨底流速、流速梯度均較大。

圖8 中軸線處湍流動能變化曲線

圖9工程反弧段破壞照片圖10 2-2斷面附近側壁破壞照片
由于臨底流速的增大,在近底處邊界層變薄,流速梯度增大,根據于洪銀等[16]的研究,過大的流速梯度是工程沖刷破壞的主要原因之一,可見該工程出現沖刷破壞的主要原因之一是局部臨底流速及其梯度過大。
(1)通過工程原型觀測數據與數值仿真結果對比表明,采用RNGk-ε和VOF模型仿真模擬溢洪道泄流是可行的。在此基礎上開展泄流臨底流速分布規律研究,可有效解決原型或模型臨底流速難以測量的問題。
(2)工程實際沖刷破壞處與臨底流速數值仿真中最大臨底流速及近底側臨底流速最大梯度出現位置是吻合的,可見引起該工程沖刷破壞的主要原因之一是局部臨底流速及其梯度過大。