泵站進水池吸入口渦旋結構及湍流特性的大渦模擬

張德勝 李普熙 趙睿杰 潘 強 施衛東

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

0 引言

泵站進水池是供安裝有水泵的吸水管直接取水的水工建筑物，一般布置在前池與泵房之間或在泵房之下。進水池或引水渠方位設置不當引起的不均勻行進水流、流速梯度較大的剪切流、沿行進水流方向設置幾何體或障礙物引起的旋轉尾流等都會誘發渦旋[1-2]。水泵進水池可能形成的渦旋主要分為自由表面渦和液下渦兩類：前者是由于喇叭口吸力拉拽自由表面的圓柱繞流產生的渦旋而形成的吸氣渦，后者主要是由于喇叭口吸力和水流受壁面形狀約束形成的局部渦流而產生的吸入渦。而這些高度不穩定的渦旋會影響進水池水流的流態，不均勻流態被水泵吸入后，將影響水泵運行效率，甚至導致水泵發生空化現象，產生振動和噪聲，嚴重時會導致整個機組不能正常工作。

進水池中的水流運動具有強烈的三維各向異性湍流特性，由于泵站進水池原型尺寸較大，為研究進水池內部流場，常按一定比例縮放進行模型實驗。文獻[3]在對水泵進水口渦旋方面的實驗研究中，通過采用傳統的染色法顯示渦旋的發生，并應用PIV(Particle image velocimetry)測量技術捕獲渦旋的數目、位置、形狀、大小及強度等定量信息。文獻[4-5]對進水池中自由表面渦和液下渦進行數值模擬，湍流模型采用的是標準k-ε雙方程模型，計算結果可顯示渦旋的大小和位置，其渦量的大小反映渦旋的強度，但對比PIV實驗發現，數值模擬計算出的渦旋偏大，強度偏弱。文獻[6-7]對比了不同湍流模型在泵站進水池渦旋模擬中的適用性。文獻[8]通過PIV實驗研究自由表面渦的形成和演化，并通過數值模擬預測渦旋的位置及結構。文獻[9-10]為更好地了解大渦模擬(LES)對進水池流態的預測性能，將剪切應力傳輸(SST)模型的模擬結果與LES的模擬結果進行比較，LES可以更準確地模擬主流和湍流統計量。文獻[11]用LES對小型水力發電站進水池中具有大渦流的局部區域三維非定常流動進行數值模擬，并用λ2(速度梯度的第二不變量)準則將水電站進水池中的渦三維可視化。文獻[12]為研究軸流泵裝置進水渦旋從生成到潰滅的過程對壓力脈動的影響，用數值模擬及高速攝影實驗對比此過程，揭示喇叭管下方激勵源為渦旋的壓力脈動與葉輪進口斷面處壓力脈動的關系。文獻[13]用VOF(Volume of fluid)方法結合LES對進水池流動進行模擬，并與ADV(Acoustic Doppler velocimetry)實驗結果進行對比。

渦旋結構與周圍流場的湍流特性很難通過實驗獲得，而數值模擬可以提供流場局部細節信息。研究表明，基于RANS湍流模型的數值模擬可以提供相對準確的時均流場特性，而LES方法可以捕獲進水池中更加精細的流動特征。此外，對LES模型的驗證還缺乏系統的分析，對進水池渦流的產生和發展機理也缺乏深入的了解。本文利用LES模型和VOF方法相結合的方法進行數值模擬，并與實驗結果進行比較，討論兩者之間異同。在此基礎上，揭示渦旋周圍各向異性湍流結構，研究渦旋核心向周圍湍流的動量傳遞過程。

1 數值方法

1.1 幾何模型

