黃先北 郭 嬙 仇寶云

(揚州大學電氣與能源動力工程學院， 揚州 225127)

0 引言

開敞式進水池(后文簡稱進水池)是軸流泵站常用的進水形式，主要特征在于水面上方與空氣直接接觸，形成典型的自由界面。受進水池結構以及水位變動的影響，喇叭管淹沒深度往往不足而導致池內存在明顯的漩渦，引起振動、噪聲以及泵效率下降[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)模型，從網格數量級以及模擬的物理時間兩方面進行分析，為進水池吸氣渦的準確模擬提供參考。

1 數學模型

1.1 非常大渦模擬

基于渦粘假設，不可壓流動的雷諾平均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

1.2 S-CLSVOF

本文采用文獻[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

2 進水池算例描述

2.1 幾何與網格

進水池的幾何如圖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)

2.2 邊界條件與數值算法

由于進水池被劃分為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。

3 結果分析

3.1 網格的影響

本文以如下步驟對兩套網格進行計算：基于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作進一步分析。

3.2 計算時間的影響

為檢驗計算時間的影響，基于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)。

3.3 VLES模型特性分析

由前文結果可見，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模式。

4 結論

(1)VLES所得結果與試驗值基本一致，且網格數量級達到O(106)時即可較好地捕捉到吸氣渦的流場結構與程度。

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

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