噴水推進泵葉輪空化渦流的數值模擬研究

郭 嬙，王 宇，黃先北, 仇寶云

（揚州大學 電氣與能源動力工程學院，江蘇 揚州 225127）

0 引 言

噴水推進是船舶的一種動力推進方式，其推進原理是利用推進泵噴水獲得反作用力。相比傳統(tǒng)的螺旋槳推進方式，噴水推進泵具有運行平穩(wěn)、適應變工況能力強、可操作性能優(yōu)等特點，在船舶、艦艇領域得到了廣泛應用[1-2]。為適應不同的工作環(huán)境和運行工況，噴水推進裝置的結構形式分為外懸式和內藏式，其中推進泵的類型常用混流式和軸流式。空化是水力機械中難以避免的問題，噴水推進器的空化性能決定了艦船的航速及其安全穩(wěn)定運行特性，其中推進泵的空化性能是影響整個推進裝置性能的重要因素。文獻[3-4]利用高速攝影及成像技術，對噴水推進泵中的空化現象進行了可視化研究，展示了推進泵中的片空化、云空化和渦空化等流動結構及演變規(guī)律。經過實驗結果驗證，數值模擬技術能夠對空化工況下的推進泵水力性能進行合理預測[5-6]，并對其內部空化流場做進一步探究。

噴水推進泵的內部湍流場充滿不同尺度和多種類型的渦流運動，旋渦與空化存在相互作用。Huang等[7]采用數值模擬方法揭示了旋渦區(qū)域與空泡邊界的關系，分析了空化對渦量生成的影響。由于旋轉葉片與泵殼間的相對運動，使得葉頂區(qū)域存在間隙，間隙泄漏流影響推進泵的水力性能[8]，間隙泄漏渦加劇了流場的復雜特性[9]。因此，辨識空化渦流結構、預測空化渦流影響是研究推進泵內流特性的重要內容。本文以某噴水推進泵模型實驗為參考，基于改進的渦空化模型[10]進行流場仿真，將不同類型的渦識別方法應用于推進泵的空化渦流分析，揭示了噴水推進泵的空化渦流特性。

1 噴水推進泵的數值模擬

1.1 推進泵計算模型及網格劃分

參照Chesnakas等[11]的軸流式噴水推進泵的模型實驗，建立推進泵內部流道數值計算模型，如圖1所示，葉輪繞Z軸旋轉，從進口側觀察為逆時針方向。相關幾何參數包括：葉輪外殼直徑D=0.304 8 m，葉頂間隙徑向寬度τ=5×10-4m，葉輪葉片數B=6，下游為8葉片的導葉，推進泵的進出水流道均為圓柱直管，直徑分別為1.0D和0.7D。網格劃分在ANSYS 軟件中完成，根據理查德森外推法進行了網格誤差分析[12-13]，針對葉頂間隙區(qū)的流動研究，網格評判過程中特別考慮了葉頂間隙渦區(qū)域的壓力特征，當網格收斂指標滿足要求時（GCI 系數小于3%），確定網格劃分方案。圖1 所示全流道計算域的總網格節(jié)點數約800萬。為了節(jié)省計算資源，本文采用含單葉片的葉輪和單導葉的流道區(qū)域，進出水段計算域也沿周向按比例縮減，得到單流道數值計算域。單流道計算域中的葉輪外殼、葉片表面及輪轂上的網格示意如圖2所示，為了保證間隙流動和葉片繞流的計算精度，在葉頂及葉片表面附近流場進行了局部網格加密，從葉輪頂端至外殼間徑向的網格數約50，沿葉輪頂端葉片剖面弦長方向的網格數約100，沿葉頂剖面厚度方向的網格數約30，葉頂間隙附近壁面法向y+最大值約為50，符合壁面函數的求解需求。

圖1 噴水推進泵的數值計算域Fig.1 Numerical calculation domain for the water-jet pump

圖2 葉輪外殼和葉片表面與輪轂上的網格Fig.2 Meshes on the rotor casing,blade surface and rotor hub

1.2 數值計算方法

基于雷諾時均方程（RANS）對噴水推進泵的內部流場進行穩(wěn)態(tài)計算，湍流模型采用含旋轉與曲率修正的SST-CC 模型，空化模型采用局部旋轉修正的Zwart模型[10]。參照推進泵的空化實驗，設置葉輪轉速n=33.33 r/s。模型實驗中的流量系數Q*和空化數N*定義為

式中，pV表示水的汽化壓力，設置為3540 Pa；D為葉輪外殼直徑。Q*和N*數值根據實驗工況獲得，可以計算得到若干組體積流量QV（m3/s）和進口總壓pt1（Pa）的數據，用于確定計算域進、出口截面上的邊界條件。將劃分每一個單流道的界面設置為旋轉周期邊界（rotational periodicity），各流道段之間的交界面網格采用GGI聯接方式，固體壁面設置壁面函數，其中葉頂及葉片表面的最大y+值約為150，滿足壁面函數計算要求。數值模擬計算在ANSYS CFX軟件中完成，收斂精度為RMS=10-5。