以PIV實驗的泵站進水池作為研究對象，圖1為進水池三維示意圖。內徑d為0.075 m的吸水管和直徑D為0.115 m的喇叭口垂直安裝在矩形進水池中，進水池長為1.22 m，寬為0.3 m(4d)。喇叭口到底部的距離為0.9d，吸水管軸線距后壁0.9d，吸水管從進水池中吸水，流量恒定為0.004 m3/s，吸水管中的平均流速Up=0.9 m/s，進水流道平均速度Um=0.12 m/s。吸水管內雷諾數Rep大約為75 000。根據實驗數據,基于吸水管直徑和速度的Froude數和Weber數分別約為1.1和840。

圖1 進水池三維示意圖Fig.1 Schematic of pump sump

1.2 控制方程

用六面體結構網格生成計算域，在商用CFD軟件ANSYS CFX 14.5中求解質量和動量守恒方程。流體假定是不可壓縮和等溫的。大渦模擬的基本思想是可解尺度(大尺度脈動)湍流由方程直接數值求解，不可解尺度(小尺度脈動)湍流的質量、動量和能量運輸對大尺度運動的作用采用建立模型(亞格子模型)的方法模擬，從而使可解尺度運動方程封閉。在LES中，過濾尺度Δ通常與局部單元尺寸相同，對小于該方程尺寸的渦流進行方程封閉建模，過濾后的控制方程表示為

(1)

(2)

t——時間，s

ρ——水密度，kg/m3

v——運動粘度，m2/s

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

xi、xj——i、j方向位移

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

對于亞格子尺度(SGS)應力，渦流粘度模型的定義為

(3)

(4)

(5)

(6)

Δ=(ΔxΔyΔz)1/3

(7)

式中τkk——亞格子尺度應力τij中角標相同的應力分量，m2/s2

δij——克羅內克函數

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

Δx、Δy、Δz——網格在3個方向間距的平均值

Cw——自由衰減的各向同性均勻湍流校準常數，取0.5

xk——k方向位移

本文采用壁面適應局部渦粘模型(WALE)求解亞格子應力，該模型可以檢測到與動能耗散相關的所有湍流結構且可以通過線性不穩定增長模式來再現層流到湍流的過渡過程[14]。

1.3 自由表面處理

本文中，液面波動較小，主要研究對象在液面以下，且對自由液面處的網格沿Z軸加密處理，使得VOF方法對界面的捕捉更精確。VOF模型是建立在固定的歐拉網格下的表面追蹤辦法，建立在兩種或者多種流體(或相)不互相混合的前提下。當需要得到一種或多種互不相融的流體交界面時，可以采用VOF模型。在此模型中，不同流體組分共用著一套動量方程，通過引進相體積分數這一變量，實現對每一個計算單元相界面的追蹤。在每個控制容積內，所有相體積分數總和為1。在單元中，若第n相流體體積分數為αn，當αn=0，單元內不存在第n相流體；當αn=1，單元內充滿了第n相流體；當0<αn<1，單元里包含了第n相流體和一相或其他多相流體的界面。跟蹤相之間的界面通過求解容積比率的連續方程來完成。對第n相，有

(8)

(9)

式中V——速度

1.4 邊界條件及設置

泵站進水池中的流態受進口邊界條件影響很大。進口速度與PIV實驗測量的進口截面時均速度相同，模型中進口距離吸水管較遠，水流通過進水池可以得到充分發展。為恒定水位，在出口處指定等于入口的質量流量邊界。在CFX中定義在整個計算域中的初始水位和水中初始靜壓。壁面設置光滑無滑移。頂部設置為開口，使空氣可流入或流出邊界。RNGk-ε湍流模型可以很好地處理帶旋流、高應變率流動及流線彎曲程度較大的流動[15]。為了節約計算時間，首先在約3 800次迭代、速度分量的殘差從初始值降低4～5個數量級后，獲得基于RNGk-ε湍流模型的結果。在CFX 14.5中RNGk-ε湍流模型默認的壁面函數為可縮比例壁面函數(Scalable wall functions)，該壁面函數對于任意細化的網格，能給出一致的解。然后，以該定常結果作為非定常大渦模擬的初始條件。瞬態控制方程的離散采用有限容積法，瞬態項采用二階迎風后插方法，對流項采用中心差分方式。計算時間共50 s，前30 s的時間步長為0.01 s，后20 s的時間步長為0.002 s。兩部分計算的庫朗數均方根分別約為0.8和0.15，說明這兩個時間步長均適用于瞬態模擬。

