基于非定常空化流動的離心泵渦旋結構數值分析

伍 杰，邱 寧，朱 涵，徐 佩，司喬瑞

（江蘇大學流體機械工程技術研究中心,江蘇 鎮江 212013）

離心泵是一種重要的流體機械[1]，內部流動復雜且不穩定，時常伴隨著流動分離、失速及回流等現象的發生[2?4]，同時泵內壓力低于汽化壓力還會造成空化的發生[5?6]。因此，泵內的復雜流動結構還需要進一步研究探討[7]。

空化是水力機械中一種常見的現象。在離心泵當中，積聚的空泡能夠堵塞葉輪流道，造成泵的性能下降。Medvitz 等[8]模擬了離心泵的空化流動，發現快速水頭下降現象的物理機制是氣泡爆炸性增長，堵塞了一半以上的葉輪流道。空泡的潰滅還會對葉片表面形成沖擊，造成表面材料的剝蝕。Qiu 等[9]以水翼的研究載體，對空泡破滅時產生的能量沖擊進行了模擬預測，發現最大能量沖擊位于空腔閉合線附近。

在傳統的流場研究中，研究人員以壓力變化作為內部流場分析的重要參數。為了更加深入地了解流場特征，基于速度變化的渦動力學從提出之后便被廣泛應用[10?11]。李等[12]將實驗與模擬相結合，運用正則化螺旋度法，完成了混流泵啟動時瞬時流場中渦核的提取，發現渦旋結構會隨著轉速的增加而出現正反交替現象，但在轉速穩定之后流動也逐漸平穩。Lin 等[13]研究了高壓渦輪的尾跡渦和非定常流動，對周期性尾跡渦和渦輪損失進行了捕捉與分析?？紤]到流量大小對于流動狀態的影響，Zhang 等[14]基于數值模擬方法研究了低比轉速泵的復雜流動，將3 種不同流量工況方案對比發現，隨著流量的降低，泵內流動變得紊亂，在小流量工況條件下已形成明顯的渦流結構。敏等[15]也得出相似結論，流道內的渦量隨著流量的增加而逐漸減小，并且隨著介質流動，旋渦發生破裂并逐漸減小。Posa 等[16]對比研究了在設計工況與非設計工況下帶有不同擴壓器的離心泵性能，發現流動分離和回流現象是由流經葉輪的壓力梯度所決定的。由于非定常流場比較復雜，為了揭示流動原理，Zhang 等[17]研究了離心泵內的非定常流場，將壓力脈動與渦旋結構演化相結合，證明了葉片尾緣脫落渦與壓力脈動存在明顯的相關性。曹等[18]研究發現局部壓力較低的吸力面脫落的旋渦逐漸向上游回流，并發展成為獨立的分離渦。Sun 等[19]研究發現旋轉方向不同的渦旋對會導致空腔周期性脫落，空化的發展能夠增加渦旋強度，而增強的渦旋又反過來促進空腔的脫落。

隨著計算流體力學的發展，有學者發現以往的渦流提取方法存在一些局限性。Liu 等[20?22]基于渦旋的物理特征，將其分為旋轉與非旋轉部分，提出了新Omega 渦識別方法。該方法能夠很好地避免傳統渦識別方法在閾值選取的不確定性。本文基于密度修正模型（DCM）對RNGk?ε湍流模型進行修正，模擬了額定流量、不同空化條件下的非定常流場結構；基于平面速度分量，分析二維渦旋結構的演變以及與空化的相互作用；采用新Omega渦識別方法對三維渦旋結構與空化進行耦合分析。

1 離心泵參數和實驗裝置

1.1 離心泵參數

本文以IS65-50-174 型單級單吸式離心泵為研究載體，該測試泵的計算域三維結構如圖1 所示，設計參數如表1 所示。

表1 測試泵設計參數Tab.1 Design parameters of test pump

圖1 離心泵計算域三維結構Fig.1 3D calculation domain of centrifugal pump

1.2 實驗裝置

本次實驗在江蘇大學流體機械工程技術研究中心進行，實驗臺結構如圖2 所示，測試設備參數如表2 所示。測試系統主要包括水循環回路和數據采集系統兩部分。其中，通過調節測試泵出口管路閥門來控制流量，并基于流量計進行流量數據采集，在進出口管路布置壓力傳感器以監測壓力。

圖2 實驗臺結構Fig.2 Experimental bench structure

表2 實驗設備參數Tab.2 Parameters of experimental equipment

