摻氣減蝕設(shè)施數(shù)值模擬研究探討

， ，，

(1.中國(guó)電建集團(tuán)貴陽(yáng)勘測(cè)設(shè)計(jì)研究院 工程勘察分院，貴陽(yáng) 550000; 2.長(zhǎng)江科學(xué)院 水力學(xué)研究所，武漢 430010)

1 研究背景

隨著計(jì)算機(jī)和計(jì)算流體力學(xué)的迅猛發(fā)展，數(shù)值模擬技術(shù)以其成本低、信息獲取完整、方案變換靈活等優(yōu)點(diǎn)逐漸成為研究復(fù)雜流動(dòng)問(wèn)題的有效工具，如鄧軍等[1]采用VOF模型對(duì)泄洪洞進(jìn)行了數(shù)值模擬研究，計(jì)算得到了空腔形態(tài)、通風(fēng)設(shè)施進(jìn)氣量及挑距等水力參數(shù)；刁明軍[2]結(jié)合κ-ε模型和VOF模型，對(duì)龍?zhí)ь^明流泄洪洞進(jìn)行了三維紊流數(shù)值模擬，并采用插值技術(shù)在加密后的網(wǎng)格中插入初算的粗網(wǎng)格的計(jì)算成果，獲得較高的計(jì)算精度;張曉東[3]采用歐拉模型研究了陡槽挑坎摻氣水流;覃昕慧[4]采用歐拉模型獲得了摻氣坎后下游斷面的摻氣濃度分布;張宏偉等[5]采用歐拉模型對(duì)陡槽挑坎的摻氣水流進(jìn)行了二維模擬，在相間作用力本構(gòu)方程的構(gòu)建上考慮了湍流擴(kuò)散的影響，詳細(xì)研究了相間作用力對(duì)挑坎下游摻氣濃度分布的影響。

由于摻氣水流問(wèn)題本身的復(fù)雜性，在理論研究方面至今也還沒(méi)有取得突破性進(jìn)展，而物理模型試驗(yàn)受量測(cè)手段、幾何比尺和技術(shù)等條件的限制，研究時(shí)間往往較長(zhǎng)，工作量也較大，并且無(wú)法避免縮尺效應(yīng)的影響。對(duì)摻氣水流進(jìn)行數(shù)值模擬時(shí)，都會(huì)對(duì)相間作用力采用一定的簡(jiǎn)化假設(shè)，因此數(shù)值模擬也只能在一定程度上模擬摻氣水流。

本文采用fluent流體軟件對(duì)某泄洪洞及摻氣設(shè)施水力特性進(jìn)行了數(shù)值模擬，并利用該泄洪洞單體1∶50和1∶30比尺水力學(xué)模型試驗(yàn)成果進(jìn)行了驗(yàn)證。

2 模型試驗(yàn)裝置及試驗(yàn)工況

某泄洪洞工作閘門(mén)后的無(wú)壓段采用“龍落尾”型式，截面為“城門(mén)洞”形式。城門(mén)洞截面為14 m×18 m(寬×高)，側(cè)邊墻高14 m，洞頂拱形半徑為8.12 m，弧度為118.98°的圓弧。無(wú)壓段底板曲線分3段：上游段長(zhǎng)約為275 m、坡度為0.03的緩坡直線段；中間長(zhǎng)約為30 m的曲線連接段，曲線方程為y=0.03x+x2/250；下游長(zhǎng)約為122 m、坡度0.268的陡坡直線段。泄洪洞出口采用挑流方式消能，挑流鼻坎高程為860.00 m。鼻坎下游設(shè)置人工水墊塘消能，塘底高程790 m。如圖1所示。

圖1 泄洪洞明流段模型的體型布置縱剖面圖