1.3 計算結果與實驗驗證

當空化數N*=4.0時，將數值模擬得到的推進泵水力特性參數與參考實驗數據對比，根據表1所列公式，計算推進泵的揚程系數H*、功率系數P*和效率η，其中，pt1和pt2分別表示泵進、出口測量斷面的總壓值（Pa），T表示力矩（N·m）。由表1可見，在設計流量Q*=0.85工況下，計算結果與實驗數據偏差小于1%。隨著空化數的進一步降低，葉輪中的空化現象逐漸加劇，設計流量工況下的推進泵水力特性和間隙流場已在文獻[13]中進行了討論。本文沿用文獻[13]驗證的數值計算方法，進行小流量工況的流動研究。如表1所示，當Q*=0.71時，H*、P*和η的偏差在5%以內，相比設計流量工況，小流量工況下的計算偏差略有增大。分析原因，一方面由于偏工況本身流動的復雜性，導致數值預測精度受限；另一方面，由于單流道計算區(qū)域的使用，在一定程度上忽略了流道間可能出現的流動差異。總體看來，數值模擬計算得到的推進泵水力特性參數，滿足用于分析工程實際問題的計算精度需求。

表1 推進泵的水力特性參數（N*=4.0）Tab.1 Parameters of the hydraulic performance for the water-jet pump(N*=4.0)

分析空化對推進泵水力特性的影響。圖3以推進泵的揚程特性為例，顯示了Q*=0.71工況下的揚程系數隨空化數的變化曲線。縱坐標為無量綱處理的揚程系數H*/H0*，其中H0*表示無空化時的揚程系數，在數值計算中采用較大空化數N*=4.0時的結果近似代替。由圖3中的參考實驗曲線可見，隨著空化數的減小，揚程系數先緩慢降低后迅速下降，曲線開始陡降的位置出現在N*=0.9附近。數值計算結果中，當N*＞0.9 時，各工況的揚程系數H*相比H0*的變化幅度不超過0.5%，相比實驗數據偏差在3%以內；當N*＜0.9 時，數值計算的揚程系數曲線陡降。這是由于空化程度的加劇制約了泵的性能，同時加劇了流場的復雜特性，使得數值預測的計算精度受限。本文將在N*＞0.9 范圍內分析推進泵的空化流場。

圖3 不同空化數下的揚程系數變化曲線Fig.3 Variation of the head coefficient with different cavitation coefficients

為了進一步驗證空化流場的計算結果，借助汽相體積分數（α）的分布揭示葉片背面的汽相特征，并與參考實驗中的空化圖像進行對比，如圖4 所示。對應實驗工況，當N*=1.54 時，α＞0 區(qū)域近似反映了葉片背面的空泡影響范圍，其邊緣輪廓呈三角形，覆蓋葉片進水邊和葉頂翼弦約四分之一區(qū)域，α數值沿葉片繞流方向逐漸增大。采用α=0.1的等值面顯示空化的空間區(qū)域，結合分布云圖可見α數值較高區(qū)域對應空泡厚度較大。當空化數減小至N*=1.05 時，空化區(qū)域保持三角形輪廓但面積增大，覆蓋葉頂翼弦約一半區(qū)域；沿葉片繞流方向，空泡尾緣的α值有所減小，這是由于小空化數工況下的葉片表面片狀空化逐漸向云空化狀態(tài)轉變，空泡尾緣受到空泡潰滅的影響，其內部水-汽狀態(tài)使得α值降低。總體看來，采用工程中常用的含汽量α=0.1的等值面，能夠反映與實驗現象一致的空化區(qū)域，數值模擬結果將進一步揭示實驗中未充分展示的空化渦流特征。

圖4 不同空化數下葉片表面的空化流場Fig.4 Cavitation flow field on blade surface with differentcavitation coefficients

2 推進泵內的空化渦流辨識

推進泵內的空化渦流現象集中于葉輪流場中，以空化較為顯著的N*=1.05 工況為例，借助渦流辨識準則研究葉輪內的空化渦流特性。本文分析的渦流辨識準則主要包括以下三類：（1）反映空化特征和速度旋度特征的壓力準則和渦量準則；（2）基于速度梯度張量的Q準則和λ2準則；（3）無量綱化的Ω準則。為了分析葉輪內的空化流場，在葉輪流道中截取垂直于Z軸的軸向截面和葉頂附近的2個展向截面。軸向截面距離葉頂翼型前緣20%弦長位置，展向截面其一位于葉輪頂部，對應100%翼展長度位置，展向截面其二位于葉頂間隙區(qū)域中部，對應葉輪外殼半徑的99.6%位置。由圖5～6可見，含汽量α=0.1的等值面和汽相體積分數云圖均顯示出葉頂空化區(qū)具有楔形輪廓，這是由葉頂翼型繞流和間隙泄漏流形成的渦流影響所致，見圖6中絕對坐標系下的流動顯示，葉頂間隙空化流場中的渦流特性將借助旋渦辨識做進一步分析。

