基于三方程VLES模型的進(jìn)水池吸氣渦數(shù)值模擬

黃先北 郭 嬙 仇寶云

(揚(yáng)州大學(xué)電氣與能源動(dòng)力工程學(xué)院， 揚(yáng)州 225127)

0 引言

開(kāi)敞式進(jìn)水池(后文簡(jiǎn)稱進(jìn)水池)是軸流泵站常用的進(jìn)水形式，主要特征在于水面上方與空氣直接接觸，形成典型的自由界面。受進(jìn)水池結(jié)構(gòu)以及水位變動(dòng)的影響，喇叭管淹沒(méi)深度往往不足而導(dǎo)致池內(nèi)存在明顯的漩渦，引起振動(dòng)、噪聲以及泵效率下降[1-4]。若漩渦將水面上方空氣吸入，軸流泵站進(jìn)水池則形成吸氣渦，進(jìn)一步加劇了泵的振動(dòng)與噪聲，同時(shí)易導(dǎo)致空化與氣蝕[5-6]。因此，研究進(jìn)水池吸氣渦意義重大。

分析吸氣渦特征及機(jī)理，關(guān)鍵在于準(zhǔn)確獲取吸氣渦流場(chǎng)結(jié)構(gòu)。隨著計(jì)算流體動(dòng)力學(xué)(Computational fluid dynamics，CFD)技術(shù)的快速發(fā)展，近年來(lái)已有不少研究采用該技術(shù)[7-13]對(duì)吸氣渦進(jìn)行數(shù)值模擬。影響模擬精度的因素主要有湍流模型、網(wǎng)格數(shù)和模擬時(shí)間[7-16]。

除上述因素之外，還應(yīng)考慮吸氣渦現(xiàn)象中水與空氣的交界面捕捉方法。根據(jù)現(xiàn)有研究，耦合水平集與流體體積(Coupled level-set and volume-of-fluid， CLSVOF)方法較優(yōu)[10,14,17-18]。

本文以文獻(xiàn)[19]的進(jìn)水池為研究對(duì)象，基于非常大渦模擬(Very large eddy simulation，VLES)模型，從網(wǎng)格數(shù)量級(jí)以及模擬的物理時(shí)間兩方面進(jìn)行分析，為進(jìn)水池吸氣渦的準(zhǔn)確模擬提供參考。

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

1.1 非常大渦模擬

基于渦粘假設(shè)，不可壓流動(dòng)的雷諾平均Navier-Stokes方程形式為

(1)

式中u——速度矢量，m/st——時(shí)間，s

ρ——密度，kg/m3p——壓力，Pa

ν——流體的動(dòng)力粘度，Pa·s

νt——渦粘系數(shù)，Pa·s

S——應(yīng)變率張量，s-1

St——?jiǎng)恿吭错?xiàng)，m/s2

考慮到吸氣渦有局部旋轉(zhuǎn)效應(yīng)，應(yīng)在模型中加以體現(xiàn)。為計(jì)算νt，文獻(xiàn)[20]提出的三方程k-ω-v2模型為

(2)

(3)

(4)

其中

α=0.44-0.116F1
σω=0.856-0.356F1
βω=0.082 8-0.007 8F1
ψ=0.162ω
σk=1-0.15F1

式中k——湍動(dòng)能，m2/s2

Pk——湍動(dòng)能生成率，m2/s3

F1——耦合函數(shù)

ω——湍動(dòng)能比耗散率，s-1

β*——模型系數(shù)，取0.09

CDkω——交叉耗散項(xiàng)，s-2

v2——垂直于壁面方向的雷諾正應(yīng)力，m2/s2

ψ——修正的比耗散率

η——渦粘系數(shù)修正率

渦粘系數(shù)定義為

(5)

式中a1——模型系數(shù)，取0.31

S——應(yīng)變率張量的模，s-1

F2——模型轉(zhuǎn)換函數(shù)

為了實(shí)現(xiàn)混合RANS/LES，基于文獻(xiàn)[15]的方法，將渦粘系數(shù)修正為

(6)

其中

(7)

式中γ——模型系數(shù)，取0.002

n——模型系數(shù)，取2

Lc——截?cái)喑叨?，m

Li——湍流積分尺度，m