1.5 網格分辨率和收斂性驗證

實驗中產生的渦旋近似椎體，渦旋中心區域尺度較小，需要比較精細的網格來捕捉渦旋的特征物理量，在會產生渦旋的區域和近壁區對網格進行加密，共劃分3套六面體結構網格，網格數分別為5.0×106、6.8×106和9.5×106，依次在邊界及自由液面處加密，本文采用網格數為6.8×106的模型。所有固體邊界和自由表面處的第1層網格為1 mm，以滿足墻壁處y+為1～3(y+表示第1層網格大小的無量綱量)的要求。網格尺寸從墻壁開始增加，增長率為1.1，以便平滑增加。特別是在渦旋發生的區域中，網格的縱橫比需接近于一致，由于網格比較精細，故在靠近墻壁處可以達到4。而精細的網格則需要大量的計算資源，粗糙的網格無法求出足夠的渦流，故有必要對網格收斂性進行驗證。LES模型和數值誤差對精度的影響相反，相互抵消，導致總誤差較低[16-18]。為了估算計算的數值精度，采用了文獻[19]推薦的既定方法。這種方法是基于Richardson外推理論[20]，并已發展成適用于實際CFD案例的通用公式。經過計算，本文網格解的數值不確定度估計為2.3%。

2 PIV實驗驗證

將模擬結果與PIV實驗作比較，包括自由表面渦、底渦和側壁渦的位置、形狀、大小及渦核強度等方面。實驗中每種渦旋均以0.21 s時間間隔拍攝85幅PIV圖像做時均處理。LES模擬中的時間步長為0.002 s，求解渦旋的瞬時結果。為避免在初期出現任何不確定性，遂將后17 s的模擬結果進行時均處理。選擇底部、側壁和自由表面上3種典型渦旋的流線和時均渦量場來驗證模擬，分別取距離底部和自由表面20 mm截面和距外側壁60 mm截面來對比。總體來說，模擬與實驗具有良好一致性。當比較渦量場時，模擬結果與實驗結果一致，如在底渦周圍會存在一個二次渦，表面渦為兩個旋向相反的渦。通過表1和圖2(ωz表示Z方向的渦量)比較底渦和表面渦渦量發現模擬渦量要略小于實驗值，模擬中渦量更加集中，渦核極值區域較小，渦形狀較實驗略小，并且模擬中的渦旋位置與實驗中的渦旋位置略有不同，這可能歸因于幾個因素。首先，LES模型在極小渦核半徑區域模擬計算會產生偏差，如果局部網格尺寸不夠小，可能會對不準確的值進行建模。另外，邊界條件也存在不確定性，例如入口處速度為PIV拍攝時均速度而實際中入口處流速是不斷變化的，且缺乏入口湍流強度數據。因此，不可能在數值和實驗之間獲得完美匹配。但LES顯示出比其他基于RANS(雷諾時均)的模擬更好的預測能力，如渦旋的位置、形狀、大小以及渦核內的渦量和湍動能的分布[9-10]。總體來說，本文LES與VOF方法模擬結果與實驗結果基本一致，模擬結果可用于深入分析。

表1 渦心渦量(ωzd/Up)Tab.1 Vorticity of vortex center

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

3 結果與討論

3.1 表面渦和底渦與湍流的相互作用