2 數值計算方法

2.1 網格劃分

圖3 所示為采用ANSYS-ICEM 軟件對葉輪流域的網格劃分。與四面體的非結構網格相比，六面體的結構網格具有更少的網格數量，可以提高計算效率。又由于結構網格更為規整，網格控制也更為容易，尤其是葉片壁面附近網格。因此，為了保證數值計算精度的同時還能提高計算效率，測試泵采用六面體結構化網格劃分方法。

圖3 葉輪網格Fig.3 Impeller grid and details

考慮到網格對計算的影響，現選取額定工況進行網格無關性驗證，并以揚程的變化作為網格適用性的評判標準。網格無關性驗證如圖4 所示。根據網格驗證計算結果，最終選擇網格數為2.54×106的方案作為計算所用網格。無關性驗證后的最優網格方案可以在數值計算時降低網格因素對模擬計算結果準確性的干擾，同時還能避免由過多網格數目造成計算資源的浪費。為能夠更好地捕捉到葉片壁面附近的流場變化，需要保證第一層網格高度在適當范圍之內。本文將第一層網格高度設為0.042 mm，網格生長比為1.5，葉片壁面Y+數值小于100，達到計算要求，其分布如圖5 所示。

圖4 網格無關性驗證Fig.4 Grid-independent verification

圖5 葉片壁面Y+分布Fig.5 Yplus distribution of blade wall

2.2 數值計算

在數值計算中，進口邊界條件為壓力進口，其值為101325 Pa，流量出口數值為50 m3/h。葉輪與蝸殼計算域的壁面粗糙度設為0.03 mm，其余壁面均設為無滑移壁面。葉輪設置為旋轉域，轉速2900 r/min，其他部分設為靜止域。將蝸殼、進口延長段分別與葉輪相接的面設為動靜交接面，其余交接面以靜靜交接的方式進行耦合。

由于非定常流動狀態較為復雜，為獲得準確的內部流道模擬結果，本文選用RNGk?ε湍流模型進行模擬計算。該湍流模型可以很好地捕捉高速流動介質的流動狀態，提高了旋渦流動的模擬精度。同時，為更加精確地獲得離心泵內的空化發展情況，考慮氣相和液相混合的可壓縮性，使用DCM 方法進行混合密度修正。湍流黏度定義如下：

式中，N根據文獻推薦取10[23]。

3 結果分析

3.1 實驗驗證與定常分析

本文采用實驗與模擬對比驗證測試泵在額定轉速下不同流量工況的揚程，如圖6 所示。實驗與模擬均呈現出隨著流量的增大而揚程逐漸降低的趨勢。在額定工況時，實驗測得揚程為35.71 m，模擬計算揚程為36.39 m，相對誤差為1.9%，計算結果較為可靠，可以進行進一步的空化計算。

圖6 不同流量下的揚程曲線Fig.6 Flow rate-head curve

圖7 所示為額定工況下離心泵模擬與實驗的空化性能曲線對比。在實驗中，保持流量不變，使用真空泵以抽真空的方式降低進口壓力，從而獲得不同進口壓力下的汽蝕余量。從圖7 中可以看出，在空化初始階段，a點的揚程略有上升，但并不影響計算結果，此時也不會有空化現象的發生，這一點在傳統意義上被定義為非空化條件。隨著進口壓力的降低，b點的揚程開始出現下降趨勢，將該點稱為初始空化條件。當進口壓力降低到c點所對應的壓力時，離心泵已經出現了較為明顯的揚程下降，葉輪當中也會出現明顯的空化。隨著進口壓力的進一步降低，揚程出現急劇下降的趨勢，當揚程下降3%，即空化性能曲線上的d點，普遍認為此時已經處于臨界空化狀態，并將該點稱為臨界空化條件。當進口壓力低于d點時，揚程陡降，并出現嚴重空化，此時也已失去工程意義。

圖7 空化性能曲線Fig.7 Cavitation performance curve

3.2 非定常空化性能分析

若流體的流動狀態不隨時間發生變化，這種流動狀態稱為定常流動。隨著時間不斷變化的非定常流動更貼近自然界中流體的流動狀態。在非定常計算時，將計算總時長設為0.4138 s，即葉輪旋轉運動20 個周期所用時間，時間步長設為2.298 9×10?4s，即葉輪每旋轉4°進行一次計算，每個時間步長內迭代1～10 次，盡可能達到收斂精度，同時葉輪每旋轉20°保存一個結果文件。將數值計算的收斂殘差精度設為1×10?4，以此保證計算結果的可靠性。

