劉 文，賀昌海，江 維，覃天林，傅少君

（1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2.廣西南寧水利電力設(shè)計院，南寧 530001；3.西京學(xué)院，西安 710123）

數(shù)值模擬方法可以在花費較小的情況下獲得全息水流流場，得到實際工程中不同方案下的流場水力參數(shù)，目前已有較多關(guān)于溢洪道內(nèi)流場的數(shù)值模擬研究。耿敬、馬世領(lǐng)等［1］建立龍頭橋水庫溢洪道三維結(jié)構(gòu)模型，結(jié)合物理模型試驗結(jié)果得到溢洪道的泄流能力、水面線、流速分布、堰面壓強等水力特性，分析了水冠和折沖水流的成因，對中墩和尾墩進行優(yōu)化，在水冠和折沖水流現(xiàn)象消弱方面取得了很好的效果。周招、王均星等［2］通過數(shù)值模擬研究不同方案消力池內(nèi)流態(tài)、紊動能分布、摻氣濃度分布、流速分布，得到非完全寬尾墩方案不僅能增強水流紊動、調(diào)整水流流態(tài)，更能提升消能效果，進而為消力池內(nèi)水流二次水躍、消能不充分等問題提出解決方案。程香菊、陳永燦等［3］基于雙流體連續(xù)介質(zhì)模型，對階梯溢流壩面兩相流進行模擬計算，得到流速、初始摻氣點的位置及負壓分布等特征參數(shù)，與模型試驗結(jié)果相比，其結(jié)果符合實際，為了解非摻氣區(qū)范圍及消除空蝕破壞提供理論依據(jù)。蘇燕，張挺等［4］采用雙方程紊流模型模擬溢洪道內(nèi)有復(fù)雜水氣交界面的流場，得到流場流速、壓強分布及摻氣特性，其中最大沖擊壓強計算誤差較大，進一步通過定性分析比選得到合適的消能方式。覃昕慧［5］采用雙流體模型，對摻氣槽附近的摻氣水流進行模擬，同時在考慮水氣相互作用的前提下，以0.5、1、2 mm 特征直徑氣泡代表氣相，得到通氣量特性、壓力分布、摻氣濃度分布。其中通氣孔處風(fēng)速和通氣量的計算值偏大，氣泡大小及變化規(guī)律還沒有明確結(jié)果。高學(xué)平、賈來飛等［6］基于雙流體歐拉法對糯扎渡開敞式溢洪道摻氣挑坎摻氣水流進行模擬，針對摻氣坎高、坡度、流速及不同單寬流量的多種組合工況，研究了空腔負壓、空腔長度及摻氣濃度分布，驗證了雙流體歐拉模型模擬高速摻氣水流的可行性。

已有研究成果表明數(shù)值模擬結(jié)果與模型試驗結(jié)果吻合較好，能夠滿足實際工程需要。但是，已有的研究成果沒有考慮不同數(shù)學(xué)模型對局部流場的影響，在溢洪道內(nèi)局部流場中的計算成果不能很好地反映實際情況。

本文綜合采用不同數(shù)學(xué)模型（卷氣模型、VOF 模型、歐拉雙流體模型），應(yīng)用不同軟件進行數(shù)值模擬，研究泄水建筑物在摻氣條件下局部流場的變化，同時與模型試驗結(jié)果進行對比分析，以期獲得更加合理的數(shù)值模擬成果，為溢洪道內(nèi)流場研究及相關(guān)設(shè)計提供參考。

1 三維模型

某工程溢洪道底孔溢流面高程485.00 m，孔口尺寸8.0 m（寬）×8.5 m（高）。表孔溢流面高程507.00 m，孔寬15 m。在實際工程二期導(dǎo)流階段，溢洪道表孔溢流堰體不施工，預(yù)留缺口作為導(dǎo)流泄水建筑物，缺口底部高程485.00 m（圖1）。