基于LES精確求解的渦流速度場，可以分析湍動能、雷諾應力等湍流信息，揭示渦旋和湍流之間的相互作用，這在實驗或基于RANS的模擬中很難實現。由于附壁渦強度較弱且難以識別，不利于分析，本文只選擇底渦和自由表面渦作為研究對象，揭示渦旋與周圍流場的作用機理。模擬中的渦旋切向速度計算公式為

(10)

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

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

定義壓力最低點為渦心，得到切向速度Vθ在半徑方向的分布，將Vθ峰值對應的半徑定義為渦核半徑，即渦旋特征半徑。

通過λ2等值面將渦旋可視化，如圖3所示，λ2定義為

(11)

圖3a中λ2閾值設為300 s-2，A為自由液面，B為自由表面渦，C為環繞在自由表面渦附近的二次渦。隨著主渦向喇叭口延伸，二次渦螺旋環繞在自由表面渦周圍并與自由表面渦相互作用。圖3b中λ2閾值設為600 s-2，可以看到底渦和環繞底渦的二次渦。λ2閾值的設定是為了使模擬中可視化的渦旋形狀與實驗中所能觀察到的渦旋一致。

圖3 λ2等值面Fig.3 Iso-surfaces of λ2

圖4(T1～T3表示3個連續時刻)為距離自由液面和底面10 mm處，表面渦和底渦在3個連續時刻的切向速度分布，展示了渦旋的減弱過程。將直角坐標系下的雷諾應力轉化為圓柱坐標，揭示渦核區域與相鄰區域的湍流特征，公式為

(12)

其中

ur=uxcosθ+uysinθ

(13)

uθ=-uxsinθ+uycosθ

(14)

式中ux、uy——x、y方向速度分量，m/s

ur——圓柱坐標中徑向速度分量，m/s

uθ——圓柱坐標中切向速度分量，m/s

θ——圓柱坐標系中的極角

圖4 切向速度分布Fig.4 Distribution of Vθ

圖5 隨時間變化的自由表面渦雷諾應力Fig.5 Temporal change of Reynolds stresses for free-surface vortex

圖6 隨時間變化的底渦雷諾應力Fig.6 Temporally changes of Reynolds stresses for floor-attached vortex

3.2 自由表面渦和底渦的空間特性

為研究渦旋的空間特性，創建距自由表面10、20、30 mm的3個平面(Z1～Z3)和距底面50、40、30、20、10 mm的5個平面(Z4～Z8)，如圖7所示。

圖7 截面示意圖Fig.7 Schematic of positions of different planes

圖8 自由表面渦Vθ和Γ的空間分布Fig.8 Spatial distributions of Vθ and Γ for free surface vortex

圖8a和圖9a為自由表面渦和底渦在不同截面處的環量分布，隨著渦旋向喇叭口靠近，環量逐漸匯聚增大。圖8b中，自由表面渦在Z2和Z3平面的渦核半徑稍小于Z1，在Z2和Z3平面渦旋更加集中，渦旋呈現出錐形。底面渦渦核半徑從底面開始先增大后減小，底面渦渦核呈現梭型，這可能是由于無滑移邊界阻礙了流動的圓周運動。

圖9 底渦Vθ 和Γ的空間分布Fig.9 Spatial distributions of Vθ and Γ for floor-attached vortex

圖10 自由表面渦雷諾應力的空間分布Fig.10 Spatial distributions of Reynolds stressses for free surface vortex

圖11 底渦雷諾應力的空間分布Fig.11 Spatial distributions of Reynolds stressses for floor-attached vortex

4 結論

(1)將LES及VOF方法的數值模擬結果與PIV實驗進行對比，結果相差不大，表明模擬結果可以用來進一步分析二階湍流參數，如湍動能和雷諾應力等。

(2)結合LES及VOF方法模擬泵站進水池內的旋渦流動，經過計算，網格解的數值不確定度估計為2.3%。根據實驗對模擬結果進行定性和定量驗證，在渦旋的位置、大小、方向和形狀方面有良好的一致性。