圖8 所示為典型時刻3 個不同空化條件下葉輪流道中截面液相速度流線–氣相體積分數分布情況。在NPSHa=3.59 m 條件下，流道截面上的空腔覆蓋面積較小，液相速度流線分布更為均勻合理，而隨著汽蝕余量的逐漸降低，空腔覆蓋面積逐漸增大，流道內的流線也逐漸紊亂起來。尤其是，當NPSHa=2.125 m（臨界空化條件）時，由于流道內的空腔尾部迅速加厚卷曲，擠壓流道空間，造成流道堵塞，從而在卷曲的空腔之后區域形成一個由于液體介質回流而產生的局部旋渦，表現出液相速度流線逐漸匯聚，形成一個或多個流線團。對比不同空化條件下的截面流線，可以看出空腔在較大程度上影響著流道內的介質流動狀態。發生較為嚴重的空化時，葉輪流道下游區域形成較大程度的回流，直接造成泵的性能下降。為進一步分析NPSHa=2.125 m 時的空腔分布情況，圖9 所示為典型時刻不同流道截面氣相體積分數分布。可以發現，更為靠近前蓋板的Span=0.8 截面上空腔覆蓋面積最大，說明空化最為嚴重的區域位于前蓋板附近，且空腔覆蓋面積向后蓋板方向逐漸減小。在Span=0.2 截面，空腔主要分布于葉片吸力面前端，很好地附著于葉片表面，而在吸力面中后部區域還存在已完全卷曲但體積分數較小的空腔。相較于Span=0.2 截面，在Span=0.5 截面上葉片表面空腔附著區域沿著流向方向有較大程度的擴大，在葉片中上游區域空腔很好地附著于吸力面表面，而空腔尾部迅速擴展變厚，并在葉輪旋轉效應作用之下空腔尾部卷起脫離葉片表面，造成流道堵塞，嚴重影響液體介質流動，造成泵的性能下降。在Span=0.8 截面，后端得到充分發展的空腔開始附著于鄰近葉片壓力面，至此，同時連接著吸力面與壓力面的空腔已完全將流道堵住，流動性能完全惡化。

圖8 不同汽蝕余量葉輪流道中截面液相速度流線–氣相體積分數分布Fig.8 Distribution of liquid velocity streamline-vapor volume fraction in the middle section of impeller passage with different NPSHa

圖9 不同葉輪流道截面氣相體積分數分布（NPSHa=2.125 m）Fig.9 Distribution of vapor volume fraction in different impeller passage sections（NPSHa=2.125 m）

3.3 二維渦旋結構與演變分析

在高速旋轉的離心泵中，內部復雜流態監測捕捉較為困難，通常依靠數值模擬的手段完成內部流場分析。渦量作為一個顯著的湍流特征，分析渦量的分布與演變可以更好地理解泵內湍流的非定常流場。式(3)給出了y方向的渦量定義[14]。

式中，vx與vz分別為x、z方向的速度分量。由于本文分析的泵模型繞y軸旋轉，故采用x、z方向的速度分量來進行渦量的定義。

圖10 所示為典型時刻不同流道截面氣相體積分數–渦量分布。同圖9 進行對比可知，圖10 白色標注的區域1 為空泡覆蓋區域，而在其后的藍色標注區域2 為渦量分布區域。從圖中能明顯地看到，更為靠近后蓋板的Span=0.2 截面中，在空腔卷起之后區域存在著許多高強度的正值渦量分布。這是因為該截面中，空腔覆蓋區域較小，流道內的擾動較小，流態相對穩定，所以截面內的高強度渦量更為集中連貫。在Span=0.5 截面中，由于空腔卷曲更為嚴重，流道內的擾動大，流態不夠穩定，集中連貫的高強度渦量破裂成為一些面積較小的分散渦量片。在Span=0.8 截面中，卷曲的空腔已附著于相鄰葉片的壓力面上，截面內的流道已被完全堵住，空腔之后的區域流動極不穩定，高強度渦量進一步分散。附著于相鄰葉片壓力面的空腔之后區域，即黑色標注3 所在區域也分布著面積可觀的高強度渦量，該渦量還有向上游區域進一步發展的趨勢。

圖10 不同葉輪流道截面氣相體積分數-渦量分布（NPSHa=2.125 m）Fig.10 Distribution of vapor volume fraction-vorticity in different impeller passage sections（NPSHa=2.125 m）