用CATIA 建立溢洪道三維幾何模型（圖2），導(dǎo)入ICEM 中，通過幾何編輯功能完成流場幾何域建立。局部流場部分計算設(shè)溢洪道長度為x軸，水流流向為負x方向，寬度方向為y軸，水深方向為z軸。

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

2.1 紊流模型

本文研究的溢洪道作為導(dǎo)流泄水建筑物，其地形起伏多變，消力池區(qū)域水流流態(tài)復(fù)雜、流線彎曲程度大，故采用RNGk-ε模型［7］?？刂品匠倘缦拢?/p>

連續(xù)方程：

動量方程：

紊動能k方程：

紊動能耗散率ε方程：

式中：ui為流速分量；Ai、gi、fi分別為三維直角坐標(biāo)方向上可流動的面積分數(shù)、重力加速度和黏滯力；VF是可流動的體積分數(shù)；ρ是流體密度；p是作用在流體微元上的壓力；μt是紊動黏滯系數(shù)；Gk是紊動能k的產(chǎn)生項：σk、σε是紊動能和耗散率對應(yīng)的Prandtl 數(shù)，均為1.39；，經(jīng)驗常數(shù)Cε1、Cε2分別取1.42、1.68。

2.2 卷氣模型

消力池中水流伴隨著摻氣現(xiàn)象，采用Flow-3D 中的卷氣模型（Air Entrainment Model）?？紤]到流體體積膨脹及卷入氣體的浮力效應(yīng)，卷氣模型中還要使用密度變化方程。卷氣模型基于3種因素：由紊流產(chǎn)生的擾動、重力及表面張力。紊流是卷氣過程產(chǎn)生的主要因素。流體表面摻氣是因為紊流渦體在流動過程中將表面空氣卷入流體中，其主要取決于紊流強度是否克服重力和表面張力組成的表面穩(wěn)定力［8］。

紊流渦體的特征尺寸：

單位體積的紊動能：

表面穩(wěn)定力：

單位時間卷氣體積量：

式中：Q為紊動能；D為耗散函數(shù)；在RNG 紊流模型中，cnu的值為0.085，用來描述表面擾動的特征；ρ是液體密度；gn是重力加速度在自由表面法線方向的分量；Lt流體單元升高高度；As是表面張力系數(shù)；As是表面面積；Cair是一個經(jīng)驗參數(shù)，表示單位時間內(nèi)有一部分表面積卷入空氣，初步假設(shè)值為0.5，文獻中已驗證了其在數(shù)值模擬中的合理性。當(dāng)紊動能Pt小于表面穩(wěn)定力Pd時，δV值為0。

2.3 歐拉模型

歐拉模型假設(shè)各相流體間可以空間共存和相互貫穿，并且每相各自滿足質(zhì)量和動量守恒定律，其中體積分數(shù)表示為，代表了每相所占據(jù)的空間。歐拉模型中，時間平均方程、空間平均方程可以通過對瞬時及局部方程求平均得到，同樣的方式也能得到相間作用表達式［9，10］。

q相的體積Vq為：

且q相的有效密度為：

各相之間體積分數(shù)滿足：

歐拉模型的守恒方程包括質(zhì)量守恒方程、動量守恒方程和能量守恒方程，本文涉及水氣二相流，只考慮質(zhì)量守恒和動量守恒。

q相的連續(xù)方程可寫為：

又由q相的動量平衡得動量方程：

式中：是第q相的壓力應(yīng)變張量，其表達式如下：

式中：ρq是q的物理密度；υq是q相的速度，mpq為從p相到q相的質(zhì)量傳遞，由質(zhì)量守恒得到mpq= -mqp，mpp= 0。是外部體積力；是升力；是虛擬質(zhì)量力；μq和λq是q相的剪切和體積黏度；P是所有相共享的壓力是相之間的相互作用力；是相間的速度，其定義如下：如果mpq＞0（也就是相p的質(zhì)量轉(zhuǎn)移到相q），；反之，則有和方程mpp= 0 需要有適當(dāng)?shù)南嚅g表達式將相間作用力封閉，其中相間作用力滿足條件

