倪誠陽，唐 戀，呂文龍，王 航

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065；2.長江三峽技術(shù)經(jīng)濟發(fā)展有限公司，重慶 401120)

1 概述

水躍在上游急流和下游緩流之間存在著流態(tài)復(fù)雜的旋滾區(qū)，通常伴隨著大尺度渦流的產(chǎn)生、自由面波動、能量耗散以及大量摻氣[1]。河床坡度和粗糙度、過水?dāng)嗝娴耐蝗蛔兓蛘呷斯そㄖ锏雀鞣N形式的流動阻力常常在天然河道和水利工程中引發(fā)水躍，為了流體混摻或污水處理，水躍也同樣被廣泛應(yīng)用在工業(yè)加工中。

國內(nèi)外研究人員針對水躍的流動形態(tài)和水氣兩相流特性開展了大量研究。Hager[2]通過實驗觀測給出了水躍長度和自由面輪廓的經(jīng)驗描述；Chanson和Brattberg[3]使用電阻式兩相流探直接測量水躍的水氣兩相流特性，提出摻氣濃度和水氣界面速度的典型垂直分布。劉沛清[4]論述了水躍擴散段內(nèi)流速分布的理論分析結(jié)果，其可靠性通過與實驗資料的對比得到了證實；張沁等[5]對多種工況下不同形式的水躍進行實驗研究，測量并分析了水躍消能率，建議在實際消能工程中多采用穩(wěn)定水躍。數(shù)值模型可以方便地改變水流和邊界條件，并以比物理模型低得多的成本詳細顯示內(nèi)部流動結(jié)構(gòu)[6]，因此被廣泛用于水工結(jié)構(gòu)設(shè)計和明渠流動研究。隨著計算能力和水氣兩相流基礎(chǔ)理論的發(fā)展，關(guān)于水躍的數(shù)值研究越來越多地出現(xiàn)。最常用的方法是求解多相形式的雷諾平均Navier-Stokes(RANS)方程、湍流模型以及界面追蹤模型。Bayon等[7]對比模擬與實驗的水躍自由面輪廓、時均水平速度，評估了兩個RANS平臺建模水躍的能力?；菘礫8]等采用FLUENT軟件分析不同坡度對水躍躍長的影響；鄭鐵剛等[9]通過RANS方法模擬了來流弗勞德數(shù)為3.19的淹沒水躍，分析內(nèi)部流場的同時發(fā)現(xiàn)水躍旋滾顯著影響著紊動的產(chǎn)生和耗散，但對壓力分布影響較??；程香菊等[10]使用同樣的方法模擬波浪形底板上的水躍，展示了水躍的演變過程、躍長以及流速分布。

摻氣濃度和水躍的部分自相似性在水環(huán)境領(lǐng)域和工程應(yīng)用中至關(guān)重要，而對此進行的數(shù)值研究卻很少見。本文基于二維混合模型NEWFLUME[11]對3個來流弗勞德數(shù)下的經(jīng)典水躍進行數(shù)值研究，對比了模擬結(jié)果與Wang等[12]的實驗數(shù)據(jù)集，包括自由面形態(tài)、時均水平速度和摻氣濃度分布、以及斷面最大流速和摻氣濃度的沿程自相似分布，系統(tǒng)地驗證了此模型的可靠性。

2 數(shù)學(xué)控制方程

水利工程中水相和氣相的流動速度遠小于聲速，混摻的氣泡和水被視為不可壓縮的混合流體，平均運動由雷諾平均Navier-Stokes方程描述，混合模型的連續(xù)性方程和動量方程如下：

?um=0

(1)

(2)

式中，um、ρm、Pm—兩相流體的混合速度、密度和壓力；g—重力加速度；t—時間；τGm—混合流體的黏性應(yīng)力張量和湍流應(yīng)力張量；τDm—氣泡和液體相對運動引起的擴散應(yīng)力張量。

應(yīng)力張量的表達公式如下：

(3)

τDm=-ρmC(1-C)urur

(4)