Lk——Kolmogorov尺度，m

Δx——網(wǎng)格在x方向的尺度，m

Δy——網(wǎng)格在y方向的尺度，m

Δz——網(wǎng)格在z方向的尺度，m

1.2 S-CLSVOF

本文采用文獻(xiàn)[21-22]簡(jiǎn)化后的版本，命名為Simple coupled level-set and volume-of-fluid(簡(jiǎn)稱S-CLSVOF)。為便于編譯該方法，本文選用開(kāi)源CFD軟件OpenFOAM-2.2.x，其中液相體積分?jǐn)?shù)的輸運(yùn)方程為

(8)

式中αl——液相體積分?jǐn)?shù)，%

uc——壓縮速度，m/s

采用均相流模型處理液相、氣相，即兩相的速度與壓力相等。

為了將體積分?jǐn)?shù)與動(dòng)量方程相耦合，將表面張力作為式(1)的動(dòng)量源項(xiàng)，定義為

(9)

其中

(10)

(11)

γ=1.5xΔ

(12)

式中σ——表面張力系數(shù)，N/m

φ——流場(chǎng)中某點(diǎn)至交界面的無(wú)量綱距離

xΔ——距離交界面最近的網(wǎng)格尺度，m

2 進(jìn)水池算例描述

2.1 幾何與網(wǎng)格

進(jìn)水池的幾何如圖1所示。為便于模擬，將進(jìn)水池分為進(jìn)水管域、空氣域以及水域，各域之間采用交界面連接。除圖1所示的幾何參數(shù)外，喇叭口直徑為0.15 m，同時(shí)，進(jìn)水管中心偏離進(jìn)水池中心線，距離兩側(cè)壁面分別為0.14 m與0.16 m，水位為0.23 m。

圖1 進(jìn)水池幾何與計(jì)算域劃分Fig.1 Geometry and calculation domain division of intake

為了分析網(wǎng)格數(shù)量級(jí)對(duì)計(jì)算結(jié)果的影響，本文采用網(wǎng)格數(shù)分別為1.61×106(M1)與1.03×107(M2)的兩套網(wǎng)格，如圖2所示，對(duì)應(yīng)的平均y+分別為13～40、9～20。

圖2 進(jìn)水管附近的網(wǎng)格分布(俯視圖)Fig.2 Grid distribution near intake pipe (bottom view)

2.2 邊界條件與數(shù)值算法

由于進(jìn)水池被劃分為3個(gè)域，不同域之間采用任意網(wǎng)格交界面(Arbitrary mesh interface，AMI)[23]進(jìn)行處理。進(jìn)口采用均勻來(lái)流速度，設(shè)定為0.241 5 m/s；出口采用流量出口，設(shè)定為0.016 7 m3/s，空氣域上邊界為總壓邊界，等于大氣壓；其余邊界為固壁，采用壁面函數(shù)處理近壁面流動(dòng)。

壓力-速度求解采用PISO(Pressure-implicit with splitting of operators)算法，梯度與散度項(xiàng)采用二階精度的“Gauss linear”格式進(jìn)行離散，時(shí)間離散采用二階精度的“backward”格式，時(shí)間步長(zhǎng)為1×10-4s。

3 結(jié)果分析

3.1 網(wǎng)格的影響

本文以如下步驟對(duì)兩套網(wǎng)格進(jìn)行計(jì)算：基于M1計(jì)算14 s(此時(shí)吸氣渦達(dá)到穩(wěn)定)；以上一步的結(jié)果作為初始條件，基于M1與M2進(jìn)行計(jì)算，2 000步(即0.2 s)后開(kāi)始求平均，采用5.5 s內(nèi)的結(jié)果作時(shí)間平均。

圖3所示為喇叭口下方15 mm的線段上的速度分布，該線段與喇叭口中心線位于同一平面，如圖3a所示。

由圖3b、3c可見(jiàn)，兩套網(wǎng)格均可捕捉到與試驗(yàn)[19]較為一致的速度分布規(guī)律。圖3b所示的速度峰值與試驗(yàn)值之間存在錯(cuò)位，這與文獻(xiàn)[19,24]中的CFD結(jié)果相似，推測(cè)該錯(cuò)位是試驗(yàn)中存在數(shù)值模擬中未體現(xiàn)的外界擾動(dòng)。