3 三維數(shù)值模擬

采用卷氣模型進行數(shù)值模擬計算過程及模擬結(jié)果在文獻中［7］已經(jīng)詳細說明，采用VOF 模型和歐拉雙流體模型模擬計算時，由于閘墩前緣和消力坎位置結(jié)構(gòu)復(fù)雜，網(wǎng)格加密至0.3 m 精度，對其他區(qū)域，運用六面體核心網(wǎng)格技術(shù)劃分邊長為2 m左右的立方體網(wǎng)格，最終得到約134 萬單元。消力池部分所畫網(wǎng)格的網(wǎng)格線與水流方向一致，減少了數(shù)值模擬計算對流場的不利影響。

邊界條件設(shè)置如下：采用VOF 模型時，在模型的上游邊界設(shè)置為流量入口，并給定上游水位和底板高程；下游邊界為壓力出口邊界，同時給定下游水位；水氣交界面設(shè)置為壓力入口，其余為壁面邊界。采用歐拉雙流體模型時，入流斷面包括水進口和空氣進口，水進口為質(zhì)量流量入口邊界，氣進口、上表面和出流斷面采用壓力出口。溢洪道底部和側(cè)面采用無滑移壁面邊界。

用VOF 模型計算時，基于有限體積法，采用RNG 紊流模型，隱格式迭代求解，為了加快求解速度，選擇在網(wǎng)格扭曲度較大時有明顯優(yōu)勢的PISO算法。用歐拉雙流體模型計算時，基于有限體積法，采用RNG 紊流模型，隱格式迭代求解，同時，壓力速度耦合求解采用SIMPLE 的擴展算法“Phase Coupled SIMPLE”方法。

計算工況見表1。

表1 數(shù)值模擬工況Tab.1 The cases of numeric simulation

4 數(shù)值模擬成果對比分析

4.1 流 態(tài)

流態(tài)對比結(jié)果如圖4、5 所示，可看出數(shù)值模擬和模型試驗?zāi)M的缺口進口水流都比較平靜。工況2 中，試驗的閘室出口流速較大，流態(tài)紊亂，出口處有明顯的水躍產(chǎn)生，水流摻氣明顯，數(shù)值模擬流態(tài)與之相符合。

4.2 流 速

由圖6、7 中可知，工況一和工況二中最大流速均出現(xiàn)在摻氣坎（樁號0-50.8）附近。在各個測點中，消力池末端處流速在兩個工況中均為最小值，進一步說明水流進入消力池后經(jīng)過摻氣后消耗大量能量，流速在消力池逐漸減小，下泄急流迅速變?yōu)榫徚鳎Τ鼐邆浜芎玫南苄Ч?/p>

一路上，我怏怏不樂，老大不情愿，幾年前那里的情景又一幕幕地浮現(xiàn)在眼前。然而， 這次故地重游，我發(fā)現(xiàn)，這里的土，這里的人，一切都變了，短短幾年時間，這里就煥然一新，再也不是以前的樣子了。

不同工況下，3 種數(shù)學(xué)模型和模型試驗結(jié)果誤差對比如表2。

表2 流速計算結(jié)果對比Tab.2 The contrast of velocity

在工況1 中，模型試驗結(jié)果、卷氣模型、VOF 模型及歐拉模型模擬結(jié)果參見圖6。與模型試驗結(jié)果相比較，卷氣模型計算結(jié)果：最大絕對誤差1.43 m/s（0-99.6），最大相對誤差0.49，絕對誤差平均值0.83 m/s，相對誤差平均值0.28；VOF 模型計算結(jié)果：最大絕對誤差3.00 m/s（0-19.6），最大相對誤差0.56，絕對誤差平均值1.54 m/s，相對誤差平均值0.37；歐拉模型計算結(jié)果最大絕對誤差2.39 m/s（0-19.6），最大相對誤差0.25，絕對誤差平均值0.77 m/s，相對誤差平均值0.12。在此工況下，歐拉模型流速模擬結(jié)果相比其他兩個模型模擬結(jié)果更好一些。