式中，I—單位張量；k—混合流體的紊動能；C—摻氣濃度，定義為每單位體積氣-水混合物的氣相占比；ur—氣泡和水之間的相對速度；νm、νTm—混合流體的運動黏度和運動渦流黏度。

根據(jù)每個計算網(wǎng)格的摻氣濃度確定混合物的密度和運動黏度，得出：

ρm=ρgC+ρl(1-C)，νm=νgC+νl(1-C)

(5)

式中，ρg、ρl—空氣和水的密度；νg、νl—空氣和水的運動黏度。

對于氣泡在靜水中的上浮速度urs，Comolet[13]強調(diào)表面張力和浮力的重要性并建議：

(6)

式中，d—氣泡的平均直徑；σt—空氣-水表面張力。

在摻氣濃度為C的水氣兩相流中，氣泡上浮速度ur的表達式變?yōu)閇14]：

(7)

采用標(biāo)準(zhǔn)k-ε湍流模型來模擬水躍的湍流特性，該湍流模型計算成本較低且不易出現(xiàn)收斂問題：

(8)

(9)

(10)

式中，ε—湍動能耗散率；Ps—紊動能生成項。

Ps表達式及方程中的系數(shù)設(shè)置如下：

Ps=τGm?um

(11)

Cd=0.09，σk=1.0，σε=1.0，C1ε=1.44，C2ε=1.92

(12)

VOF(流體體積)模型的有效性在廣泛的流動建模中得到了證明，在模擬劇烈水氣混合現(xiàn)象時，既需要計算水體中摻氣濃度的分布，又需要追蹤紊動自由面，相較于單相流和存在清晰界面的兩相流復(fù)雜得多。本研究在傳統(tǒng)VOF方法的基礎(chǔ)上引入額外的對流項以控制數(shù)值耗散并提供分辨率更高的相界面，通過求解該優(yōu)化后的模型提高運動自由面和摻氣濃度的模擬精度[15]：

(13)

通過求解方程(13)獲得計算域中每個網(wǎng)格單元的摻氣濃度，然后根據(jù)摻氣濃度分布使用二階精度的Youngs算法重建自由面。Youngs算法允許每個網(wǎng)格中重建的兩相界面為斜線，而不局限為垂線或者水平線，可以有效減少數(shù)值耗散、數(shù)值色散性和不穩(wěn)定性，從而更準(zhǔn)確地獲得貼合實驗現(xiàn)象的清晰自由面[16]。

3 物理模型及算例設(shè)置

Wang等[12]發(fā)表于2015年的論文描述了一組經(jīng)典水躍實驗。橫斷面為矩形的水平平底水槽，其長3.2m，寬0.5m，深0.41m，由光滑的高密度聚乙烯底板和玻璃側(cè)壁構(gòu)成。上游高位水箱具有恒定水頭，通過半圓形水閘形成水平入流。對于給定的流量，調(diào)節(jié)渠道末端尾水閘門的高度改變下游水位，從而控制水躍的發(fā)生位置。實驗裝置示意圖及相關(guān)符號如圖1所示，圖1中x是自閘門進口處沿入流方向的水平坐標(biāo)，y是自渠道底板向上的垂直坐標(biāo)。固定上游閘門開度h=0.03m，躍趾位于縱向位置Xt=1.25m。水躍的來流邊界層厚度δ

表1 實驗流動條件總結(jié)

圖1 經(jīng)典水躍模型實驗概化圖

如圖2所示，實驗使用多個超聲測距儀(ADMs)測量水躍自由面高程，并使用雙針的電阻式兩相流探針(PDPs)測量水氣兩相流特性。超聲測距儀朝向檢測水面，發(fā)射聲束并接收水面反射的聲束，由聲束行進時間獲得傳感頭與被檢測水面之間的距離。固定于中心線上的探針有兩個平行但長度不同的針式傳感器，正對來流方向。探針信號是與內(nèi)部電極導(dǎo)電率成比例的電壓樣本，根據(jù)電壓變化來區(qū)分氣相和水相。采用氣相和水相電壓水平之和的50%為臨界值對信號進行二值化處理，其中0代表水，1代表空氣，從而獲得時均摻氣濃度。兩個傳感器的原始信號之間的互相關(guān)分析進一步提供氣泡-水界面的水平速度。