圖3 喇叭口下方的速度分布Fig.3 Velocity distribution below bell mouth

圖4所示為基于兩網(wǎng)格預(yù)測(cè)的瞬時(shí)流線與水的體積分?jǐn)?shù)等值面，其中流線的起點(diǎn)位于水面，等值面根據(jù)文獻(xiàn)[12]的研究取αl=95%。為便于分析，圖中還用箭頭顯示了吸氣渦附近流線的總體趨勢(shì)。

圖4 瞬時(shí)吸氣渦形態(tài)與流線對(duì)比(t=19.7 s)Fig.4 Comparison of instantaneous air-entrained vortex shape and streamlines

顯然，盡管位置略有差異，但兩套網(wǎng)格均預(yù)測(cè)到了水面附近的吸氣渦。此外，與M1不同的是，M2所得吸氣渦更加破碎，這是因?yàn)殡S著網(wǎng)格密度的增加，原本表征同一體積分?jǐn)?shù)的網(wǎng)格被分成多個(gè)更小的網(wǎng)格，在迭代與插值算法的作用下，這些更小網(wǎng)格上的體積分?jǐn)?shù)不再是同一數(shù)值，從而導(dǎo)致等值面變得更加零碎。為了更好地顯示吸氣渦，應(yīng)發(fā)展一種可同時(shí)體現(xiàn)水面凹陷(體積分?jǐn)?shù))與漩渦的吸氣渦識(shí)別準(zhǔn)則。

為了更好地表征吸氣渦的程度，采用相對(duì)吸氣率，定義為

(13)

式中αl,ave——計(jì)算域出口的平均液相體積分?jǐn)?shù)，%

圖5所示為相對(duì)吸氣率隨時(shí)間的變化，為便于顯示，圖中曲線上兩點(diǎn)之間的時(shí)間間隔為200時(shí)間步，即0.02 s。由圖5可見(jiàn)，M1所得β的變化幅度高于M2，但二者的變化規(guī)律相近。計(jì)算時(shí)間平均相對(duì)吸氣率，M1為2.59×10-4，M2為2.44×10-4，二者差異不大。

圖5 不同網(wǎng)格計(jì)算所得相對(duì)吸氣率隨時(shí)間的變化曲線Fig.5 Variation of relative air-entrainment rate along with time on different meshes

綜合以上分析，網(wǎng)格數(shù)量級(jí)達(dá)到O(106)時(shí)可以較好地捕捉到吸氣渦流場(chǎng)的速度分布，同時(shí)也可以體現(xiàn)吸氣渦的程度。后文將基于M1作進(jìn)一步分析。

3.2 計(jì)算時(shí)間的影響

為檢驗(yàn)計(jì)算時(shí)間的影響，基于M1網(wǎng)格計(jì)算至60 s。圖6所示為β從t=0 s至t=60 s的變化情況。從圖中可見(jiàn)，吸氣率在0～3 s內(nèi)逐漸增大，這是因?yàn)槲鼩鉁u逐步形成，具體見(jiàn)文獻(xiàn)[12]。此外，由圖6可知，盡管將計(jì)算延長(zhǎng)至60 s，吸氣率的變化從t=3 s開(kāi)始并未出現(xiàn)明顯的上升或下降，而是呈現(xiàn)出類似周期性的變化規(guī)律。

圖6 M1計(jì)算所得相對(duì)吸氣率隨時(shí)間的變化曲線Fig.6 Variation of relative air-entrainment rate along with time on M1

