基于LES的泵站前池表面渦及液下渦流瞬態特性分析

潘 強 施衛東，2 趙瑞杰 張德勝

(1.江蘇大學流體機械工程技術研究中心， 鎮江 212013; 2.南通大學機械工程學院， 南通 226019)

0 引言

泵站前池或泵站進水池是供水泵或吸水管直接吸水的水工建筑物，搭配大流量軸流泵常用于防洪抗旱、工農業用水以及大型電廠或核電站的冷卻系統。實驗表明，在泵站進水池內部不同工況的流動中存在多種旋渦，根據誘發位置可分為自由表面渦和液下渦，這些高度不穩定的旋渦會影響水泵吸入口的流態，不僅會造成葉輪載荷的不均勻分布，影響運行效率，甚至引起水泵汽蝕，產生噪聲及震動，嚴重時導致水泵不能正常運行[1-2]。而流態良好的進水池可以保證機組的運行效率和泵站的穩定性。

為研究吸水池的內部流場，對吸水池進行性能評估或優化設計，常按照原型縮放后進行模型實驗。LI等[3]采用二維PIV(Particle image velocimetry)測量了開式進水池的流態，結果表明，進水管周圍流場受來流流速及進水管淹沒深度的影響很大；MANSA等[4]用PIV測量了裝有“T”型消旋器的閉式進水池流場，得到了良好的流態；SUERICH-GULICK等[5-6]通過實驗研究了自由表面渦在水電站進水池內的形成機理，提出了一個半經驗模型去預測渦特性。這些結果與理論的Burger渦模型十分相似，但模型實驗經過縮放后雖然保證了Fr數相似，卻無法滿足Re數相似，在預測湍流及湍流粘度上存在誤差。同時，渦結構與周圍流場的相互作用很難通過實驗來獲得，而CFD(Computational fluid dynamics)提供了更高分辨率的流場信息。CONSTANTINESCU等[7-9]、RAJENDRAN等[10]先后采用標準k-ε、RNGk-ε和SST模型對泵站進水池的內部流動進行了數值模擬，提出了數值計算的可行性，但在瞬態流場及湍流的捕捉上存在誤差；OKAMURA等[11]對比了不同數值計算軟件、網格數、湍流模型的計算結果并與實驗值對比；CHUANG等[12]用流體體積函數(Volume of fluid，VOF)處理自由表面，模擬進水池的流動，并與ADV(Acoustic doppler velocimetry)的實驗結果對比；LUCINO等[13]采用FLOW- 3D模擬出表面渦、壁面渦和附底渦，并在表面渦觀測到了液面凹陷；AKIHIKO等[14]采用λ2準則的等值面將水電站進水池中的渦結構可視化；王福軍[15]對比了不同湍流模型在旋轉流數值計算中的適用性問題；宋希杰等[16]分析了水流壓力脈動的時域特性、頻域特性及進水池底部喇叭管下方的壓力分布，揭示了水流壓力脈動規律及壓力分布與漩渦之間的關系。以上研究結果表明，雷諾時均模型可以得到較為準確的時均流場信息，如速度、渦量與實驗值的吻合度較高，但在瞬態參數的預測上與實驗值有很大出入。相比于RANS(Reynolds-averaged Navier- Stokes equations)方法，大渦模擬(Large eddy simulation，LES)可以提供更為準確、精細的瞬態流場信息[8-9,14,17]，但在泵站進水池的數值研究中仍缺乏LES系統的驗證及分析。

本文采用LES及VOF方法研究泵站進水池內附底渦的時空特性，對數值計算結果進行系統驗證，包括近壁面網格及體網格尺寸、渦流區SGS(Sub-grid scale)模型求解湍動能比例以及與實驗結果對比，并分析數值計算與實驗結果差異產生的原因。基于LES的非定常計算結果，采用λ2準則觀測自由表面渦及附底渦形態，討論旋渦的時均特性和瞬態特性。

1 數值方法

1.1 幾何模型

