基于IBM-LES的貫流泵湍流特性

闞闞，張清瀅，黃佳程，李昊宇，鄭源，陳龍

(1. 河海大學能源與電氣學院，江蘇 南京，211100； 2. 南通河海大學海洋與近海工程研究院, 江蘇 南通，226019)

貫流泵具有結構緊湊、過流能力強等諸多優(yōu)點，被廣泛應用于中國平原地區(qū)的防洪排澇、農業(yè)灌溉等水利工程和低揚程水力系統(tǒng).近年來，許多學者通過數(shù)值模擬和試驗的方法，對其復雜水力特性進行了深入的研究.石麗建等[1]對考慮回流間隙的全貫流泵裝置進行了計算流體動力學(CFD)數(shù)值模擬，研究發(fā)現(xiàn)，回流間隙對全貫流泵的水力性能的負面影響顯著.謝榮盛等[2]通過數(shù)值模擬研究分析了豎井貫流泵裝置和軸伸貫流泵裝置內部三維流動的特點，總結了不同流道形式下水力性能差異產生的原因.周亞軍等[3]比較了不同進出水流道的水力特性差異，進而提出了豎井式進水流道和直管式出水流道的優(yōu)化方案.楊帆等[4]研究不同葉輪轉速下貫流泵裝置出水流道的流場特征，分析了出水流道進口的平均速度環(huán)量隨流量系數(shù)的變化規(guī)律，得出了轉速對出水流道水力損失的影響.李四海等[5]研究了豎井貫流泵站機組的振動特性，分析了貫流泵裝置在不同工況下不同監(jiān)測點處的壓力脈動特性.ZHU等[6]對導葉葉片固定和可調的貫流泵內部流場進行了數(shù)值模擬，提出了一種導葉可調的貫流泵設計方案.

現(xiàn)有研究關于貫流泵水力性能、時均流動、壓力脈動等宏觀外特性以及平均流動的特點與規(guī)律較多，對貫流泵內部復雜三維湍流的流動特性和湍流結構的深入研究較少.文中以某貫流泵為例，采用大渦模擬(LES)和高精度CFD求解器，以湍動能的角度分析圓柱坐標系下該貫流泵設計工況的湍流特性和湍動能各向異性組成結構，擬為深入了解復雜水力機械內部的湍流流動提供依據(jù).

1 數(shù)值計算方法

計算所用算法描述不可壓縮流體流動的控制方程包括Navier-Stokes方程及連續(xù)性方程[7]，即

(1)

?·u=0,

(2)

式中：u為速度；ρ為密度；p為壓強；μ為動力黏度；S為剪切張量；g為重力加速度；在大渦模擬中，亞格子應力τsgs使用動態(tài)Smagorinsky模型[8]，而在直接數(shù)值模擬中此項為0.

固體在流體中與流體的相互作用通過浸入邊界法(Immersed Boundary Method，IBM)來表達.相對于流體本身的流動，流體中存在的固體會表現(xiàn)出額外的力，此時可以視為動量方程右邊增加了相應的附加項f，即

ρg-?·τsgs]+f.

(3)

文中所使用的為基于水平集(level-set，LS)法的IBM[9-10].LS函數(shù)為一個距離函數(shù).這里定義為每個網格節(jié)點至固體壁面的距離.定義位于流體域中的網格點LS值為正值，位于固體域中的網格點LS值為負值.LS函數(shù)值為0的等值面即描述為固體壁面.為了使得速度場在流固交界面處滿足無滑移邊界條件，固體對流體所施加的額外的力需要準確計算.如圖1所示為IBM在二維情況下的示意圖.圖中紅色點為IB點，即為施加力所位于的網格點.通過IB點在其沿壁面法向投影方向獲得相應固體邊界處的投影點，結合流場中IB點相鄰流體點的速度值，可以插值得到固體邊界處IB點的速度.

圖1 浸入邊界法二維示意圖

圖2所示為壁面函數(shù)示意圖.為保證施加在IB點的力使得IB點切向速度滿足壁面速度規(guī)律，考慮到文中算例雷諾數(shù)較高，結合文中網格分辨率，采用對數(shù)律壁面函數(shù)法進行插值計算IB點的壁面切向速度[7]

圖2 對數(shù)律壁面函數(shù)示意圖

(4)

式中：UIB，UPP分別為IB點和外部流體點切向速度；uτ為壁面摩擦速度；κ為馮卡門系數(shù)；yIB和yPP分別為IB點和外部流體點至壁面的法向距離.而IB點法向速度通過流體點速度與壁面速度進行二次插值獲得

(5)

式中：Un,IB和Un,PP分別為IB點和外部流體點法向速度.