在工況2 中，模型試驗結(jié)果、卷氣模型、VOF 模型及歐拉模型模擬結(jié)果參見圖7。卷氣模型計算結(jié)果最大絕對誤差1.31 m/s（0-99.6），最大相對誤差0.15，絕對誤差平均值0.75 m/s，相對誤差平均值為0.08；VOF 模型計算結(jié)果最大絕對誤差0.77 m/s（0+38.2），最大相對誤差0.10，絕對誤差平均值0.40 m/s，相對誤差平均值為0.04。歐拉模型計算結(jié)果最大絕對誤差1.93 m/s（0+38.2），最大相對誤差0.26，絕對誤差平均值1.18 m/s，相對誤差平均值為0.14。VOF 模型模擬結(jié)果相比其他兩個模型流速模擬結(jié)果更好一些。總體而言，在工況1 和工況2 中采用卷氣模型、VOF 模型及歐拉模型模擬局部流場內(nèi)水流速度大小及分布差別較小，三者均能較好地模擬實際情況。

4.4 壓強分布

用卷氣模型、VOF 模型及歐拉模型模擬得到消力池內(nèi)各測點壓強，對比結(jié)果如圖8、9 所示。上游壓強較小，無劇烈變化，閘室出口和消力池內(nèi)壓力變化較劇烈且消力池底部壓力最大。與模型試驗結(jié)果相比較（表3），工況1 中，卷氣模型計算結(jié)果：最大絕對誤差1.51 m（0+12.9），最大相對誤差0.34，絕對誤差平均值1.16 m，相對誤差平均值0.15；VOF 模型計算結(jié)果：最大絕對誤差2.3 m（0-12.0），最大相對誤差0.72，絕對誤差平均值0.96 m，相對誤差平均值0.16；歐拉模型計算結(jié)果：最大絕對誤差1.51 m，最大相對誤差0.34，絕對誤差平均值1.07 m，相對誤差平均值0.14。三種數(shù)學(xué)模型計算結(jié)果與模型試驗符合較好。

表3 壓強計算結(jié)果對比Tab.3 The contrast of pressure

工況2 中，卷氣模型計算結(jié)果：最大絕對誤差3.6 m（0-12.0），最大相對誤差0.55，絕對誤差平均值1.36 m，相對誤差平均值0.16；VOF 模型計算結(jié)果：最大絕對誤差10.7 m（0-12.0），最大相對誤差1.62，絕對誤差平均值6.0 m，相對誤差平均值0.59。歐拉模型計算結(jié)果：最大絕對誤差10.0 m（0-12.0），最大相對誤差1.52，絕對誤差平均值6.93 m，相對誤差平均值0.65。卷氣模型計算結(jié)果與模型試驗結(jié)果符合較好，VOF 模型和歐拉模型計算結(jié)果誤差較大。

對比VOF 模型和歐拉模型計算結(jié)果，工況1 中兩種模型的最大誤差值-1.4 m 水柱，誤差平均值0.85 m，相對誤差平均值10%。工況2 中兩種模型的最大誤差值3.7 m 水柱，誤差平均值1.4 m，相對誤差平均值9%。同時，通過與模型試驗值比較，考慮水流摻氣時，在消力池外模擬結(jié)果較好，消力池中壓強有較大偏差。

4.5 紊動能分布