本文采用的泵站進水池物理模型是RAJENDRAN等[10]進行PIV實驗的模型，如圖1所示。進水池長1.22 m，寬0.3 m，高0.46 m，圖中吸水管內徑d為0.075 m，喇叭口直徑D為0.115 m，對稱布置在與后壁距離為1.4d的位置，泵管內流量0.004 m3/s，管內流速Up為0.9 m/s，水位高度4.5d，進水池橫截面平均流速Um為0.03 m/s。管內Re數75 000左右，進水池內Re數10 000左右，Fr數和Wb數分別為1.1和840。

圖1 進水池三維示意圖Fig.1 3-D schematic of pump sump

1.2 控制方程

在整個計算域生成六面體結構化網格，不可壓縮流體介質的質量和動量守恒方程在ANSYS CFX中求解，LES中的濾波尺度與局部網格尺度一致，比濾波尺度小的渦對流場的影響通過亞格子模型求解，濾波后的控制方程為

(1)

(2)

ρ——密度，kg/m3

ν——運動粘度，m2/s

gi——重力加速度，m/s2

τij——亞格子尺度應力，m2/s2

根據Boussinesq假定τij表達式為

(3)

(4)

(5)

(6)

式中τkk——同向的亞格子應力，m2/s2

δij——克羅內克函數

vsgs——亞格子湍流粘度，m2/s

Δ——網格尺度，m

Cw——WALE常數，取0.5

采用WALE模型求解亞格子應力，該模型克服了Smagorinsky模型耗散過大的問題，可以合理地重現層流及湍流過渡。

1.3 自由表面處理

VOF模型是一種在歐拉網格下的表面追蹤方法，可以有效地模擬出具有明顯界面的兩相流動[12]。通過在水- 空氣交界面求解含有某一相體積分數的控制方程來追蹤運動界面，公式為

(7)

式中x、y、z——坐標軸方向

u、v、w——速度分量，m/s

φ為各相體積分數，當φ為1時則網格中充滿該相流體，當φ為0時則網格中充滿另一相流體，當φ在0～1之間時則含有界面，如此通過φ函數實現對運動界面的追蹤。

1.4 邊界條件及設置

泵站進水池中的流態受進口邊界條件影響很大[2,12]，采用速度進口并保持速度分布與實驗值一致[10]。雖然模擬與實驗難以保證進口處的湍流度一致，但進水池入口與泵管之間的距離足夠湍流的充分發展。為了保證水位恒定，采用流量出口并保證進出口流量差為零。進水池內的初始壓力分布通過CFX軟件中的表達式語言CEL設置。此外，壁面設置為無滑移且光滑，空氣域的頂面設置為opening，允許空氣流出或流入邊界。瞬態控制方程的離散采用有限體積法，對流項采用中心差分，瞬態項采用二階隱式后插法。非定常計算以穩態結果為初始值，時間間隔為0.002s，可以捕捉更為精細的流場信息。

1.5 網格分辨率

對于LES方法而言，網格尺度對計算精度的影響很大，粗糙的網格會導致大部分湍動能通過亞格子模型求解而不是直接求解，精細的網格則需要大量的計算資源。本文所用計算域網格單元數為6.2×106，對近壁面區網格進行了加密，使y+值(近壁面第1層網格的無量綱厚度)在1～3之間，并以1.1的增長率逐漸增加。CELIK等[18]采用相同物理模型、不同網格數的兩套模型來計算LES方法中所用網格尺度的分辨率指數，得到兩套網格中直接求解的湍動能占比，細網格和粗網格中直接求解的湍動能占比公式為

(8)

(9)

p——數值格式精度，取2

α——網格尺度比

采用一套數量為1.2×107的網格用以驗證直接求解的湍動能占比，如表1所示，兩套網格都直接求解了大部分湍動能。

表1 渦流區10個取樣點上求解的湍動能占比Tab.1 Proportion of resolved turbulence kinetic energyat 10 points in representative regions

2 實驗對比驗證

圖2 模擬結果與實驗對比Fig.2 Comparisons of LES and experiment results for three kinds of vortices