文中貫流泵模型由三維軟件建模后儲存為三角形網格(STL)格式[11]，如圖3所示.為構建精確的LS函數(shù)場，計算網格節(jié)點至固體壁面距離即為計算網格節(jié)點至三角形網格上所有三角形單元的最短距離.

圖3 以STL離散格式網格表達的貫流泵模型

為區(qū)分固體點與流體點上LS函數(shù)值的符號，文中使用射線法[12]辨識網格點的位置位于固體邊界內外.即從每個網格節(jié)點發(fā)出一條隨機射線，計算其穿過固體邊界的次數(shù)，如為奇數(shù)次，則此點為固體點，LS函數(shù)值為負；若為偶數(shù)次，此點為流體點，LS函數(shù)值為正.

2 研究對象及模擬設置

選擇某貫流泵泵段作為研究對象，如圖4所示.流動區(qū)域包括進水管、前置導葉、葉輪、后置導葉、軸及出水管.前置導葉、葉輪、后置導葉的葉片數(shù)分別為5，4，7.運行參數(shù)為設計流量10 m3/s, 額定轉速250 r/min, 最優(yōu)工況點揚程 2.2 m.

圖4 計算域示意圖

建立了基于笛卡爾坐標系(x，y，z)下的三維正交網格，x，y，z方向計算域長度分別為6D，1.1D，1.1D，其中進水管直徑D=1.7 m.葉輪中心位于(2.85D×0.55D×0.55D)，輪轂比為0.388，葉頂間隙為0.01D.3個方向的網格數(shù)分別為1 024，240，240，共計約5 900萬個網格節(jié)點.考慮到固體域并不參與N-S方程求解，實際有效的流體計算域網格節(jié)點約為4 400萬.其中x方向在導葉及葉輪段進行了網格加密，y，z方向使用均勻網格.雷諾數(shù)以進口速度和圓管直徑作為量綱一化后可得，約為7.15×106.量綱一化后的葉輪轉速為11.99.

x方向進口采用均勻入流進口[13-14]，出口采用對流出口.由于包括圓管及葉輪室壁面的流道均由IBM生成，因此y，z方向邊界條件并不對流動過程產生影響，本算例簡單定義為無滑移邊界條件.由于文中泵內流動雷諾數(shù)較高，采用基于對數(shù)律壁模型的大渦模擬方法進行近壁區(qū)的流動求解.

文中計算網格采用交錯網格，N-S方程采用投影法求解，空間上使用二階中心差分格式，時間推進上采用二階龍格庫塔格式.泊松方程在科學計算可移植拓展工具包PETSc上求解完成.程序借助數(shù)據(jù)并行計算平臺信息通過界面(MPI)進行并行計算.計算工作完成于美國圣安東尼瀑布實驗室高性能計算中心平臺，模擬由768核3.06 GHz的英特爾X5675處理器利用MPI并行計算完成，算例總耗時約15 d.外特性的計算結果表明，設計流量工況下模擬的揚程與效率與試驗數(shù)據(jù)吻合良好，誤差分別約為5%和2%，驗證了求解器針對文中算例數(shù)值模擬方法的可行性與計算精度的有效性.

3 結果分析

圖5所示為貫流泵軸向速度(U)和展向速度(V)分布云圖.由圖可見，在均勻來流的條件下，進水管和前置導葉段流動較為均勻，但當旋轉葉輪對水流作用后，尾流表現(xiàn)出了復雜的三維湍流流動特性.顯著的平均速度梯度和湍動能的產生與耗散是尾流中水力損失的主要組成部分.隨著流動中的對流和擴散，水流瞬時速度波動和不均勻度在流動方向下表現(xiàn)出逐漸下降的趨勢.

圖5 貫流泵軸向速度和展向速度分布云圖

考慮到流體機械的軸對稱流動特點，模擬結果轉換至圓柱坐標系下進行分析.計算共有24個葉輪旋轉周期，取后12個葉輪旋轉周期作為后處理數(shù)據(jù).葉輪葉片為4片,文中選取48個葉片旋轉周期進行數(shù)據(jù)處理.這里選擇對靠近輪轂(r=0.25D)、中等跨度(r=0.35D)、靠近流道內壁(r= 0.45D)3個不同徑向距離的圓周截面的展開平面進行分析.湍動能作為湍流強度的度量，對其進行分析可以得到貫流泵湍流流場的強度分布及組成.湍動能公式為

(6)

式中：u′x，u′r，u′θ分別為湍動能k在軸向、徑向和圓周方向的脈動速度，〈 〉代表統(tǒng)計下的時間平均物理量.

