經典水躍水氣兩相流特性的數值模擬

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

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

1 概述

水躍在上游急流和下游緩流之間存在著流態復雜的旋滾區，通常伴隨著大尺度渦流的產生、自由面波動、能量耗散以及大量摻氣[1]。河床坡度和粗糙度、過水斷面的突然變化或者人工建筑物等各種形式的流動阻力常常在天然河道和水利工程中引發水躍，為了流體混摻或污水處理，水躍也同樣被廣泛應用在工業加工中。

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

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

2 數學控制方程

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

?um=0

(1)

(2)

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

應力張量的表達公式如下：

(3)

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

(4)

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

根據每個計算網格的摻氣濃度確定混合物的密度和運動黏度，得出：

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

(5)

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

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

(6)

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

在摻氣濃度為C的水氣兩相流中，氣泡上浮速度ur的表達式變為[14]：

(7)

采用標準k-ε湍流模型來模擬水躍的湍流特性，該湍流模型計算成本較低且不易出現收斂問題：

(8)

(9)

(10)

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

Ps表達式及方程中的系數設置如下：

Ps=τGm?um

(11)

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

(12)

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

(13)

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

3 物理模型及算例設置

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

表1 實驗流動條件總結

圖1 經典水躍模型實驗概化圖

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

圖2 水躍實驗中的儀器

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

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

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

圖3 使用不同尺寸網格模擬的摻氣濃度分布與實驗數據的比較

4 模擬結果與分析

4.1 水躍形態

經典水躍共軛水深比d2/d1的理論關系式：

(14)

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

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

(15)

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

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

(16)

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

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

4.2 時均水平速度

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

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

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

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

(17)

4.3 時均摻氣濃度

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

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

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

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

(18)

(19)

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

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

5 結語

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

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

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

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