數值計算結果與RAJENDRAN等[10]的實驗數據進行了對比，包括自由表面渦及液下渦的位置、形狀和強度。實驗數據包括瞬態及時均流場信息，時均值通過在16 s內平均32幅PIV圖像得到，而非定常數值計算的時間步長為0.002 s，足夠進行實驗對比及流場信息的捕捉。由于自由表面渦和側壁渦對稱分布在泵管兩側，具有相似性，因此只選擇其中一個作為對比。圖2(圖中ωx、ωy、ωz表示渦量沿x、y、z方向的分量)為后壁渦、左側壁渦和自由表面渦瞬態流線圖及時均渦量圖的對比，可以發現，所有渦的瞬態流線與實驗值十分吻合，包括渦的位置及形狀。在時均渦量圖中，模擬出的渦核數與實驗值一致，如后壁渦和左側壁渦的模擬結果與實驗結果中都可看到3個渦核。然而，在渦量時均圖中，渦核位置與實驗值也存在出入，分析原因如下：首先，實驗和數值計算難以保證進口條件完全相同，尤其是進口湍流度無法統一，在RAJENDRAN的實驗中也可以發現，流場極度不穩定；其次，16 s內的32幅PIV圖像只能得到相對的時均信息，在數值計算中也是如此，時間間隔越大得到的時均信息才越準確。總體來看，本文采用LES及VOF方法模擬進水池中的流動，與實驗值結果吻合良好，與RANS方法相比，在渦尺度、位置、形狀、渦量及渦核內湍動能預測方面優勢明顯[9,17]。

3 結果與討論

3.1 渦流時均特性

采用LES方法模擬出旋渦結構，將壓力最低點定義為渦心，得到平均切向速度Vθ沿半徑方向的分布，并將Vθ的極值處半徑定義為渦核半徑，即旋渦特征半徑，公式為

(10)

式中R——以渦心為圓心的任意半徑，m

Γ——在半徑為R的圓內的速度環量，m2/s

圖5 自由表面渦λ2等值面Fig.5 Iso-surface of λ2 for free-surface vortex

由于前水池內液流的高度不穩定性，為得到3種旋渦的時均切向速度及環量分布，在距離自由表面或壁面10 mm的平面上，對16 s內的32個瞬態結果取平均值，該時均方法與實驗中所用方法一致[10]。如圖3和圖4所示，3種渦的切向速度分布均展現出與理論渦模型一致的分布特點。其中，自由表面渦特征半徑及環量較大，而附底渦在渦核內的速度梯度較大，其切向速度沿半徑方向快速增加至極值，表現出較小的旋渦半徑。這是由于附底渦在喇叭口的正下方，較強的抽吸力導致旋渦軸向速度拉伸，形成更加凝聚的渦結構。側壁渦距離喇叭口較遠，抽吸力較弱，表現出十分平緩的速度分布。由于粘性耗散，切向速度在渦核外逐漸降低，速度環量趨向平穩。基于RANS方法的數值計算會假設湍流的各項同性，從而過度地預測渦在半徑方向上的耗散，往往得到較大的渦核半徑及較小的旋渦強度[8]，相比于RANS方法，LES方法得到的旋渦結構更加符合真實流動。

圖3 平均切向速度分布曲線Fig.3 Distribution of tangential velocity for three vortices

圖4 速度環量分布曲線Fig.4 Distribution of circulation for three vortices

3.2 自由表面渦瞬態特性

由于自由表面渦對稱分布在泵管兩側，因此選擇其中之一的左側壁渦，討論其在3個時刻、0.2 s間隔下的瞬態特性。常見的旋渦結構定義方法有Q準則、λ2方法和Δ方法等，本文采用λ2等值面來對旋渦結構進行可視化處理，公式為

(11)

式中λ2——速度梯度張量的二階不變量