圖2 水躍實驗中的儀器

二維計算域與水槽長度相同，高度大于躍后水深。在計算域左邊界的閘門出口處設(shè)定與實驗條件相同的均勻水平速度，對固邊界施加無滑移條件，在計算域右邊界的過沖閘門處施加自由出流條件，對自由面邊界施加零剪應(yīng)力條件。邊界層內(nèi)的湍流發(fā)展未被解析，固邊界附近的湍流場用壁面函數(shù)處理。

計算域被離散為統(tǒng)一尺寸的網(wǎng)格系統(tǒng)，采用4種不同尺寸的網(wǎng)格進行收斂性分析。網(wǎng)格的長和高(Δx×Δy)分別為0.025m×0.003m、0.03m×0.004m、0.03m×0.005m、0.035m×0.005m，最大與最小網(wǎng)格尺寸的比值大于1.3，滿足Celik等[17]的建議。選取Fr1=5.1水躍的摻氣濃度分布作為收斂指標(biāo)，圖3比較了相對位置(x-Xt)/d1=11.6處的模擬與實驗結(jié)果。同時斷面平均摻氣濃度Cmean的相對誤差在表2中顯示，為滿足相對誤差小于10%的要求[18]，采用0.025m×0.003m的網(wǎng)格尺寸。本文根據(jù)實驗經(jīng)驗采用45s的模擬時長，以確保特征變量的收斂。

表2 模擬的斷面平均摻氣濃度與實驗結(jié)果的相對誤差

圖3 使用不同尺寸網(wǎng)格模擬的摻氣濃度分布與實驗數(shù)據(jù)的比較

4 模擬結(jié)果與分析

4.1 水躍形態(tài)

經(jīng)典水躍共軛水深比d2/d1的理論關(guān)系式：

(14)

如圖4所示，模擬得到的弗勞德數(shù)為3.8、5.1和7.5水躍的共軛水深比分別為4.95、6.62和9.97，此結(jié)果與公式(14)的比較顯示出良好的一致性。模擬的無量綱水躍長度Lj/d1分別為18.7、27.0和45.2，共軛水深比和水躍長度的模擬結(jié)果與實驗值的相對誤差均小于3%。水躍消能率與共軛水深比密切相關(guān)，所以能夠預(yù)見本模型對經(jīng)典水躍消能率的準(zhǔn)確預(yù)測。

在水躍長度范圍內(nèi)(0<(x-Xt)/Lj<1)，實驗測量的自由面形態(tài)表現(xiàn)出自相似性[19]：

(15)

式中，η—水躍長度內(nèi)的時均自由面高程，接近特征高程Y50(摻氣濃度C=0.5對應(yīng)的高程)。

兩相流探針還記錄了定義摻氣水流上邊界的Y90(C=0.9對應(yīng)的高程)，同樣具有自相似輪廓[19]：

(16)

從計算域的時均摻氣濃度分布中提取模擬的特征高程Y90和Y50。圖5將模擬結(jié)果與擬合公式進行比較，并計算出Y50和Y90的決定系數(shù)分別為R2=0.980、R2=0.984，表明此數(shù)值模型能夠在模擬時長內(nèi)準(zhǔn)確再現(xiàn)水躍的流動形態(tài)。

圖5 模擬自由面輪廓和擬合曲線的比較

4.2 時均水平速度

