某工程溢洪道摻氣水流流場(chǎng)數(shù)值模擬研究

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

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

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

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

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

1 三維模型

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

用CATIA 建立溢洪道三維幾何模型（圖2），導(dǎo)入ICEM 中，通過(guò)幾何編輯功能完成流場(chǎng)幾何域建立。局部流場(chǎng)部分計(jì)算設(shè)溢洪道長(zhǎng)度為x軸，水流流向?yàn)樨?fù)x方向，寬度方向?yàn)閥軸，水深方向?yàn)閦軸。

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

2.1 紊流模型

本文研究的溢洪道作為導(dǎo)流泄水建筑物，其地形起伏多變，消力池區(qū)域水流流態(tài)復(fù)雜、流線彎曲程度大，故采用RNGk-ε模型［7］。控制方程如下：

連續(xù)方程：

動(dòng)量方程：

紊動(dòng)能k方程：

紊動(dòng)能耗散率ε方程：

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

2.2 卷氣模型

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

紊流渦體的特征尺寸：

單位體積的紊動(dòng)能：

表面穩(wěn)定力：

單位時(shí)間卷氣體積量：

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

2.3 歐拉模型

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

q相的體積Vq為：

且q相的有效密度為：

各相之間體積分?jǐn)?shù)滿足：

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

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

又由q相的動(dòng)量平衡得動(dòng)量方程：

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

式中：ρ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表達(dá)式將相間作用力封閉，其中相間作用力滿足條件

3 三維數(shù)值模擬

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

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

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

計(jì)算工況見(jiàn)表1。

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

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

4.1 流 態(tài)

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

4.2 流 速

由圖6、7 中可知，工況一和工況二中最大流速均出現(xiàn)在摻氣坎（樁號(hào)0-50.8）附近。在各個(gè)測(cè)點(diǎn)中，消力池末端處流速在兩個(gè)工況中均為最小值，進(jìn)一步說(shuō)明水流進(jìn)入消力池后經(jīng)過(guò)摻氣后消耗大量能量，流速在消力池逐漸減小，下泄急流迅速變?yōu)榫徚?，消力池具備很好的消能效果?/p>

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

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

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

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

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

4.4 壓強(qiáng)分布

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

表3 壓強(qiáng)計(jì)算結(jié)果對(duì)比Tab.3 The contrast of pressure

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

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

4.5 紊動(dòng)能分布

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

4.6 水流摻氣率

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

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

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

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

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

5 結(jié) 論

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