通過設定λ2的閾值，可將自由表面渦與周圍湍流區分開。如圖5d(圖中T1、T2、T3表示時間間隔為0.2 s的3個時刻)所示為自由表面渦漸弱過程，λ2等值面的閾值為200 s-2，可以清晰看到泵管邊的自由表面渦(結構B)，自由液面為結構A，可見液面凹陷位置與渦位置保持一致，自由表面渦主渦周邊環繞的二次渦為結構C，來源于泵管壁面的流動分離。隨著液流下沉，結構C呈螺旋狀環繞在主渦周邊并與主渦相互作用，在圖5b和圖5c中可以看到主渦發生了彎曲變形，原因可歸結于液面的波動或二次渦自旋的影響。此外，由于前水池內部流動的高度不穩定性，在圖5中可以看到存在很多雜亂無章的渦流。

在自由液面下10 mm處，將笛卡爾坐標下的渦量換算為圓柱坐標，觀察二維視角下，主渦與周圍二次渦的演變，圓柱坐標的圓心為不同時刻下自由表面的渦心。圖6所示為z軸方向的渦分量ωz，圖7為沿圓周方向的渦分量ωθ，黑色十字表示渦心位置，其長度表示渦核大小。圖6中可以看到，ωz變化不大，而渦核半徑逐漸增大，說明旋渦強度在逐漸減弱，旋渦切向速度變得更為平坦，與圖5一致。圖7展示了主渦附近環繞著的二次渦，在三維視角下呈螺旋狀逐漸靠近主渦并與之相互作用，一方面，二次渦的自旋會引起主渦震蕩，造成主渦彎曲[19]，另一方面，這種相互作用會增加主渦和二次渦的動量耗散，圖5中可以看到，二次渦很難隨主渦下沉向喇叭口運動。由于自由表面的波動往往使渦流難以匯聚，這種二次渦環繞現象在附底渦中更為明顯。

圖6 3個時刻自由表面渦ωz云圖Fig.6 ωz contours of free-surface vortex at three times

圖7 3個時刻自由表面渦ωθ云圖Fig.7 ωθ contours of free-surface vortex at three times

為進一步認識主渦及二次渦的演變，本文通過渦量方程中的對流項及彎曲、拉升項來分析渦量變化[20-21]，公式為

(12)

(13)

式中ω——渦量矢量，s-1

V——速度矢量，m/s

式(12)左邊第2項為對流項，右邊第1項為斜壓作用項，第2項為拉伸、彎曲項，第3項為散度項，第4項為粘性耗散項。由于數值計算是基于不可壓縮流動，式(12)左邊第1項和第3項為零，并且在二維視角下討論渦量的的輸運過程，因此在xy平面上可將式(12)簡化為式(13)。對流和彎曲、拉升可以造成渦量的重新分布，粘性力則會耗散渦量，包括分子粘性及亞格子粘性，但在高雷諾數流動中，粘性力耗散對渦量輸運的影響很小[21]，因此在以下的分析中忽略不計。圖8和圖9分別為對流項及拉升、彎曲項對渦量變化的影響，從式(13)可以看出，對流項正值造成渦量減少，彎曲、拉升項正值使渦量增加。在圖8中，速度對流在渦核內和渦核外對渦量的重新分布都有影響，這表明宏觀流場對自由表面渦強度的變化起一定的作用，而圖9中可以看到，拉升、彎曲項的作用主要存在于渦核內及邊界處，兩者的作用區域不同。對比圖7和圖9可以發現，在二次渦存在的區域，拉伸、彎曲對主渦的渦量改變明顯，說明二次渦的自旋在一定程度上可以引起主渦軸向拉升或者彎曲，從而改變主渦沿z軸方向的渦量。從T1到T3時刻，渦核內部對流項及彎曲、拉伸項的極值區域減小并分散，這與主渦渦核半徑逐漸增大，渦強度減弱相一致。

圖8 3個時刻自由表面渦對流項云圖Fig.8 Convection term of free-surface vortex at three times

圖9 3個時刻自由表面渦拉伸、彎曲項云圖Fig.9 Stretching/tilting term of free-surface vortex at three times

圖10 附底渦λ2等值面Fig.10 Iso-surface of λ2 for floor-attached vortex

3.3 附底渦瞬態特性

