三維寬頂堰紊流場數值模擬與分析

李春成

（吉林省通化縣英額布水庫工程管理處，吉林 通化 134118）

利用三維k-ε紊流模型和VOF方法追蹤自由表面，對寬頂堰流進行數值模擬，對寬頂堰流水流特性進行了研究，得到了寬頂堰流流場、中心剖面處水面線、流量系數等結果，其中比較數值模擬的中心剖面處水面線是前面學者沒有考慮過的因素，通過數值模擬結果和二維情況及試驗實測數據相比較，說明數值模擬得出的結果是合理的。

1 數學模型

1.1 控制方程

1.1.1 質量守恒方程（連續性方程）

任何流動問題都必須滿足質量守恒方程。該定律可表述為：單位時間內流體微元體中質量的增加，等于同一時間間隔內流入該微元體的凈質量。

三維直角坐標系下的連續性方程為：對于不可壓縮流體，其流體密度為常數，連續性方程簡化為：

1.1.2 動量守恒方程

動量守恒方程也是任何流動系統都必須滿足的基本定律。該定律可表述為：微元體中流體的動量對時間的變化率等于外界作用在該微元體上的各種力之和。該定律實際上是牛頓第二定律的數學表達形式，即N～S方程為：

式中：ρ是流體微元體上的壓力；μ是動力粘度Su，Sv，Sw和是動量守恒方程的廣義源項，Su=Fx+Sx，Sv=Fy+Sy，Sw=Fz+Sz，其中 Sx，Sy和 Sz的表達式如下：

式中：λ是第二粘度，一般可取λ=-2/3。

一般來講Sx，Sy和Sz是小量，對于粘性為常數的不可壓流體，Sx=Sy=Sz=0。

1.1.3 時均方程

為了考察脈動的影響，目前廣泛采用的方法是時間平均法，即把紊流運動看作由2個流動疊加而成，一是時間平均流動，二是瞬時脈動流動。

時均形式的N～S方程為：

（注：為方便起見，除脈動值的時均值外，上式中去掉了表示時均值的上劃線符號“—”。）

1.1.4 標準k－ε紊流模型

模型中采用三維k－ε紊流數學模型。標準kε模型的輸運方程為：

紊動粘度μt可表示成：

式中：Cμ為經驗常數。

Ck是由于平均速度梯度引起的紊動能k的產生項，由下式計算：

Gb是由于浮力引起的紊動能k的產生項，對于不可壓流體，Gb=0。對于可壓流體，有：

式中：Prt是紊動Prandtl數，在該模型中可取Prt=0.85，gi是重力加速度在第i方向的分量，β是熱膨脹系數，可由可壓流體的狀態方程求出，其定義為：

YM代表可壓紊流中脈動擴張的貢獻，對于不可壓流體，YM=0。對于可壓流體，有：

在標準k－ε模型中，根據Launder等的推薦值及后來的試驗驗證，模型常數 C1ε=1.44，C2ε=1.92，Cμ=0.99，σk=1.0，σε=1.3。

1.2 自由表面的處理

自由表面的跟蹤采用流體體積分數法，即VOF法，該方法計算模型適用于兩種或多種互不穿透流體間界面的跟蹤計算。模型對每一相引入體積分數變量，通過求解每一控制單元內體積分數確定相間界面，所以，對于水氣兩相流場，如果αw表示水的體積分數,則氣的體積分數可表示αa為：

在一個單元里，當αw=0時，控制單元內沒有水；αw=1時，控制單元內充滿著水；0<αw<1時，控制單元內包含水氣界面。

水氣界面的跟蹤通過求解下面的連續方程來完成：

根據αw的值就可以得知自由表面的位置。

在VOF模型中，由于水和氣體具有相同的速度場和壓力場，水氣兩相流場可以象單相流場那樣用一組方程來描述。因此引入VOF模型的k－ε紊流模型與單相流的k－ε模型形式完全相同，只是密度ρ和分子粘性系數μ的具體表達式不同，它們是由體積分數的加權平均值給出，即ρ和μ是體積分數的函數，而不是常數。它們可由下式表示：