圖5 軸向截面位置及空化顯示（α=0.1）Fig.5 Location of the axial section and display of the cavitation with α=0.1

圖6 展向截面上的汽相體積分數分布和流動顯示Fig.6 Distribution of the vapor volume fraction and flow field on spanwise sections

2.1 壓力準則和渦量準則

空化現象常位于局部壓力降至臨界空化壓力以下的流場區(qū)域，而在旋渦的局部旋轉運動作用下，旋渦中心形成低壓區(qū)，易引發(fā)渦空化現象。旋轉平衡方程[14]（ω×u=-?(p+ρu2/2)/ρ）揭示了高渦量（ω）與低壓力（p）之間的關聯。在選定的軸向截面、葉頂展向截面以及間隙中部展向截面上，分別顯示壓力和渦量云圖，如圖7～8所示。

圖7 葉輪流道內軸向和展向截面上的壓力分布Fig.7 Pressure distribution on axial and spanwise sections in impeller flow passage

在軸向截面圖中，低于常溫水臨界空化壓力（3 540 Pa）的區(qū)域，主要集中于葉片背面，其范圍隨半徑增大而增大，尤其在葉頂附近擴大，這是受到葉頂間隙空化流動的影響。由對應渦量云圖可見，空化區(qū)域的渦量較高，但渦量峰值主要在展向截面云圖中體現。葉頂展向截面的高渦量區(qū)域具有楔形輪廓，并在葉頂翼型中部弦長附近呈“V”形結構，這與圖5 中的空化現象相對應。間隙中部展向截面的高渦量區(qū)域集中在葉頂間隙內部，結合圖6中圓圈區(qū)域內的流動顯示可見，間隙泄漏流近似沿垂直于葉頂剖面翼弦方向流動，促使空化渦流向葉片低壓側發(fā)展。綜合壓力準則和渦量準則的顯示結果可見，空化末端“V”形區(qū)域的壓力值并未全部達到臨界汽化壓力以下，說明空化流場中的渦流特征顯著，旋渦對空化形成有促進影響。

2.2 Q準則和λ2準則

眾多理論渦模型（如Vatistas 模型，Lamb-Oseen 模型等[15]）常用于描述理想旋渦的速度分布，也為水力機械中的渦流研究提供了借鑒。由于旋轉渦核區(qū)沿渦半徑方向存在較高速度梯度，因此，基于速度梯度張量的渦識別準則有助于揭示旋渦區(qū)域。速度梯度張量（?V）可分解為應變率張量（S）和旋轉率張量（Ω），滿足關系式S≡(?V+?VT)/2 和Ω≡(?V-?VT)/2，上標T 表示張量轉置。將速度梯度張量寫成矩陣形式，在矩陣的特征值方程λ3+Pλ2+Qλ+R=0 中，P、Q和R表示速度梯度張量的三個伽利略不變量，Q準則即用第二個不變量Q＞0 的方法表示渦結構。λ2準則是在張量S2+Ω2的三個特征值（λ1≤λ2≤λ3）中，采用λ2＜0的區(qū)域表示渦。

Q準則和λ2準則都在權衡旋轉率與應變率之間的關系，兩準則在大量研究中表現出很強的關聯性，如圖9和圖10所示。雖然兩準則在數值上具有正、負差異，但其絕對值的分布規(guī)律保持一致，絕對值較高的位置表示局部旋轉較強的流動區(qū)域。理論上，Q準則和λ2準則修正了渦量等價于渦的概念，在一定程度上排除了圖8 高渦量區(qū)對剪切流的辨識影響，使得渦流結構特征更為清晰。葉頂展向截面云圖顯示了葉片背面的間隙泄漏渦區(qū)域，泄漏渦沿楔形空化區(qū)的外輪廓發(fā)展，與葉片背面存在夾角。間隙內部較強的渦流特征主要在間隙中部截面云圖中體現。在兩個展向截面中，“V”形渦流區(qū)域相比圖5的空化范圍有所延長，說明空化現象出現于渦流結構中，結合圖7所示的“V”形分支內高于臨界空化壓力的現象分析，充分證實了渦流對空化形成的促進作用。在圖9和圖10的軸向截面上，流場大部分非空化區(qū)域的Q和λ2值都接近于零，僅在葉片背面和葉頂空化區(qū)域包含較強的渦旋運動，也表明了空化對旋渦生成的促進作用。