由于前水池的幾何對稱性，往往會在底面產生交替出現且旋向相反的兩個渦[10]，在本文中，只分析沿Z軸旋轉的附底渦在3個時刻的瞬態特性。如圖10(圖中Ta、Tb、Tc表示時間間隔為0.2 s的3個時刻)所示，為了更好地呈現附底渦主渦與二次渦的結構，λ2的閾值隨著時間逐漸減小，這也表明了主渦強度在逐漸減弱，與圖10d一致。沒有自由液面波動的影響，呈螺旋狀環繞主渦的二次渦十分明顯，隨著時間的推進，二次渦向喇叭口方向運動，且旋轉軸逐漸轉變為與主渦一致的z軸方向。附底渦周圍的二次渦與主渦的相互作用可以促進主渦的動量耗散，并引起主渦渦量的重新分布，從圖10c可以看到主渦的彎曲。

圖11展示了在距底面10 mm平面上的渦量分布圖，從Ta到Tc渦量極值區域減小，渦強度減弱，在主渦渦核外側，可以看到二次渦沿z軸的渦量ωz，由于二次渦向喇叭口方向運動，其方位角也在不斷變化。圖12為圓周方向的渦量分量ωθ，箭頭所指的二次渦與圖11中的位置相對應。可以看到在渦核外側存在大量的環繞渦，并且旋向呈順時針或逆時針，這個現象與自由表面渦相一致，但環量更加集中的附底渦會壓縮ωθ呈現細長型的分布。圖13和圖14分別為對流項及拉伸、彎曲項對渦量輸運的影響，可以看出在對流作用的影響下，渦核尤其是渦心附近的渦量ωz增加，而彎曲、拉升項引起渦量ωz的減弱，這說明在Ta到Tc過程中，二次渦的作用導致主渦彎曲、拉升是造成主渦ωz降低的主要原因。圖14中的黑色虛線圓為二次渦位置，由于二次渦的旋向逐漸朝向z軸，其沿z軸方向的渦量不斷增加，從Ta到Tc，二次渦逐漸遠離主渦，主渦渦核內由于彎曲、拉升項引起的渦量降低逐漸減弱。

圖11 3個時刻附底渦ωz云圖Fig.11 ωz contours of floor-attached vortex at three times

圖12 3個時刻附底渦ωθ云圖Fig.12 ωθ contours of floor-attached vortex at three times

圖13 3個時刻附底渦對流項云圖Fig.13 Convection term of floor-attached vortex at three times

圖14 3個時刻附底渦拉伸、彎曲項云圖Fig.14 Stretching/tilting term of floor-attached vortex at three times

4 結論

(1)結合LES及VOF方法模擬泵站前水池內的旋渦流動，求解了兩套不同數量的網格，結果表明細網格中直接求解的湍動能占比超過80%，粗網格中直接求解的湍動能占比超過60%。此外，將模擬結果與實驗結果對比，3種典型渦的位置、大小、形狀吻合度良好。

(2)相比于RANS方法會過度預測渦在半徑方向上的耗散，LES方法預測得到的自由表面渦、附底渦及左側壁渦的時均圓周速度、特征半徑和環量分布更加符合真實流動特性。

(3)自由表面渦及附底渦的瞬態特性表明，在主渦渦核邊界附近螺旋環繞著二次渦，且由于自由表面的波動使渦流難以匯聚，這種二次渦環繞現象在附底渦中更為明顯。一方面，二次渦與主渦相互作用，增強主渦動量的向外耗散，另一方面，二次渦的自旋在一定程度上可以引起主渦軸向拉升或者彎曲，導致主渦渦量改變。

1 TSOU J L, MELVILLE B W, ETTEMA R, et al. Review of flow problems at water intake pump sumps[R]. New York: American Society of Mechanical Engineers, 1994.

2 ANSAR M, NAKATO T. Experimental study of 3D pump-intake flows with and without cross flow[J]. Journal of Hydraulic Engineering, 2001, 127(10): 825-834.

3 LI Y, WU Y, MANSA K, et al. The flow research in an open type pump sump by PIV experiments[C]∥ASME 2004 Heat Transfer/Fluids Engineering Summer Conference, 2004: 111-119.