該泄洪洞的正常水位泄流量為3 388 m3/s，設(shè)計(jì)水位泄流量為3 535 m3/s，校核水位泄流量為3 723 m3/s。無(wú)壓段長(zhǎng)度將近500 m，最大流速近40 m/s，高速水流產(chǎn)生空化空蝕的可能性較大，為此在無(wú)壓段設(shè)置了跌坎摻氣設(shè)施。通氣孔直徑為1.5 m,為圓形斷面孔。如圖2所示。

圖2 摻氣設(shè)施體型詳圖

本文擬采用CFD方法對(duì)該泄洪洞的無(wú)壓段進(jìn)行數(shù)值模擬計(jì)算，以便獲得與跌坎摻氣水流有關(guān)的水力特性，并將計(jì)算結(jié)果與明流段的1∶30單體水工模型試驗(yàn)數(shù)據(jù)對(duì)比。

3 VOF模型

VOF模型[6]假定兩相或多相流態(tài)之間不發(fā)生質(zhì)量交換，在一個(gè)計(jì)算單元內(nèi)，定義了每相流體之間的體積分?jǐn)?shù)，單元內(nèi)所有流體相的體積分?jǐn)?shù)之和為1。只要知道每一位置、每相流體的體積分?jǐn)?shù)，那么流場(chǎng)內(nèi)的所有特征參數(shù)就可用流體相數(shù)和體積分?jǐn)?shù)的平均值來(lái)描述。在任一給定的單元上，特征參數(shù)代表一相流體還是多相流體，取決于體積分?jǐn)?shù)。如果在一個(gè)單元內(nèi)定義第q相流體的體積分?jǐn)?shù)為αq，則：

αq=0，表示單元內(nèi)第q相流體不存在；

αq=1，表示單元內(nèi)只有第q相流體；

0<αq<1，表示單元內(nèi)存在兩相或多相流體之間的交界面。

根據(jù)局部的αq值，就可以確定計(jì)算域內(nèi)每一控制體上的特征參數(shù)。自由面位置可通過(guò)求解多相流體積分?jǐn)?shù)的連續(xù)方程來(lái)確定。對(duì)于第q相流體，其方程形式[6]為：

求解連續(xù)方程可跟蹤水氣界面。對(duì)于體積分?jǐn)?shù)方程可以用隱式或顯示時(shí)間離散格式求解。經(jīng)過(guò)試算，本文對(duì)體積分?jǐn)?shù)采用隱式時(shí)間離散格式求解。

4 歐拉模型

4.1 體積分?jǐn)?shù)

歐拉模型[6]也稱雙流體模型，可以模擬多相分離流及相互作用的相，相可以是液體、氣體、固體。在多相流模型中，歐拉模型用于流動(dòng)中的每一相。互相貫穿連續(xù)的多相流動(dòng)的描述組成了相體積分?jǐn)?shù)的概念，表示為αq。體積分?jǐn)?shù)代表了每相所占據(jù)的空間，并且每相獨(dú)自地滿足質(zhì)量和動(dòng)量守恒定律。

4.2 守恒方程

守恒方程的獲得可以使用全體平均每一相的局部瞬態(tài)平衡或者使用混合理論方法。在計(jì)算中，可以不考慮相變，同時(shí)忽略表面張力、假設(shè)水氣兩相均不可壓縮、水氣界面局部靜壓平衡，通用的守恒方程形式[7-8]如下：

連續(xù)方程

動(dòng)量方程

4.3 相間質(zhì)量交換

相間質(zhì)量交換Γk是由于界面上的相變而引起的向k相的質(zhì)量傳輸率。對(duì)于水氣二相流而言，水氣兩相間是通過(guò)蒸發(fā)、凝結(jié)、溶解或水中溶解氣體的釋放來(lái)實(shí)現(xiàn)質(zhì)量交換的，在常溫常壓的情況下，這種質(zhì)量交換過(guò)程與高速摻氣水流的空氣混摻過(guò)程相比，完全可以忽略不計(jì)[9]，即Γk=0。

4.4 相間動(dòng)量交換

根據(jù)Ishii[9]的研究，相間動(dòng)量交換Mk可進(jìn)而表示為