式中：αw為水的體積分數，ρw和ρa分別為水和氣的密度，μw和μa分別為水和氣的分子粘性系數。通過水的體積分數αw的求解，ρ和μ的值可由式（16）、（17）求出。

2 邊界條件與求解方法

2.1 邊界條件和初始條件

2.1.1 進出口邊界

在實驗室做寬頂堰物理模型試驗，試驗在矩形大型自循環多功能電控變坡試驗水槽中進行，水槽總長12 m，槽凈寬0.3 m，槽總高1.8 m，堰高為0.15 m，長為0.6 m，進口為圓角，1/4圓的半徑為0.04 m，本次計算區域長度為4.5 m，堰前長度為1.0 m。

入口邊界分別由上部的氣體入口和下面的水入口兩部分組成。水的入口采用速度入口邊界條件，為均布分布，流量一定，給定初始水深，計算得到初始速度，通過迭代計算改變水深和速度。

所有的氣體邊界都設為壓力邊界條件，此邊界條件適用于邊界壓力大小已知，但邊界上的通量未知的情況。氣體邊界處的壓力都為大氣壓，但氣體流入或流出邊界的流量未知，所以都采用壓力邊界條件。

在數值模擬寬頂堰水流時，自由出流情況下出口邊界設為壓力出口。由于出口水流為自由出流，與大氣相通，則認為出口壓力為大氣壓力。又由于出口水深未知，因此水和氣體的邊界分不開，只能作為同一個出口邊界，則采用壓力邊界比較合適，一方面氣體可以任意流動，另一方面水也可以自由出流。

在淹沒出流情況下，出口邊界條件則由上部的氣體出口和下面的水出口兩部分組成。氣體出口為壓力出口邊界條件，水的出口也是壓力邊界條件，但壓力值與水深有關，利用用戶自定義函數給出。

對于紊流參數，如紊動能k和紊動耗散率ε的值采用一定的經驗公式：

式中：Um為進口流速，H0為進口水深，m。

2.1.2 壁面邊界條件

壁面采用無滑移邊界，即u=0，v=0，近壁區，由于雷諾數較小，在粘性低層用壁函數法處理。

2.1.3 初始條件

在初始流場中，首先對水槽中水入口以下部分的水的體積分數賦值為1（一直延伸到出口），即從水入口一直到出口處下面部分充滿了水。除此之外，計算區域內水的體積分數都賦值為0，即除了水的部分其余都充滿了空氣。對于速度場，整個計算區域的初始速度都賦值為0。

2.2 數值求解方法

采用交錯網格與控制體積法相結合來離散計算區域，壓力—速度的偶合采用PISO算法（壓力的隱式算子分割算法）進行計算。為了避免出現波動的壓力場和速度場，壓力從控制體中心到面上的插值采用Body Force Weighted格式，動量方程離散采用一階上風格式。

3 數值模擬結果分析

3.1 三維流場的速度分布

利用三維k－ε紊流模型，VOF方法確定自由表面，數值計算了各種流量下的水流情況，得到了流場速度分布。圖1～3為流量Q=70.144 m3/h、不同下游水位條件下、不同位置縱剖面上的速度分布圖，圖中z表示寬度方向的坐標，b為水槽寬度。

圖1為自由出流時的流場速度圖，從兩個縱剖面圖來看基本上是差不多，這是由于試驗裝置是規則的矩形玻璃水槽，邊界對兩剖面的影響不大。從圖1可以看出，水流入口處的流速很均勻，為均布入口條件，但是水流在經過寬頂堰堰頂時速度增加，斷面收縮，水流狀態為急流，流過寬頂堰后，水流再次跌落，速度增加，并且在寬頂堰后面形成了順時針的漩渦。

圖2是下游水位提高時的情況，可以看到，水流入口處的流速很均勻，為均布入口條件，但是水流經過寬頂堰堰頂時，斷面收縮，速度變大，流過寬頂堰后在水面上出現了波狀水流，同時可以看到在堰后下方有順時針的漩渦，并且較前面一種情況產生的漩渦長度要長。