4 MANSA K, ZHANG B, LI X, et al. PIV experimental investigation on the flow in a model of closed pump sump[J]. Tsinghua Science and Technology, 2003, 8(6): 681-686.

5 SUERICH-GULICK F, GASKIN S J, VILLENEUVE M, et al. Free surface intake vortices: theoretical model and measurements[J]. Journal of Hydraulic Research, 2014, 52(4): 502-512.

6 SUERICH-GULICK F, GASKIN S J, VILLENEUVE M, et al. Free surface intake vortices: scale effects due to surface tension and viscosity[J]. Journal of Hydraulic Research, 2014, 52(4): 513-522.

7 CONSTANTINESCU G S, PATEL V C. Numerical model for simulation of pump-intake flow and vortices[J]. Journal of Hydraulic Engineering, 1998, 124(2): 123-134.

8 RAJENDRAN V P, CONSTANTINESCU S G, PATEL V C. Experimental validation of numerical model of flow in pump-intake bays[J]. Journal of Hydraulic Engineering, 1999, 125(11): 1119-1125.

9 TOKYAY T E, CONSTANTINESCU S G. Validation of a large-eddy simulation model to simulate flow in pump intakes of realistic geometry[J]. Journal of Hydraulic Engineering, 2006, 132(12):1303-1315.

10 RAJENDRAN V P, PATEL V C. Measurement of vortices in model pump-intake bay by PIV[J]. Journal of Hydraulic Engineering, 2000, 126(5):322-334.

11 OKAMURA T, KAMEMOTO K, MATSUI J. CFD prediction and model experiment on suction vortices in pump sump[C]∥ The 9th Asian International Conference on Fluid Machinery, 2007:1-10.

12 CHUANG W L, HSIAO S C. Three-dimensional numerical simulation of intake model with cross flow[J]. Journal of Hydrodynamics, Ser. B, 2011, 23(3): 314-324.

13 LUCINO C, GONZALO DUR S. Vortex detection in pump sumps by means of CFD[C]∥XXIV Latin American Congress on Hydraulics, 2010: 21-25.

14 AKIHIKO N, NOBUYUKI H. Large eddy simulation of vortex flow in intake channel of hydropower facility[J]. Journal of Hydraulic Research, 2010, 48(4):415-427.

15 王福軍. 流體機械旋轉湍流計算模型研究進展[J/OL]. 農業機械學報, 2016, 47(2): 1-14.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160201&journal_id=jcsam.DOI: 10.6041/j.issn.1000-1298.2016.02.001.

WANG Fujun. Research progress of computational model for rotating turbulent flow in fluid machinery[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(2): 1-14. (in Chinese)

16 宋希杰,劉超,楊帆,等. 水泵進水池底部壓力脈動特性實驗[J/OL]. 農業機械學報,2017,48(11):196-203.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20171124&journal_id=jcsam.DOI: 10.6041/j.issn.1000-1298.2017.11.024.

SONG Xijie, LIU Chao, YANG Fan, et al. Experiment on characteristics of pressure fluctuationat bottom of pumping suction passage[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2017,48(11):196-203. (in Chinese)

17 TOKYAY T, CONSTANTINESCU S G. Large eddy simulation and Reynolds averaged Navier Stokes simulations of flow in a realistic pump intake: a validation study[C]∥World Water and Environmental Resources Congress, 2005.

18 CELIK I B, CEHRELI Z N, YAVUZ I. Index of resolution quality for large eddy simulations[J]. ASME Journal of Fluids Engineering, 2005, 127(5):949-958.

19 MARSHALL J S, BENINATI M L. External turbulence interaction with a columnar vortex[J]. Journal of Fluid Mechanics, 2005, 540: 221.

20 JI B, LUO X, ARNDT R E A, et al. Numerical simulation of three dimensional cavitation shedding dynamics with special emphasis on cavitation-vortex interaction[J]. Ocean Engineering, 2014, 87: 64-77.

21 JI B, LUO X, ARNDT R E A, et al. Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a NACA66 hydrofoil[J]. International Journal of Multiphase Flow, 2015, 68: 121-134.