式中：上標(biāo)D,L,V,B分別表示靜態(tài)阻力、升力、虛擬質(zhì)量力和Basset力。結(jié)合文獻(xiàn)[5]，在摻氣水流中，阻力和虛擬質(zhì)量力是重要的，本文僅考慮這2種作用力。

5 紊流模型

在標(biāo)準(zhǔn)κ-ε模型中κ和ε是2個(gè)基本未知量，與之對(duì)應(yīng)的輸運(yùn)方程[10-11]為：

連續(xù)方程

動(dòng)量方程

κ方程

Gk+Gb-ρε-YM+Sk;

ε方程

式中：t為時(shí)間；ui和xi分別為速度分量和坐標(biāo)分量；ρ和μ為密度和分子黏性系數(shù)；P為修正的壓力；Gk是由于平均速度梯度引起的紊動(dòng)能κ的產(chǎn)生項(xiàng)；Gb是由于浮力引起的紊動(dòng)能κ的產(chǎn)生項(xiàng)；YM代表可壓縮紊流中脈動(dòng)擴(kuò)張的貢獻(xiàn)；C1ε，C2ε，C3ε為經(jīng)驗(yàn)常數(shù)；σk和σε分別是與紊動(dòng)能κ和耗散率ε對(duì)應(yīng)的Prandtl數(shù)；Sk和Sε是用戶定義的源項(xiàng)。對(duì)于不可壓流體，且不考慮用戶自定義的源項(xiàng)時(shí)，Gb=0，YM=0，Sk=0，Sε=0。在標(biāo)準(zhǔn)κ-ε模型中，根據(jù)Launder等的推薦值及后來(lái)的實(shí)驗(yàn)驗(yàn)證，模型常數(shù)的取值為：C1ε=1.44，C2ε=1.92，C3ε=0.09，σk=1.0，σε=1.3。

選擇紊流模型時(shí)，一方面要測(cè)定其用于各種不同流動(dòng)時(shí)能在不調(diào)整其中的常數(shù)項(xiàng)前提下以多大精度描述流動(dòng)，另一方面還要測(cè)定其處理問(wèn)題所需的時(shí)間及計(jì)算所需的費(fèi)用，前者對(duì)工程應(yīng)用極為重要。目前在工程中應(yīng)用最廣泛的紊流模型為雷諾時(shí)均模型和大渦模擬，雷諾時(shí)均模型又以雙方程κ-ε模型最為成熟，以雷諾應(yīng)力模型最為精確。本文選擇標(biāo)準(zhǔn)κ-ε模型進(jìn)行計(jì)算。

6 計(jì)算方法

6.1 邊界條件

計(jì)算過(guò)程中，入口分為上部空氣入口和下部水入口。其中上部空氣入口為壓力邊界條件，為大氣壓力;下部水入口為由模型率定的入口流量與入口過(guò)水?dāng)嗝婷娣e比值計(jì)算入口速度，依據(jù)計(jì)算工況的不同分別取:V正常=25.45 m/s;V設(shè)計(jì)=26.30 m/s;V校核=27.19 m/s。考慮到實(shí)際的水流運(yùn)動(dòng)過(guò)程在初始狀態(tài)時(shí)通氣孔風(fēng)速和壓力均為0，隨著摻氣的持續(xù)進(jìn)行，風(fēng)速發(fā)生變化，但通氣孔入口壓力變化不大。故兩側(cè)通氣孔處均取為大氣壓力入口。出口邊界也是由上部的氣體和下部的水組成，水流出口一般較平順，本文對(duì)出口邊界采用壓力邊界條件，為大氣壓。整個(gè)過(guò)水道邊界都為固壁邊界，在固壁邊界上取為無(wú)滑移邊界，對(duì)黏性底層采用標(biāo)準(zhǔn)壁面函數(shù)處理。

6.2 網(wǎng)格劃分