圖8 葉輪流道內軸向和展向截面上的渦量分布Fig.8 Vorticity distribution on axial and spanwise sections in impeller flow passage

圖9 葉輪流道內軸向和展向截面上的Q分布Fig.9 Invariant Q distribution on axial and spanwise sections in impeller flow passage

圖10 葉輪流道內軸向和展向截面上的λ2分布Fig.10 λ2 distribution on axial and spanwise sections in impeller flow passage

2.3 無量綱的渦識別準則

如圖7～10所示，各準則在辨識葉頂空化渦流區(qū)域時，所取數值相差若干個數量級，閾值問題不可避免。為了克服人為確定合適閾值的問題，Liu 等[17]提出一種無量綱的渦識別方法，命名為Ω渦識別方法。基本原理是將渦量分解為旋轉和剪切部分，通過變量Ω表征旋轉渦量在總渦量中的比例，Ω的計算公式為

式中，a=trace(ATA)，b=trace(BTB)，a和b分別表示速度梯度張量的對稱張量A和反對稱張量B的Frobenius范數的平方，ε為防止分母為零的小正數（本算例中取為1×10-3[s^-2]）。據式可見，Ω的數值介于0到1之間，實現了參數的無量綱化。當反對稱張量B相對于對稱張量A占優(yōu)時，則有Ω＞0.5的常用渦判據。

根據Ω的定義式，采用CEL語句寫入CFX-Post軟件，結果如圖11所示。為凸顯旋轉占優(yōu)區(qū)域，以經驗值Ω=0.52 為參照，將圖中標尺范圍調至0.5～0.55。由兩個展向截面所示的葉頂附近渦流區(qū)域可見，Ω準則顯示出葉頂楔形空化的外輪廓渦流區(qū)，以及空化尾流的“V”形渦流結構，與Q準則和λ2準則所示截面辨識效果基本一致。為了分析渦流的空間結構特征，采用渦準則等值面進行顯示，圖12 比較了Q準則和Ω方法。根據圖9 中的Q準則分布，采用Q=1×106[s^-2]等值面顯示，受到葉輪外壁及間隙附近旋轉剪切影響，葉頂空化外緣出現一個圓柱環(huán)面，圖中進行了透明處理。提高Q值可以消除旋轉剪切對渦流辨識的干擾，局部強旋轉區(qū)域集中在空化尾緣“V”形處，此處壓力場（圖7）特征揭示了旋轉角動量對空化的影響。相比之下，Ω等值面未受到葉頂旋轉剪切的影響，Ω=0.52的經驗判據涵蓋了葉頂空化渦流區(qū)域，同時顯示出葉輪流道內更加豐富的渦流結構。根據眾多算例驗證反饋[18]，Ω經驗判據能夠揭示流場中的強渦和弱渦，借助更大數值可以尋找渦核結構。由圖12 中Ω=0.8 的等值面可見，渦流結構從葉片正面中部區(qū)域向下游發(fā)展，強渦區(qū)域集中在葉片之間。盡管葉道渦附近的壓力未達到空化臨界值，但從對流動影響的角度考慮，葉道渦的影響不可忽略。

圖11 葉輪軸向和展向截面上的Ω 分布Fig.11 Ω distribution on axial and spanwisesections in impeller flow passage

圖12 葉輪流道內渦準則的等值面圖Fig.12 Iso-surface of vortex identification criteria in impeller flow passage

3 結 論

本文以某噴水推進泵模型為例，采用數值模擬計算方法，關注葉頂空化渦流特性，辨識空化流場中的渦流特征，得出以下結論：

（1）針對旋渦空化的特點，采用考慮旋轉修正的湍流模型和空化模型，對推進泵的水力特性和空化區(qū)域進行數值預測，經過與參考實驗結果對比，驗證了數值預報的可行性。在流量系數Q*=0.71 的工況下，分析了空化程度對水力特性的影響，即使空化尚未導致推進泵的性能陡降，其內部流場中的空化現象仍不可忽略。

（2）推進泵內的空化區(qū)域，主要集中在葉輪葉片背面和葉頂間隙處，基于空化區(qū)域的含汽量分析，揭示了葉片表面附著空化和葉頂間隙渦空化的特性差異。采用不同類型的渦流辨識方法，針對葉頂間隙空化渦流場進行分析，綜合壓力、渦量、Q和λ2等常用準則的辨識效果，揭示了葉頂間隙泄漏渦對空化的促進作用。針對渦識別準則普遍存在的閾值選取問題，將無量綱的Ω渦識別方法應用于葉頂間隙空化渦流辨識，體現了Ω經驗值的實用性。基于Ω方法顯示出更為豐富的渦流結構，后續(xù)將對渦流影響空化流場的穩(wěn)定性做進一步探究。