圖3為下游水位繼續抬高，下游水位差不多和上游相當了，水流入口處的流速很均勻，為均布入口條件，但是水流在經過寬頂堰堰頂時，斷面收縮，速度變大，在寬頂堰上方的水流狀態為緩流，寬頂堰堰后下方同樣有順時針的漩渦，而它的長度較前面兩種情況的漩渦的長度要長。

3.2 中心剖面處自由水面線位置比較

從t=0時開始通過非穩態的數值計算，達到進口和出口流量平衡時結束，對流量Q=70.144 m3/h時計算的二維、三維實測的自由水面線位置進行比較，從比較結果可以得出，自由出流時，二維計算和三維計算的結果，與試驗測得的水面線非常接近；接近淹沒出流時，水流經過寬頂堰后出現了波狀水波，由于試驗時較難測出此處的水面線，所以與計算的有點差別，但是由計算所得的自由水面線與實際的相差不大；淹沒出流時，三維計算的自由水面線在最上方，二維計算的自由水面線在中間，而試驗得到的自由水面線在下方，不過3種在數值上相差不大。同時可以得出，二維計算和三維計算得到的自由水面線幾乎是吻合的，這是因為整個試驗水槽是規則的矩形水槽的緣故，二維計算和三維計算相差不大。

3.3 流量系數比較

表1為二維、三維情況下數值計算得到的流量系數與實測流量系數結果的比較，表中相對誤差指相對于實測值計算得到的誤差。

圖4為流量系數數據點圖，從中可以看到數值計算得到的流量系數與實測流量系數值基本上是相近的。

同時，由于本試驗是在規則的玻璃水槽中進行的，并且寬頂堰是規則的圓角進口的寬頂堰，所以二維和三維計算所得到的流量系數大致相同，這個從上面的數據和圖都可以看出。

表1 流量系數的比較

4 結語

1）本文利用三維k～ε紊流模型和VOF自由面調整方法，較好地模擬了帶有曲線自由表面的寬頂堰的流場，數值模擬的結果是合理的。

2）重點介紹了三維數值模擬流場和中心剖面處自由水面線位置比較結果。結果表明，利用三維k～ε紊流模型，結合VOF方法，可以較好地計算帶有復雜自由表面的泄水建筑物紊流場。同時對二維和三維數值計算結果也作了比較，結果表明，在規則對稱的水工建筑物上，利用二維計算來代替三維計算是可以的，兩者相差很小。

3）比較了數值計算得到的流量系數和試驗實測的流量系數，兩者基本一致，二維計算流量系數的最大相對誤差等于3%，三維計算流量系數的最大相對誤差等于4%。

4）數值模擬圓角進口寬頂堰的水流情況，這是參照實驗室里面的物理模型，所以它缺少了邊墩、岸墻及其他水工建筑物的影響，在實際工程問題計算時，則需要把這些影響因素考慮進去，這樣才能更好地模擬真實流動，才能得到更精確的研究結果。

［1］袁新明，毛根海，吳壽榮，唐錦春.閘后水躍水流的數值模擬［J］.水利學報，1999（5）：44-48.

［2］陳群，戴光清，劉浩吾.帶有曲線自由表面的階梯溢流壩面流場的數值模擬［J］.水利學報，2002（9）：20-26.

［3］李志勤，李洪，李嘉，李然.溢流丁壩附近自由表面的試驗研究與數值模擬［J］.水利學報，2003（8）：53-57.

［4］刁明軍，楊永全，王玉蓉，劉善均.挑流消能水氣二相流數值模擬［J］.水利學報，2003（9）：77-82.

［5］戴會超,王玲玲.淹沒水躍的數值模擬［J］.水科學進展，2004，15（2）：184-188.

［6］陳娓，陳大宏.溢流堰過堰流動的數值計算［J］.人民長江，2005（1）：40-42.