紊動能是衡量下泄水流湍流耗散劇烈程度的物理量。觀察水流沿著溢洪道中心線各斷面的紊動能分布，可以發(fā)現(xiàn)以下規(guī)律。工況1 中，卷氣模型計算結(jié)果表明（圖10）：其最大值發(fā)生在在摻氣坎（樁號0-50.8）前部，消力池內(nèi)紊動能在0.6～3.90 m2/s2之間。VOF 模型計算結(jié)果表明：其最大值發(fā)生在在摻氣坎前部、尾坎處（樁號0-160.0），為0.28 m2/s2，且在消力池底部紊動能很小。歐拉雙流體模型計算結(jié)果表明：斷面紊動能最大值發(fā)生在摻氣坎附近、尾坎處，為0.41 m2/s2，在消力池內(nèi)紊動能很小，消力池底部紊動能比消力池中部略大，在0.16～0.24 m2/s2之間?？傮w而言，同模型試驗相比，本工況流速較低，紊動能較小，水流較平靜。

4.6 水流摻氣率

數(shù)值模擬在一定程度上能反映摻氣過程及摻氣分布，在工程實際中，確定初始摻氣點的準(zhǔn)確位置，將有助于了解非摻氣區(qū)的范圍，并及時采取適當(dāng)措施避免非摻氣區(qū)內(nèi)可能產(chǎn)生的空蝕破壞。研究摻氣水流及不同模型模擬的摻氣結(jié)果異同，可以為相關(guān)工程實際提供參考。

用VOF 模型和歐拉雙流體模型計算摻氣率時，以空氣所占體積分數(shù)作為衡量摻氣濃度的參數(shù)。不同數(shù)學(xué)模型計算求得的消力池底部摻氣率結(jié)果如表4所示。

由表4可以看出：①無論大流量還是小流量工況，卷氣模型模擬得到的溢洪道中心線測點底部摻氣率明顯小于歐拉模型。由于卷氣模型基于由紊流產(chǎn)生的擾動、重力及表面張力三種因素，其摻氣原理是因為紊流渦體在流動過程中將表面空氣卷入流體中，主要取決于紊流強度是否克服重力和表面張力組成的表面穩(wěn)定力，在這個過程中紊流是卷氣過程產(chǎn)生的主要因素，所以消力池底部摻氣模擬數(shù)值較小。②在大流量工況下，VOF模型摻氣濃度計算值基本保持在0.24 左右，在消力池內(nèi)摻氣濃度值變化較小，不符合消力池中實際摻氣濃度沿程變化的特征，而歐拉模型則較好地反映了摻氣濃度沿程變化的情況。

表4 底部摻氣率Tab.4 Air entrainment ratio on the bottom

由于在模型試驗中未具體測定消力池的水流摻氣率，所以，只能通過搜集不同參考文獻中水流摻氣模擬數(shù)據(jù)，來側(cè)面驗證模擬結(jié)果的合理性。對于歐拉模型，由參考文獻［5，11］可知，溢洪道摻氣槽臨底0.25m 摻氣濃度均在10%以上；對于卷氣模型，由參考文獻［12，13］可知，底部摻氣率大部分在0～10%左右。根據(jù)多種文獻不同工程中摻氣濃度的計算值，可知本文中計算得到的摻氣濃度值是合理的，采用不同數(shù)學(xué)模型（卷氣模型、VOF 模型、歐拉雙流體模型）的計算結(jié)果可以為相關(guān)工程提供參考。

5 結(jié) 論

本文分別采用CATIA 建立三維模型、ICEM 建立求解所需流體域，應(yīng)用卷氣模型、VOF 模型和歐拉雙流體模型，基于FLOW-3D 和Fluent 軟件對某分期導(dǎo)流工程溢洪道內(nèi)流場進行了三維數(shù)值模擬，計算得到了水流流態(tài)、流速、壓強、紊動能以及摻氣濃度分布規(guī)律?？傮w而言，對于局部流場內(nèi)的流速、壓強、紊動能分布，采用卷氣模型可以較好地模擬實際情況，VOF模型及歐拉模型對于局部壓強和紊動能的模擬結(jié)果存在較大誤差；對于水流摻氣率，歐拉模型的模擬結(jié)果更符合實際。 □