三維寬頂堰紊流場(chǎng)數(shù)值模擬與分析

李春成

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

利用三維k-ε紊流模型和VOF方法追蹤自由表面，對(duì)寬頂堰流進(jìn)行數(shù)值模擬，對(duì)寬頂堰流水流特性進(jìn)行了研究，得到了寬頂堰流流場(chǎng)、中心剖面處水面線、流量系數(shù)等結(jié)果，其中比較數(shù)值模擬的中心剖面處水面線是前面學(xué)者沒(méi)有考慮過(guò)的因素，通過(guò)數(shù)值模擬結(jié)果和二維情況及試驗(yàn)實(shí)測(cè)數(shù)據(jù)相比較，說(shuō)明數(shù)值模擬得出的結(jié)果是合理的。

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

1.1 控制方程

1.1.1 質(zhì)量守恒方程（連續(xù)性方程）

任何流動(dòng)問(wèn)題都必須滿足質(zhì)量守恒方程。該定律可表述為：?jiǎn)挝粫r(shí)間內(nèi)流體微元體中質(zhì)量的增加，等于同一時(shí)間間隔內(nèi)流入該微元體的凈質(zhì)量。

三維直角坐標(biāo)系下的連續(xù)性方程為：對(duì)于不可壓縮流體，其流體密度為常數(shù)，連續(xù)性方程簡(jiǎn)化為：

1.1.2 動(dòng)量守恒方程

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

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

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

一般來(lái)講Sx，Sy和Sz是小量，對(duì)于粘性為常數(shù)的不可壓流體，Sx=Sy=Sz=0。

1.1.3 時(shí)均方程

為了考察脈動(dòng)的影響，目前廣泛采用的方法是時(shí)間平均法，即把紊流運(yùn)動(dòng)看作由2個(gè)流動(dòng)疊加而成，一是時(shí)間平均流動(dòng)，二是瞬時(shí)脈動(dòng)流動(dòng)。

時(shí)均形式的N～S方程為：

（注：為方便起見(jiàn)，除脈動(dòng)值的時(shí)均值外，上式中去掉了表示時(shí)均值的上劃線符號(hào)“—”。）

1.1.4 標(biāo)準(zhǔn)k－ε紊流模型

模型中采用三維k－ε紊流數(shù)學(xué)模型。標(biāo)準(zhǔn)kε模型的輸運(yùn)方程為：

紊動(dòng)粘度μt可表示成：

式中：Cμ為經(jīng)驗(yàn)常數(shù)。

Ck是由于平均速度梯度引起的紊動(dòng)能k的產(chǎn)生項(xiàng)，由下式計(jì)算：

Gb是由于浮力引起的紊動(dòng)能k的產(chǎn)生項(xiàng)，對(duì)于不可壓流體，Gb=0。對(duì)于可壓流體，有：

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

YM代表可壓紊流中脈動(dòng)擴(kuò)張的貢獻(xiàn)，對(duì)于不可壓流體，YM=0。對(duì)于可壓流體，有：

在標(biāo)準(zhǔn)k－ε模型中，根據(jù)Launder等的推薦值及后來(lái)的試驗(yàn)驗(yàn)證，模型常數(shù) C1ε=1.44，C2ε=1.92，Cμ=0.99，σk=1.0，σε=1.3。

1.2 自由表面的處理

自由表面的跟蹤采用流體體積分?jǐn)?shù)法，即VOF法，該方法計(jì)算模型適用于兩種或多種互不穿透流體間界面的跟蹤計(jì)算。模型對(duì)每一相引入體積分?jǐn)?shù)變量，通過(guò)求解每一控制單元內(nèi)體積分?jǐn)?shù)確定相間界面，所以，對(duì)于水氣兩相流場(chǎng)，如果αw表示水的體積分?jǐn)?shù),則氣的體積分?jǐn)?shù)可表示αa為：

在一個(gè)單元里，當(dāng)αw=0時(shí)，控制單元內(nèi)沒(méi)有水；αw=1時(shí)，控制單元內(nèi)充滿著水；0<αw<1時(shí)，控制單元內(nèi)包含水氣界面。

水氣界面的跟蹤通過(guò)求解下面的連續(xù)方程來(lái)完成：

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

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

式中：αw為水的體積分?jǐn)?shù)，ρw和ρa(bǔ)分別為水和氣的密度，μw和μa分別為水和氣的分子粘性系數(shù)。通過(guò)水的體積分?jǐn)?shù)αw的求解，ρ和μ的值可由式（16）、（17）求出。

2 邊界條件與求解方法

2.1 邊界條件和初始條件

2.1.1 進(jìn)出口邊界

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

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

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

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

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

對(duì)于紊流參數(shù)，如紊動(dòng)能k和紊動(dòng)耗散率ε的值采用一定的經(jīng)驗(yàn)公式：

式中：Um為進(jìn)口流速，H0為進(jìn)口水深，m。

2.1.2 壁面邊界條件

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

2.1.3 初始條件