圖11 所示為非定常流動狀態下，單枚葉片旋轉1/6 個周期，NPSHa=2.125 m 時離心泵葉輪與蝸殼4 個不同時刻的渦量分布。以圖11（a）中灰白色葉片所在位置t=0 時刻，且葉輪旋轉20°為一個時刻進行截面分析。從圖中可以明顯地看到，葉輪流道內部分布著許多速度為正的正向旋轉渦量，而葉片尾緣為反向旋轉的負值渦量。

圖11 NPSHa=2.125 m 葉輪-蝸殼截面渦量分布及演變Fig.11 Distribution and evolution of vorticity in the impellervolute section of NPSHa=2.125 m

在不同時刻，同蝸殼流道相交接的葉片尾緣附著有尺寸不一的尾跡渦。由于葉片尾緣與蝸殼隔舌的動靜干涉作用，葉片尾跡渦在隔舌附近變化較大。當t=0 時，灰色葉片正好旋轉到隔舌前端，可以很明顯地看到，在動靜干涉作用之下，該枚葉片尾緣處的大尺寸負值尾跡渦開始發生脫落。在t=1/18T時刻，灰色葉片已完全掃過蝸殼隔舌，尾跡渦脫落完成，并附著于隔舌前端的蝸殼壁面。當t=2/18T時，灰色葉片尾跡渦重新附著生長，而隔舌前端渦量也基本消散，在其之后的葉片也攜帶著大量尾跡渦向隔舌靠近。在t=3/18T時刻，灰色葉片之后的葉片已旋轉到隔舌前端，并開始重復t=0時刻在動靜干涉作用下而發生的尾跡渦脫落現象。

對于4 個不同時刻，葉片吸力面中上游區域存在一個緊貼葉片壁面的負值渦量帶。其后區域還有著正值渦量的分布，且該正值渦量帶很快便發生卷曲而脫離葉片表面，并隨著葉輪的旋轉開始出現渦帶脫落，然后在流道內逐漸消散。

觀察圖12、13 可以發現，從演變起始時刻t=0開始，流道內的渦量隨著葉輪的旋轉而呈現出逐漸減少的趨勢，這與NPSHa=2.125 m 的渦量演變相對應。同時，有一個比較有意思的現象是尾跡渦附近還存在有些許正值高強度渦，并隨著葉片的旋轉而逐漸生長發展，最后脫落附著在蝸殼壁面。將3 個不同的空化條件對比可知，空化程度較小的工況，流道內的高強度渦量分布更少，流態更為簡單合理。

圖12 NPSHa=2.624 m 葉輪-蝸殼截面渦量分布及演變Fig.12 Distribution and evolution of vorticity in the impellervolute section of NPSHa=2.624 m

圖13 NPSHa=3.59 m 葉輪-蝸殼截面渦量分布及演變Fig.13 Distribution and evolution of vorticity in the impellervolute section of NPSHa=3.59 m

圖14 給出了NPSHa=2.125 m 時，一個旋轉周期內典型葉片不同展長基于 ?y的渦量分布。從圖中可以看到，在葉片起始部分均存在一個渦量驟變的過程，并且不同展長位置的起始渦量還存在一定差異。這是由于來流首先沖擊葉片前緣，從而在葉片前緣區域形成局部渦旋，而span=0.8 相對于span=0.2 與span=0.5 受到的沖擊更大，所以在span=0.8上形成最大起始渦量。葉片吸力面前半段有高強度渦量分布，說明在流道上游區域流動變化較大，并且流動情況復雜，同時空化也容易在該區域內生長并發展。

圖14 NPSHa=2.125 m 一個旋轉周期內葉片吸力面上渦量分布Fig.14 Distribution of vortex on the suction surface during one rotation cycle with NPSHa=2.125 m

隨著時間的演變，葉片前半段的渦量出現正負交替，形成一個周期性的變化過程，而后半段的渦量數值在小范圍內變化。由于空腔在葉片中段位置卷起，而渦旋隨著卷曲的空腔脫離葉片表面，在葉片中段位置呈現出驟變的過程。在L=0.3～0.6 區間內，不同展長的空腔體積存在較大差異，所以在此區域的渦量出現較為明顯的波動。

3.4 三維空化– 渦結構分析

