黃先北 郭 嬙 仇寶云
(揚州大學電氣與能源動力工程學院, 揚州 225127)
開敞式進水池(后文簡稱進水池)是軸流泵站常用的進水形式,主要特征在于水面上方與空氣直接接觸,形成典型的自由界面。受進水池結構以及水位變動的影響,喇叭管淹沒深度往往不足而導致池內存在明顯的漩渦,引起振動、噪聲以及泵效率下降[1-4]。若漩渦將水面上方空氣吸入,軸流泵站進水池則形成吸氣渦,進一步加劇了泵的振動與噪聲,同時易導致空化與氣蝕[5-6]。因此,研究進水池吸氣渦意義重大。
分析吸氣渦特征及機理,關鍵在于準確獲取吸氣渦流場結構。隨著計算流體動力學(Computational fluid dynamics,CFD)技術的快速發(fā)展,近年來已有不少研究采用該技術[7-13]對吸氣渦進行數值模擬。影響模擬精度的因素主要有湍流模型、網格數和模擬時間[7-16]。
除上述因素之外,還應考慮吸氣渦現象中水與空氣的交界面捕捉方法。根據現有研究,耦合水平集與流體體積(Coupled level-set and volume-of-fluid, CLSVOF)方法較優(yōu)[10,14,17-18]。
本文以文獻[19]的進水池為研究對象,基于非常大渦模擬(Very large eddy simulation,VLES)模型,從網格數量級以及模擬的物理時間兩方面進行分析,為進水池吸氣渦的準確模擬提供參考。
基于渦粘假設,不可壓流動的雷諾平均Navier-Stokes方程形式為
(1)
式中u——速度矢量,m/st——時間,s
ρ——密度,kg/m3p——壓力,Pa
ν——流體的動力粘度,Pa·s
νt——渦粘系數,Pa·s
S——應變率張量,s-1
St——動量源項,m/s2
考慮到吸氣渦有局部旋轉效應,應在模型中加以體現。為計算νt,文獻[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——湍動能,m2/s2
Pk——湍動能生成率,m2/s3
F1——耦合函數
ω——湍動能比耗散率,s-1
β*——模型系數,取0.09
CDkω——交叉耗散項,s-2
v2——垂直于壁面方向的雷諾正應力,m2/s2
ψ——修正的比耗散率
η——渦粘系數修正率
渦粘系數定義為
(5)
式中a1——模型系數,取0.31
S——應變率張量的模,s-1
F2——模型轉換函數
為了實現混合RANS/LES,基于文獻[15]的方法,將渦粘系數修正為
(6)
其中
(7)
式中γ——模型系數,取0.002
n——模型系數,取2
Lc——截斷尺度,m
Li——湍流積分尺度,m
Lk——Kolmogorov尺度,m
Δx——網格在x方向的尺度,m
Δy——網格在y方向的尺度,m
Δz——網格在z方向的尺度,m
本文采用文獻[21-22]簡化后的版本,命名為Simple coupled level-set and volume-of-fluid(簡稱S-CLSVOF)。為便于編譯該方法,本文選用開源CFD軟件OpenFOAM-2.2.x,其中液相體積分數的輸運方程為
(8)
式中αl——液相體積分數,%
uc——壓縮速度,m/s
采用均相流模型處理液相、氣相,即兩相的速度與壓力相等。
為了將體積分數與動量方程相耦合,將表面張力作為式(1)的動量源項,定義為
(9)
其中
(10)
(11)
γ=1.5xΔ
(12)
式中σ——表面張力系數,N/m
φ——流場中某點至交界面的無量綱距離
xΔ——距離交界面最近的網格尺度,m
進水池的幾何如圖1所示。為便于模擬,將進水池分為進水管域、空氣域以及水域,各域之間采用交界面連接。除圖1所示的幾何參數外,喇叭口直徑為0.15 m,同時,進水管中心偏離進水池中心線,距離兩側壁面分別為0.14 m與0.16 m,水位為0.23 m。
圖1 進水池幾何與計算域劃分Fig.1 Geometry and calculation domain division of intake
為了分析網格數量級對計算結果的影響,本文采用網格數分別為1.61×106(M1)與1.03×107(M2)的兩套網格,如圖2所示,對應的平均y+分別為13~40、9~20。
圖2 進水管附近的網格分布(俯視圖)Fig.2 Grid distribution near intake pipe (bottom view)
由于進水池被劃分為3個域,不同域之間采用任意網格交界面(Arbitrary mesh interface,AMI)[23]進行處理。進口采用均勻來流速度,設定為0.241 5 m/s;出口采用流量出口,設定為0.016 7 m3/s,空氣域上邊界為總壓邊界,等于大氣壓;其余邊界為固壁,采用壁面函數處理近壁面流動。
壓力-速度求解采用PISO(Pressure-implicit with splitting of operators)算法,梯度與散度項采用二階精度的“Gauss linear”格式進行離散,時間離散采用二階精度的“backward”格式,時間步長為1×10-4s。
本文以如下步驟對兩套網格進行計算:基于M1計算14 s(此時吸氣渦達到穩(wěn)定);以上一步的結果作為初始條件,基于M1與M2進行計算,2 000步(即0.2 s)后開始求平均,采用5.5 s內的結果作時間平均。
圖3所示為喇叭口下方15 mm的線段上的速度分布,該線段與喇叭口中心線位于同一平面,如圖3a所示。
由圖3b、3c可見,兩套網格均可捕捉到與試驗[19]較為一致的速度分布規(guī)律。圖3b所示的速度峰值與試驗值之間存在錯位,這與文獻[19,24]中的CFD結果相似,推測該錯位是試驗中存在數值模擬中未體現的外界擾動。
圖3 喇叭口下方的速度分布Fig.3 Velocity distribution below bell mouth
圖4所示為基于兩網格預測的瞬時流線與水的體積分數等值面,其中流線的起點位于水面,等值面根據文獻[12]的研究取αl=95%。為便于分析,圖中還用箭頭顯示了吸氣渦附近流線的總體趨勢。
圖4 瞬時吸氣渦形態(tài)與流線對比(t=19.7 s)Fig.4 Comparison of instantaneous air-entrained vortex shape and streamlines
顯然,盡管位置略有差異,但兩套網格均預測到了水面附近的吸氣渦。此外,與M1不同的是,M2所得吸氣渦更加破碎,這是因為隨著網格密度的增加,原本表征同一體積分數的網格被分成多個更小的網格,在迭代與插值算法的作用下,這些更小網格上的體積分數不再是同一數值,從而導致等值面變得更加零碎。為了更好地顯示吸氣渦,應發(fā)展一種可同時體現水面凹陷(體積分數)與漩渦的吸氣渦識別準則。
為了更好地表征吸氣渦的程度,采用相對吸氣率,定義為
(13)
式中αl,ave——計算域出口的平均液相體積分數,%
圖5所示為相對吸氣率隨時間的變化,為便于顯示,圖中曲線上兩點之間的時間間隔為200時間步,即0.02 s。由圖5可見,M1所得β的變化幅度高于M2,但二者的變化規(guī)律相近。計算時間平均相對吸氣率,M1為2.59×10-4,M2為2.44×10-4,二者差異不大。
圖5 不同網格計算所得相對吸氣率隨時間的變化曲線Fig.5 Variation of relative air-entrainment rate along with time on different meshes
綜合以上分析,網格數量級達到O(106)時可以較好地捕捉到吸氣渦流場的速度分布,同時也可以體現吸氣渦的程度。后文將基于M1作進一步分析。
為檢驗計算時間的影響,基于M1網格計算至60 s。圖6所示為β從t=0 s至t=60 s的變化情況。從圖中可見,吸氣率在0~3 s內逐漸增大,這是因為吸氣渦逐步形成,具體見文獻[12]。此外,由圖6可知,盡管將計算延長至60 s,吸氣率的變化從t=3 s開始并未出現明顯的上升或下降,而是呈現出類似周期性的變化規(guī)律。
圖6 M1計算所得相對吸氣率隨時間的變化曲線Fig.6 Variation of relative air-entrainment rate along with time on M1
為驗證吸氣渦在t=14 s以內是否已達到穩(wěn)定,分析不同時刻下吸氣渦的位置,如圖7所示,吸氣渦位置取自其中心坐標。從圖7可見,預測得到的吸氣渦存在兩個集中區(qū)域,其一位于(-0.06 m,0.08 m)附近,與試驗值一致,但CFD預測到了另一個位于(-0.05 m, 0.15 m)附近的吸氣渦。文獻[19,24]中的CFD結果同樣出現了該現象,推測與上文一致,是因為試驗中存在數值模擬中未體現的外界擾動。對比3~14 s以及14~60 s之間的吸氣渦位置可知,計算時間的延長并未改變吸氣渦的分布規(guī)律。因此,吸氣渦從t=3 s開始處于穩(wěn)定變化的狀態(tài),與圖6吸氣率的變化規(guī)律一致,而延長計算時間不會改變吸氣渦的程度(相對吸氣率)以及位置。
圖7 自由水面上吸氣渦的位置示意Fig.7 Sketch of air-entrained vortex locations on free surface
根據進水池設計標準[25],特定形式的表面渦需至少觀測10 s。因此,對于本算例而言,吸氣渦穩(wěn)定之后再計算10 s即可滿足要求(即t=13 s)。
由前文結果可見,VLESk-ω-v2預測結果與試驗及文獻[19,24]基本一致。為分析該模型的特性,有必要分析吸氣渦附近的Fr與Li,這兩個變量分別代表了解析程度與湍流積分尺度。
圖8所示為自由水面上Li與Fr的分布。由圖8a可見,湍流尺度沿流動方向有降低的趨勢,且渦心附近的值略低于周圍區(qū)域。由圖8b可見,進水管附近的Fr為1,說明此處為RANS模式。對于VLES或其它混合RANS/LES方法而言,降低計算量最有效的方法就是保證近壁區(qū)為RANS模式,因為RANS對近壁面網格要求較低,可以極大地降低近壁面網格密度。
對比圖8a與圖8b可見,二者的分布規(guī)律基本相反,即Li較高的位置則Fr較低。根據式(6),Li越高,則Fr越低。圖8所示的規(guī)律表明,在既定網格下,Fr主要受Li的影響,此時湍流積分尺度越大,則Fr越低,表明湍流的解析度越高。觀察圖8b可見,Fr在流動的上游較低(除壁面附近以外),說明此時湍流基本被解析,而在進水管附近,Fr總體高于上游,這是因為存在流動分離而存在小尺度湍流結構,模型難以完全解析。
圖8 自由水面上瞬時Fr、Li與流線分布(t=19.7 s)Fig.8 Instantaneous Fr and Li and streamline distributions on free surface
總體而言,本文采用的VLES模型在Li的主導下,保證了近壁區(qū)的RANS模式以及湍流核心區(qū)的混合RANS/LES模式。
(1)VLES所得結果與試驗值基本一致,且網格數量級達到O(106)時即可較好地捕捉到吸氣渦的流場結構與程度。
(2)數值計算可以由相對吸氣率隨時間的變化規(guī)律判斷吸氣渦是否達到穩(wěn)定,根據標準要求,在此之后計算10 s即可,進一步延長計算時間并不影響吸氣渦的程度與位置。
(3)既定網格下,VLES模型的解析度主要受湍流積分尺度的影響,近壁面為RANS模式,保證了相對較低的計算量,而在湍流核心區(qū)為混合RANS/LES。