該泄洪洞斷面形狀為城門(mén)洞形，由于洞頂部位為弧線，靠近洞頂?shù)木W(wǎng)格質(zhì)量不好控制。因此，為了簡(jiǎn)化和網(wǎng)格劃分的方便，計(jì)算模型用等高的矩形斷面代替。經(jīng)過(guò)試算，將泄洪洞往入口方向延伸30.01 m，即將入口斷面往上游延伸30.01 m，當(dāng)水流運(yùn)動(dòng)到實(shí)際入口斷面時(shí)，水力特性參數(shù)與模型率定的比較接近。

該計(jì)算模型全部采用六面體結(jié)構(gòu)化網(wǎng)格，并在摻氣坎的上下游局部進(jìn)行加密，見(jiàn)圖3。

圖3 網(wǎng)格劃分

6.3 計(jì)算方法

采用有限體積法，隱格式迭代求解。根據(jù)文獻(xiàn)[3]，歐拉模型的氣泡直徑取為1 mm，其相間動(dòng)量交換系數(shù)采用對(duì)稱模型[5]，壓力校正采用Phase Coupled SIMPLE 方法，其它方程采用QUICK離散格式。VOF模型動(dòng)量方程的離散采用二階迎風(fēng)格式，其它為默認(rèn)格式。模擬計(jì)算采用非恒定流算法逼近恒定流穩(wěn)定解，時(shí)間步長(zhǎng)取決于網(wǎng)格尺寸和流速大小，本次計(jì)算取0.000 5～0.005 s。

表1 摻氣設(shè)施的水力特性

7 計(jì)算結(jié)果與分析

計(jì)算過(guò)程有3種工況，依次為工況：Ⅰ(Q=3 388 m3/s)，Ⅱ(Q=3 535 m3/s)，Ⅲ(Q=3 723 m3/s)。

7.1 水面線

泄洪洞無(wú)壓段的水面線見(jiàn)圖4、圖5。

圖4 計(jì)算穩(wěn)定時(shí)水面線沿程分布Q=3 535 m3/s

圖5 實(shí)測(cè)水面線與計(jì)算水面線比較Q=3 535 m3/s(VOF和Eulerian都比較接近)

計(jì)算結(jié)果以水相體積分?jǐn)?shù)為0.5的等值面作為水氣分界面。從實(shí)測(cè)水面線與計(jì)算水面線比較來(lái)看，兩者吻合較好。從計(jì)算結(jié)果可以看出，歐拉模型和VOF模型都能很好地預(yù)測(cè)水面線的沿程分布。

7.2 摻氣坎上水力特性

定義摻氣坎后摻氣水流的氣水比β和摻氣坎上水流的弗勞德數(shù)Fr為：

式中：Qa為通氣量；Qw為水流量；V為摻氣坎上的斷面平均流速；h為摻氣坎上水深。

從表1可以看出，在對(duì)摻氣設(shè)施的水力特性模擬中，2種計(jì)算模型的計(jì)算結(jié)果與實(shí)驗(yàn)值均吻合較好，只是VOF模型計(jì)算的通氣孔風(fēng)速明顯偏大，其原因可能是VOF模型在描述水氣兩相間的作用力和運(yùn)動(dòng)特性時(shí)過(guò)于簡(jiǎn)化。

7.3 空腔特性

此處以Lk表示空腔的長(zhǎng)度；Ln表示水舌內(nèi)緣的挑距。其中空腔長(zhǎng)度以水相體積分?jǐn)?shù)為0.5的等值面作為分界面。

從表2可以看出，歐拉法和VOF方法都能在一定程度上預(yù)測(cè)空腔特性，尤其是水舌內(nèi)緣挑距與實(shí)驗(yàn)值吻合較好。對(duì)空腔負(fù)壓的預(yù)測(cè)也能和實(shí)驗(yàn)值大致吻合。而對(duì)空腔長(zhǎng)度的計(jì)算值均有一定程度的偏大，其原因可能是計(jì)算模型中把空氣作為不可壓縮流體處理，而實(shí)際的空腔內(nèi)是可壓縮的水氣混合物。再者，實(shí)際的摻氣水流各相之間的作用力往往比較復(fù)雜，而2種計(jì)算模型對(duì)相間作用力進(jìn)行估計(jì)時(shí)都做了不同程度的簡(jiǎn)化。