為驗(yàn)證吸氣渦在t=14 s以內(nèi)是否已達(dá)到穩(wěn)定，分析不同時(shí)刻下吸氣渦的位置，如圖7所示，吸氣渦位置取自其中心坐標(biāo)。從圖7可見(jiàn)，預(yù)測(cè)得到的吸氣渦存在兩個(gè)集中區(qū)域，其一位于(-0.06 m,0.08 m)附近，與試驗(yàn)值一致，但CFD預(yù)測(cè)到了另一個(gè)位于(-0.05 m, 0.15 m)附近的吸氣渦。文獻(xiàn)[19,24]中的CFD結(jié)果同樣出現(xiàn)了該現(xiàn)象，推測(cè)與上文一致，是因?yàn)樵囼?yàn)中存在數(shù)值模擬中未體現(xiàn)的外界擾動(dòng)。對(duì)比3～14 s以及14～60 s之間的吸氣渦位置可知，計(jì)算時(shí)間的延長(zhǎng)并未改變吸氣渦的分布規(guī)律。因此，吸氣渦從t=3 s開(kāi)始處于穩(wěn)定變化的狀態(tài)，與圖6吸氣率的變化規(guī)律一致，而延長(zhǎng)計(jì)算時(shí)間不會(huì)改變吸氣渦的程度(相對(duì)吸氣率)以及位置。

圖7 自由水面上吸氣渦的位置示意Fig.7 Sketch of air-entrained vortex locations on free surface

根據(jù)進(jìn)水池設(shè)計(jì)標(biāo)準(zhǔn)[25]，特定形式的表面渦需至少觀測(cè)10 s。因此，對(duì)于本算例而言，吸氣渦穩(wěn)定之后再計(jì)算10 s即可滿足要求(即t=13 s)。

3.3 VLES模型特性分析

由前文結(jié)果可見(jiàn)，VLESk-ω-v2預(yù)測(cè)結(jié)果與試驗(yàn)及文獻(xiàn)[19,24]基本一致。為分析該模型的特性，有必要分析吸氣渦附近的Fr與Li，這兩個(gè)變量分別代表了解析程度與湍流積分尺度。

圖8所示為自由水面上Li與Fr的分布。由圖8a可見(jiàn)，湍流尺度沿流動(dòng)方向有降低的趨勢(shì)，且渦心附近的值略低于周圍區(qū)域。由圖8b可見(jiàn)，進(jìn)水管附近的Fr為1，說(shuō)明此處為RANS模式。對(duì)于VLES或其它混合RANS/LES方法而言，降低計(jì)算量最有效的方法就是保證近壁區(qū)為RANS模式，因?yàn)镽ANS對(duì)近壁面網(wǎng)格要求較低，可以極大地降低近壁面網(wǎng)格密度。

對(duì)比圖8a與圖8b可見(jiàn)，二者的分布規(guī)律基本相反，即Li較高的位置則Fr較低。根據(jù)式(6)，Li越高，則Fr越低。圖8所示的規(guī)律表明，在既定網(wǎng)格下，F(xiàn)r主要受Li的影響，此時(shí)湍流積分尺度越大，則Fr越低，表明湍流的解析度越高。觀察圖8b可見(jiàn)，F(xiàn)r在流動(dòng)的上游較低(除壁面附近以外)，說(shuō)明此時(shí)湍流基本被解析，而在進(jìn)水管附近，F(xiàn)r總體高于上游，這是因?yàn)榇嬖诹鲃?dòng)分離而存在小尺度湍流結(jié)構(gòu)，模型難以完全解析。

圖8 自由水面上瞬時(shí)Fr、Li與流線分布(t=19.7 s)Fig.8 Instantaneous Fr and Li and streamline distributions on free surface

總體而言，本文采用的VLES模型在Li的主導(dǎo)下，保證了近壁區(qū)的RANS模式以及湍流核心區(qū)的混合RANS/LES模式。

4 結(jié)論

(1)VLES所得結(jié)果與試驗(yàn)值基本一致，且網(wǎng)格數(shù)量級(jí)達(dá)到O(106)時(shí)即可較好地捕捉到吸氣渦的流場(chǎng)結(jié)構(gòu)與程度。

(2)數(shù)值計(jì)算可以由相對(duì)吸氣率隨時(shí)間的變化規(guī)律判斷吸氣渦是否達(dá)到穩(wěn)定，根據(jù)標(biāo)準(zhǔn)要求，在此之后計(jì)算10 s即可，進(jìn)一步延長(zhǎng)計(jì)算時(shí)間并不影響吸氣渦的程度與位置。

(3)既定網(wǎng)格下，VLES模型的解析度主要受湍流積分尺度的影響，近壁面為RANS模式，保證了相對(duì)較低的計(jì)算量，而在湍流核心區(qū)為混合RANS/LES。