水躍速度場包括渠底上方形成的邊界層、湍流剪切層中的正向速度以及自由表面回流區(qū)中的負(fù)速度。圖6對比了3個工況下的時均水平速度的模擬結(jié)果和實驗值，每個工況包含4個不同位置的垂直截面。注意混合模型的模擬速度是混合流體速度，不區(qū)分氣相和水相，而兩相流探針測量的實驗數(shù)據(jù)反映氣泡-水界面的速度，更近似于氣泡速度。模擬結(jié)果和測量值在高度較低(相對高程0

圖6 不同截面上水平速度的模擬與實驗結(jié)果比較

圖7展示了模擬的無量綱最大速度Umax/U1的流向衰減。先前的實驗研究表明，在不同流動條件下，Umax/U1在水躍長度范圍內(nèi)遵循相似的衰減曲線，可被方程(17)所擬合[20]，模擬結(jié)果與此非常一致。

圖7 模擬無量綱最大速度的流向衰減以及與經(jīng)驗公式的比較

(17)

4.3 時均摻氣濃度

摻氣過程是水躍最顯著的特征之一，且摻氣濃度對于滿足生化需氧量、防止空化空蝕破壞、確定水工建筑物的安全超高至關(guān)重要。雖然此模型無法模擬單個氣泡的形成及動力學(xué)過程，但通過追蹤宏觀水氣界面以及截留空氣的對流擴散，能夠很好地預(yù)測水躍的摻氣濃度分布。

圖8比較了時均摻氣濃度分布的模擬與實驗結(jié)果，在所有流動條件下，摻氣濃度及其特征值的分布具有相似性。局部最小摻氣濃度將分布曲線清晰地分為兩部分，在其下方和上方分別是摻氣濃度呈鐘形分布的射流剪切區(qū)和單調(diào)增加的回流區(qū)。圖8的對比反映模型低估了渠底附近的摻氣濃度，模擬的鐘形分布較窄，與實驗相比，摻氣濃度在垂向上的擴散受到了限制。目前的混合模型不能模擬出渦旋結(jié)構(gòu)攜帶的分散氣泡對渠底貢獻的摻氣濃度，因此需要進一步研究空氣擴散模型來改進強摻混過程的建模。模型以較高的精度模擬出摻氣濃度在上方回流區(qū)的增長趨勢，以及由于摻氣而導(dǎo)致的水體膨脹效應(yīng)，雖然實驗記錄的回流區(qū)由于自由面的液滴飛濺而更厚。

圖8 不同截面上模擬和實驗測量的摻氣濃度比較

隨著空氣的垂向擴散，湍流剪切層中的最大摻氣濃度Cmax在流向上減小，其相應(yīng)高程yCmax隨著流動深度的增加而增加。Wang等[20]根據(jù)實驗數(shù)據(jù)擬合出它們的變化曲線：

(18)

(19)

圖9將模擬的Cmax和yCmax與擬合公式進行了比較。0<(x-Xt)/Lj<0.3時，Cmax的模擬值偏大，歸因于模型對垂向擴散系數(shù)的低估?？傮w而言，模型很好地預(yù)測了Cmax的指數(shù)下降和yCmax的線性增加，并且成功地再現(xiàn)了它們在不同弗勞德數(shù)下的自相似性。

圖9 摻氣濃度分布的自相似性以及與經(jīng)驗公式的比較

5 結(jié)語

本文通過NEWFLUME混合模型模擬了3個弗勞德數(shù)的經(jīng)典水躍，憑借模擬結(jié)果與實驗數(shù)據(jù)的全面對比得出主要結(jié)論如下：

(1)使用RANS方法建模水躍時，不同尺寸的計算網(wǎng)格將得到差異顯著的模擬結(jié)果，因此網(wǎng)格敏感性分析必不可少。

(2)混合模型可以準(zhǔn)確再現(xiàn)水躍的時均自由面輪廓、消能率、基本的水氣兩相流特性，以及不同弗勞德數(shù)下速度和摻氣濃度表現(xiàn)出的自相似性。

(3)計算結(jié)果與高質(zhì)量數(shù)據(jù)的交叉比較對于模型驗證和分析各研究方法的可靠性不可或缺。兩相流探針在水躍回流區(qū)對負(fù)速度的測量出現(xiàn)高估，而混合模型會低估湍流剪切層中垂向上的空氣擴散速率，需要進一步改進本模型對強水氣摻混過程的模擬。