表2 摻氣坎后形成的空腔特性特征參數(shù)

7.4 摻氣濃度

實(shí)驗(yàn)測(cè)量了坎后9個(gè)斷面在渠底中心線的摻氣濃度，因此計(jì)算結(jié)果也取自相應(yīng)位置的摻氣濃度進(jìn)行比較(見(jiàn)圖6，圖7和表3)。

表3 摻氣坎下游臨底水流的摻氣濃度

圖6 摻氣濃度歐拉模型計(jì)算值與實(shí)驗(yàn)值比較

圖7 摻氣濃度VOF模型計(jì)算值與實(shí)驗(yàn)值比較

圖中顯示歐拉模型和VOF模型計(jì)算的坎后摻氣濃度與實(shí)驗(yàn)值均有較大的偏差。

歐拉模型計(jì)算的摻氣濃度沿程變化趨勢(shì)與實(shí)驗(yàn)值大致吻合，造成計(jì)算摻氣濃度偏大的原因可能有：數(shù)值模擬中沒(méi)有采用類似正態(tài)分布的氣泡直徑而是單一氣泡直徑,導(dǎo)致各氣泡的滑移速度相同, 造成氣泡過(guò)于集中；計(jì)算過(guò)程中沒(méi)有考慮升力的影響，因此氣泡的豎向滑移速度偏小, 氣泡上浮速度過(guò)小，致使氣泡過(guò)于集中在水體底部。

VOF模型計(jì)算的摻氣濃度沿程變化趨勢(shì)與實(shí)驗(yàn)值部分吻合，且計(jì)算濃度也偏大，其原因可能有：VOF模型在捕捉氣液界面方面有明顯優(yōu)勢(shì)，但在描述水氣兩相間的力學(xué)和運(yùn)動(dòng)特性時(shí)過(guò)于簡(jiǎn)化，對(duì)擴(kuò)散模型進(jìn)行相間滑移速度本構(gòu)關(guān)系推導(dǎo)時(shí)，由于假設(shè)局部動(dòng)量平衡，因而會(huì)有一定的局限性。

8 結(jié) 語(yǔ)

通過(guò)歐拉模型和VOF隱式時(shí)間離散模型并結(jié)合紊流模型對(duì)某泄洪洞摻氣坎的數(shù)值模擬計(jì)算，可以得到以下結(jié)論：在對(duì)摻氣坎上的水力特性的模擬中，歐拉模型與實(shí)驗(yàn)值非常接近。歐拉模型和VOF隱式時(shí)間離散模型計(jì)算的水舌內(nèi)緣挑距和空腔負(fù)壓與實(shí)驗(yàn)值大致相符。但是計(jì)算的空腔長(zhǎng)度均偏大。 歐拉模型和VOF隱式時(shí)間離散模型都能模擬一定程度的沿程摻氣濃度分布，但是都比實(shí)驗(yàn)值偏大。歐拉模型模擬的沿程摻氣濃度分布與實(shí)驗(yàn)值沿程變化趨勢(shì)比較一致，而VOF隱式時(shí)間離散模型只能部分模擬摻氣濃度實(shí)驗(yàn)值沿程變化趨勢(shì)。

相比于VOF模型，歐拉模型在模擬水氣兩相間的相間作用力和運(yùn)動(dòng)特性時(shí)有很大程度的改善，可以在一定程度上更精確地模擬摻氣水流的坎上水力特性和空腔特性。但是依然可以看到這2種模型在模擬空腔長(zhǎng)度和摻氣濃度沿程變化方面的不足之處，今后要加強(qiáng)對(duì)這方面的研究。

參考文獻(xiàn)：