在初始流場(chǎng)中，首先對(duì)水槽中水入口以下部分的水的體積分?jǐn)?shù)賦值為1（一直延伸到出口），即從水入口一直到出口處下面部分充滿了水。除此之外，計(jì)算區(qū)域內(nèi)水的體積分?jǐn)?shù)都賦值為0，即除了水的部分其余都充滿了空氣。對(duì)于速度場(chǎng)，整個(gè)計(jì)算區(qū)域的初始速度都賦值為0。

2.2 數(shù)值求解方法

采用交錯(cuò)網(wǎng)格與控制體積法相結(jié)合來(lái)離散計(jì)算區(qū)域，壓力—速度的偶合采用PISO算法（壓力的隱式算子分割算法）進(jìn)行計(jì)算。為了避免出現(xiàn)波動(dòng)的壓力場(chǎng)和速度場(chǎng)，壓力從控制體中心到面上的插值采用Body Force Weighted格式，動(dòng)量方程離散采用一階上風(fēng)格式。

3 數(shù)值模擬結(jié)果分析

3.1 三維流場(chǎng)的速度分布

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

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

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

圖3為下游水位繼續(xù)抬高，下游水位差不多和上游相當(dāng)了，水流入口處的流速很均勻，為均布入口條件，但是水流在經(jīng)過(guò)寬頂堰堰頂時(shí)，斷面收縮，速度變大，在寬頂堰上方的水流狀態(tài)為緩流，寬頂堰堰后下方同樣有順時(shí)針的漩渦，而它的長(zhǎng)度較前面兩種情況的漩渦的長(zhǎng)度要長(zhǎng)。

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

從t=0時(shí)開(kāi)始通過(guò)非穩(wěn)態(tài)的數(shù)值計(jì)算，達(dá)到進(jìn)口和出口流量平衡時(shí)結(jié)束，對(duì)流量Q=70.144 m3/h時(shí)計(jì)算的二維、三維實(shí)測(cè)的自由水面線位置進(jìn)行比較，從比較結(jié)果可以得出，自由出流時(shí)，二維計(jì)算和三維計(jì)算的結(jié)果，與試驗(yàn)測(cè)得的水面線非常接近；接近淹沒(méi)出流時(shí)，水流經(jīng)過(guò)寬頂堰后出現(xiàn)了波狀水波，由于試驗(yàn)時(shí)較難測(cè)出此處的水面線，所以與計(jì)算的有點(diǎn)差別，但是由計(jì)算所得的自由水面線與實(shí)際的相差不大；淹沒(méi)出流時(shí)，三維計(jì)算的自由水面線在最上方，二維計(jì)算的自由水面線在中間，而試驗(yàn)得到的自由水面線在下方，不過(guò)3種在數(shù)值上相差不大。同時(shí)可以得出，二維計(jì)算和三維計(jì)算得到的自由水面線幾乎是吻合的，這是因?yàn)檎麄€(gè)試驗(yàn)水槽是規(guī)則的矩形水槽的緣故，二維計(jì)算和三維計(jì)算相差不大。

3.3 流量系數(shù)比較

表1為二維、三維情況下數(shù)值計(jì)算得到的流量系數(shù)與實(shí)測(cè)流量系數(shù)結(jié)果的比較，表中相對(duì)誤差指相對(duì)于實(shí)測(cè)值計(jì)算得到的誤差。

圖4為流量系數(shù)數(shù)據(jù)點(diǎn)圖，從中可以看到數(shù)值計(jì)算得到的流量系數(shù)與實(shí)測(cè)流量系數(shù)值基本上是相近的。

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

表1 流量系數(shù)的比較

4 結(jié)語(yǔ)

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

2）重點(diǎn)介紹了三維數(shù)值模擬流場(chǎng)和中心剖面處自由水面線位置比較結(jié)果。結(jié)果表明，利用三維k～ε紊流模型，結(jié)合VOF方法，可以較好地計(jì)算帶有復(fù)雜自由表面的泄水建筑物紊流場(chǎng)。同時(shí)對(duì)二維和三維數(shù)值計(jì)算結(jié)果也作了比較，結(jié)果表明，在規(guī)則對(duì)稱(chēng)的水工建筑物上，利用二維計(jì)算來(lái)代替三維計(jì)算是可以的，兩者相差很小。

3）比較了數(shù)值計(jì)算得到的流量系數(shù)和試驗(yàn)實(shí)測(cè)的流量系數(shù)，兩者基本一致，二維計(jì)算流量系數(shù)的最大相對(duì)誤差等于3%，三維計(jì)算流量系數(shù)的最大相對(duì)誤差等于4%。

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

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

［2］陳群，戴光清，劉浩吾.帶有曲線自由表面的階梯溢流壩面流場(chǎng)的數(shù)值模擬［J］.水利學(xué)報(bào)，2002（9）：20-26.

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

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

［5］戴會(huì)超,王玲玲.淹沒(méi)水躍的數(shù)值模擬［J］.水科學(xué)進(jìn)展，2004，15（2）：184-188.

［6］陳娓，陳大宏.溢流堰過(guò)堰流動(dòng)的數(shù)值計(jì)算［J］.人民長(zhǎng)江，2005（1）：40-42.