圖6所示為靠近輪轂(r=0.25D)圓周截面的湍動能及其分量分布展開圖.圖6a為湍動能k的總量分布，圖6b—6d為k在流向、徑向和圓周方向的組成部分，即為3個方向上的雷諾正應力(下同).由圖可見，靠近輪轂位置附近，k的組成主要由流向方向的分量〈u′xu′x〉貢獻.因為此處流動本身受流向的對流以及輪轂處壁面的剪切影響，都會使流體產生較大的湍流能量，而這部分能量表現(xiàn)為湍流流向方向的劇烈脈動.而k在徑向方向與周向方向的雷諾正應力分量〈u′ru′r〉和〈u′θu′θ〉較小，表現(xiàn)為其脈動強度低于湍動能k在軸向的分量.

圖6 r=0.25D圓周截面湍動能分布

圖7所示為中等跨度(r= 0.35D)圓周截面的湍動能及其分量分布展開圖.此處流動代表了水泵葉輪內主流的湍流特性規(guī)律.可以發(fā)現(xiàn)，k的組成仍主要由流向方向的分量〈u′xu′x〉貢獻.最大值區(qū)域出現(xiàn)在葉輪葉片后方的尾流中.此處的尾流是在相對運動的葉輪葉片后方形成，說明葉片繞流后的尾流里包含了較大的湍流能量且主要由流向方向的分量貢獻.由于在設計工況下運行，流體質點光滑流過固體葉輪與后置導葉葉片表面，葉輪葉片及后置導葉區(qū)域徑向及圓周方向脈動強度較低，雷諾正應力〈u′ru′r〉、〈u′θu′θ〉也隨之較??；當主流流過后置導葉，圓周方向雷諾正應力〈u′θu′θ〉增加，表明后置導葉對主流圓周方向流動的控制同時增加了圓周方向的湍流脈動強度.

圖7 r=0.35D圓周截面湍動能分布

圖8所示為靠近流道內壁(r=0.45D)圓周截面的湍動能及其分量分布展開圖.由圖可見，由于考慮了葉頂間隙的影響，此處的流動較為混亂，湍流具有很大的能量，且3個方向上強度相似.流向的湍動能分量〈u′xu′x〉較大，這是因為間隙處流體在逆壓梯度作用下產生回流和間隙泄漏渦，回流、泄漏渦與主流區(qū)域的速度梯度很大，在此截面表現(xiàn)出很強的湍流運動.而葉頂處流體的平均圓周速度很大，間隙處葉輪葉片無法繼續(xù)對流體在圓周方向做功，此處流體的圓周方向轉動速度無法與主流速度保持一致.同時受間隙渦的影響，湍動能分量〈u′ru′r〉和〈u′θu′θ〉幅值較大.

圖8 r=0.45D圓周截面湍動能分布

圖9所示為后置導葉尾流中流向x=1D，2D，3D(葉輪在流向方向的中心為起點)不同位置的湍流能及其組成周向平均后沿半徑方向的變化規(guī)律，橫坐標方向為沿圓管中心線至圓管壁面.由圖可見，在后置導葉后方的尾流中，湍動能k在3個方向的幅值大小較為平均.湍動能整體及分量沿主流x方向逐漸減弱，表示湍動能能量的逐步耗散與減弱.湍動能k整體沿徑向方向有下降趨勢.在x=2D，3D的靠近壁面處，會較快上升再下降為0，說明隨著流動中速度混合，尾流內湍動能能量耗散，流動逐漸趨向均勻的圓管流動，且壁面處主要由〈u′xu′x〉貢獻.這是因為邊界層內流體正應力主要由黏性正應力提供，而湍流正應力很小(壁面處為0).沿壁面從徑向向圓管中心線看，由于流動雷諾數(shù)較高，邊界層以外流體速度增加很快.一般在均勻的圓管流動中，雷諾應力最大會出現(xiàn)在流場中平均速度梯度最大的位置，因此在靠近壁面處發(fā)展后的湍流與壁面過渡區(qū)域湍動能k會出現(xiàn)最大值.靠近壁面處k主要由〈u′xu′x〉貢獻，這是因為此處主流方向平均速度梯度最大，其次為圓周方向.

圖9 后置導葉尾流不同位置湍流能及其組成的徑向變化規(guī)律

4 結 論

1) 湍動能在靠近葉輪葉片根部及中等跨度的圓周面主要由流向的分量提供.而由于后置導葉對圓周方向平均速度的控制，主要引起了湍流在后置導葉段的圓周方向正應力高幅值區(qū)域.

2) 葉輪葉片葉頂間隙處，由于間隙回流和泄漏旋渦，在間隙渦區(qū)域內流體速度梯度大且流動混亂，湍動能3個方向分量強度均很高.

3) 后置導葉尾流中，湍動能在3個方向幅值相當，且隨著中心線向管壁的徑向和主流方向湍動能強度均隨著距離增大而下降.同時尾流隨著流向距離增加，流動中速度逐漸混合而使得湍動能特性趨向均勻的圓管流動.