[1] 鄧 軍，許唯臨.高水頭岸邊泄洪洞水力特性的數(shù)值模擬[J].水利學(xué)報(bào)，2005，(10) :1209-1218．(DENG Jun, XU Wei-lin.Numerical Simulation of Hydraulic Characteristics of High Head Spillway Tunnel[J] .Journal of Hydraulic Engineering, 2005,(10):1209-1218.(in Chinese))

[2] 刁明軍.高壩大流量泄洪消能數(shù)值模擬及實(shí)驗(yàn)研究[D].成都: 四川大學(xué),2004 .(DIAO Ming-jun.Numerical Simulation and Experiment Study on Flood Discharge and Energy Dissipation of High Dam With Large Discharge[D].Chengdu：Sichuan University, 2004.(in Chinese))

[3] 張曉東.泄洪洞高速水流三維數(shù)值模擬[D] .北京：中國(guó)水利水電科學(xué)研究院，2004 .(ZHANG Xiao-dong.Three Dimensional Numerical Simulation of High Velocity Flow in Spillway Tunnels[D].Beijing: China Institute of Water Resources and Hydropower Research,2004.(in Chinese))

[4] 覃昕慧，徐玲君，李國(guó)棟,等.溢洪道摻氣坎二相流數(shù)值模擬[J] .水力學(xué)與水利信息學(xué)進(jìn)展，2009，(3):590-594.(QIN Xin-hui, XU Ling-jun, LI Guo-dong,etal.Numerical Modelling of Two-phase Flows on Aerator in Spillway[J].Advances of Hydraulics and Water Conservancy Information, 2009, (3): 590-594.(in Chinese))

[5] 張宏偉，劉之平，張 東,等.挑坎下游高速摻氣水流的數(shù)值模擬[J].水利學(xué)報(bào),2008,39(12) :1302-1308.(ZHANG Hong-wei, LIU Zhi-ping, ZHANG Dong,etal.Numerical Simulation of Aerated High Velocity Flow Downstream of an Aerator[J] .Journal of Hydraulic Engineering, 2008, 39(12): 1302-1308.(in Chinese))

[6] 江 帆，黃 鵬.Fluent高級(jí)應(yīng)用與實(shí)例分析[M].北京：清華大學(xué)出版社，2008.(JIANG Fan, HUANG Peng.Senior Application and Instance Analysis of Fluent[M].Beijing: Tsinghua University Press, 2008.(in Chinese))

[7] 周立行.湍流兩相流動(dòng)與燃燒數(shù)值模擬[M].北京: 清華大學(xué)出版社, 1991.(ZHOU Li-xing.Numerical Simulation of Two-phase Turbulent Flow and Combustion[M].Beijing: Tsinghua University Press,1991.(in Chinese))

[8] 譚立新.水氣二相流數(shù)值模擬研究[D].成都: 四川大學(xué), 1999.(TAN Li-xin.Mathematical Modelling of Air-water Two-Phase Flows[D].Chengdu：Sichuan University, 1999.(in Chinese))

[9] 張宏偉 .摻氣水流雙流體模型數(shù)值模擬研究[D] .西安：西安理工大學(xué)，2002.(ZHANG Hong-wei.Study on Numerical Simulation of Aerated Water Flow with a Two-fluid Model[D] .Xi’an：Xi’an University of Technology, 2002.(in Chinese))

[10] 付祥釗，龍?zhí)煊?計(jì)算流體力學(xué)[M].重慶:重慶大學(xué)出版社，2007.(FU Xiang-zhao, LONG Tian-yu.Computational Fluid Dynamics[M].Chongqing: Chongqing University Press ,2007.(in Chinese))

[11] 江春波，張永良.計(jì)算流體力學(xué)[M].北京:中國(guó)電力出版社，2007.(JIANG Chun-bo, ZHANG Yong-liang.Computational Fluid Dynamics[M].Beijing: China Electric Power Press, 2007.(in Chinese))