渦是一種流體在旋轉過程中產生的物理現象，是一個具有方向與大小的矢量，但以往常用的識別方法得到的渦是一個只有大小而沒有方向的標量，不能完整地表現出渦的具體參數與物理特征。Liu 等[20]提出了Omega 渦識別方法，主要是將流場中的渦分解為旋轉部分與非旋轉部分來進行渦特征的體現。

式中：R代表旋轉部分渦量；S代表非旋轉部分渦量，在一般情況下R與S具有不同方向。

為了更好地表征具有旋轉效應的渦量場，將表征旋轉渦量與流場總渦量之比定義為 ?，而該定義參數 ?在0～1 范圍內取值?？梢钥闯觯S著數值的增大，所表征的旋轉渦具有更高的旋轉效應，計算公式如下：

為更好地識別流場中的渦旋結構，選取反對稱張量B大于對稱張量A，即 ?>0.5 時作為流場渦識別的判據條件。同時，Liu 等[20]在文獻中給出了?=0.52的參考值來進行識別渦旋結構邊界。與該Omega 渦識別方法相比，Q準則等一些常用的渦識別方法具有更大的取值范圍，合適的渦識別參數選取較為不易。

圖15 所示為NPSHa=2.125 m 時葉輪流道內氣相體積分數為0.1 的空腔與?=0.52的渦旋結構分布。從圖中可以看到，在流道上游區域，渦旋結構較為穩定，且被空腔覆蓋區域很好地包裹起來。由于臨界空化狀態的空腔結構比較穩定，空腔內部流動擾動較小，從而空腔內部形成的渦旋與空腔有著相似的穩定結構。在葉輪旋轉效應的作用下，葉片中部位置的空腔與渦旋結構均卷曲脫離葉片表面。在葉輪流道下游區域，由于沒有了穩定的空腔作為渦旋的外部覆蓋結構，渦旋結構變得不再穩定，從而形成許多隨著介質流動而不斷變化的小體積渦。

圖15 NPSHa=2.125 m 空腔-渦旋結構分布Fig.15 Distribution of cavity-vortex structure with NPSHa=2.125 m

三維空腔與渦旋結構的提取，能夠更好地對上述二維渦量的演變做出進一步的解釋。在空化的作用影響之下，葉片吸力面的渦量表現出同空腔相似的變化規律，沿著葉片弦長方向，吸力面的渦量大小隨著表面附著空腔的增厚而逐漸變大。在葉片中段位置，由于空腔卷曲引起渦旋結構脫離葉片表面，從而產生葉片表面渦量驟變的現象。由于流道下游區域的渦旋結構零散分布于其中，并沒有完整地附著于葉片表面，因而在吸力面后半段的渦量沒有太大波動，數值較小。

圖16—17 分別對應NPSHa 為2.624 m 及3.59 m時的空腔– 渦旋結構。不同于臨界空化工況，其余空化條件下渦旋結構較為穩定集中，在流道下游區域沒有出現零散分布的現象。隨著NPSHa 的增加，空腔體積減小，沒有占據過多的葉輪流道空間，擾動減小，流態更為穩定且合理，流道內的渦旋結構也隨之逐漸減少。

圖16 NPSHa=2.624 m 空腔-渦旋結構分布Fig.16 Distribution of cavity-vortex structure with NPSHa=2.624 m

圖17 NPSHa=3.59 m 空腔-渦旋結構分布Fig.17 Distribution of cavity-vortex structure with NPSHa=3.59 m

4 結論

本文利用數值模擬方法對離心泵內的非定?？栈蜏u旋結構進行了耦合分析，得到以下結論。

1）隨著汽蝕余量的降低，離心泵內空腔逐步發展，臨界空化條件下空腔結構卷曲并堵塞流道，在葉片吸力面出口區域產生回流，從而造成流動紊亂，泵的性能下降。

2）基于x、z速度分量，對葉輪與蝸殼截面進行渦量提取分析，可知隨著葉輪的旋轉，在葉片尾緣同蝸殼隔舌的動靜干涉作用之下，葉片尾跡渦發生周期性脫落，形成的脫落渦附著在隔舌前端壁面，同時還能捕捉到在臨界空化條件下由于渦量卷曲而產生的驟變過程。

3）采用新Omega 渦識別方法能對非旋轉部分渦量進行篩選過濾，能很好地捕捉流道上游區域空腔內部穩定清晰的三維渦旋結構。對比不同空化條件可知，隨著汽蝕余量的降低，流道內的渦旋結構逐漸增加。