非線性空泡模型在軸流泵空化模擬中的評估分析

洪 鋒 ，張?zhí)煊?，蘇華山 ，袁建平 ，張 帆 2，

（1.三峽大學(xué) 機械與動力學(xué)院，湖北宜昌 443002；2.三峽大學(xué) 水電機械設(shè)備設(shè)計與維護湖北省重點實驗室，湖北宜昌 443002；3.江蘇大學(xué) 國家水泵及系統(tǒng)工程技術(shù)研究中心，江蘇鎮(zhèn)江 212013）

0 引言

空化是軸流泵內(nèi)部不可避免的一種水動力學(xué)現(xiàn)象，空化嚴重時會破壞過流部件表面和產(chǎn)生振動噪聲，嚴重影響機組的高效安全運行。近年來，基于Navier-Stokes方程和相間質(zhì)量傳輸模型框架的空化數(shù)值計算方法，由于能描述空化的動量傳輸和質(zhì)量傳輸兩個過程，成為空化流動研究手段的主要趨勢。Athavale等［1］應(yīng)用Singhal模型對3種不同比轉(zhuǎn)速船用軸流泵空化流場進行了數(shù)值模擬，發(fā)現(xiàn)不可凝結(jié)氣體（non-condensable gas，NCG）能夠降低泵的揚程及影響空化的發(fā)生和發(fā)展，還指出了該空化模型中蒸發(fā)與凝結(jié)系數(shù)限制了其應(yīng)用。此外，Zwart、Schnerr-Sauer空化模型在模擬單個水翼、螺旋槳與葉片泵葉片表面上的空穴形態(tài)均被報道存在不足之處［2-5］。上述幾種空化模型的相變率表達式均基于單一球形空泡動力學(xué)中常用的Rayleigh-Plesset（R-P）方程的線性形式，該方程作了過多假設(shè)，如忽略空泡二階運動項、液相表面張力等，使得該類空化模型還不能完全且較好地勝任軸流泵葉輪內(nèi)空化流動問題的研究。

為了改善這種基于線性R-P方程的空化模型模擬空化流動的能力，本文提出了一種非線性空泡模型，聯(lián)立兩相-3種組分的均相流模型，得到了一種非線性空化模型。運用用戶自定義語言，把該模型嵌入到ANSYS FLUENT 14.5平臺，并對某一型號的軸流泵葉輪內(nèi)空化流動進行數(shù)值計算。通過與傳統(tǒng)的Schnerr-Sauer模型的數(shù)值計算結(jié)果及實驗結(jié)果作對比，驗證非線性空化模型在預(yù)測軸流泵空化流動中的可行性與準確性。

1 數(shù)值計算模型

1.1 非線性空泡模型

非線性空泡模型基于廣義的球形空泡動力學(xué)R-P方程建立。液態(tài)水及空泡均屬于低黏度工質(zhì)，其運動黏度系數(shù)均為 10-6數(shù)量級，Ivany等［6］和Knapp等［7］研究單一空泡的增長與潰滅過程發(fā)現(xiàn)，水的粘性引起的阻尼效應(yīng)可以忽略不計。同時，Schnerr等［8］在研究繞三維 Twist-11扭曲水翼空化流場時發(fā)現(xiàn)，云空泡脫落、潰滅及隨后的沖擊脈動動力學(xué)特征均是由空泡的慣性效應(yīng)決定，黏性的阻尼作用十分微小。此外，本文作者在前期工作［9］中還證明了表面張力系數(shù)對空泡半徑增長的影響明顯高于粘性力的影響。因此，非線性空泡模型在線性R-P方程基礎(chǔ)上還考慮了空泡半徑二階導(dǎo)數(shù)及液相表面張力的影響，線性、非線性R-P方程分別為：

式中pv（t）——空泡內(nèi)壓力；

p∞（t）——無窮遠處壓力；

t——時間；

ρl——液相密度；

R——空泡半徑；

S——表面張力系數(shù)。

式（2）中，等號左邊為壓力驅(qū)動項，等號右邊依次為二階導(dǎo)數(shù)項、一階導(dǎo)數(shù)項及表面張力項。

為了對比2種不同空泡動力學(xué)模型對空泡半徑變化的差異，利用MATLAB軟件，采用4-5階RUNGE-KUTTA算法分別對式（1）和（2）進行求解，采用的初始條件如下：空泡初始半徑R0=4.0 μm、液相表面張力系數(shù)S=0.073 N/m、驅(qū)動項pv（t）-p∞（t）=105Pa。計算可以得到該空化核流過低壓區(qū)時R（t）一階導(dǎo)數(shù)的一個典型解如圖1所示，從圖中可以看到，線性R-P方程預(yù)測的空泡半徑增長率為一恒定值，且當空泡半徑小于30 μm時，非線性R-P方程預(yù)測的空泡半徑增長率大于線性R-P方程計算的結(jié)果。由此可見，線性R-P方程在空泡初生階段（R＜30 μm）會過小預(yù)測空泡半徑的增長速度。

圖1 根據(jù)線性、非線性R-P方程計算得到的空泡半徑增長率

基于兩相-3種組分的均相流模型及非線性空泡模型，建立了一種非線性R-P空化模型，該空化模型考慮了不可凝結(jié)氣體、湍流脈動壓力對空化初生的影響。相變率推導(dǎo)詳細過程，以及模型系數(shù)無關(guān)性分析參見文獻［9］，相變率基本方程組為：

式中Se——蒸發(fā)率；

N0——空化核密度，m-3，取N0=1×1013m-3［8］；

ρg——NCG密度，按照密度形式理想氣體狀態(tài)方程計算；

R0——空泡初始半徑，μm，R0=1 μm［10］；

Sc——凝結(jié)率；

fg——NCG質(zhì)量分數(shù)，取1.5×10-6［10］；

psat——飽和蒸汽壓，Pa，取psat=3 169 Pa；

k——湍動能。

1.2 湍流模型

傳統(tǒng)的雙方程湍流模型會過度預(yù)測閉合空穴尾部的湍流粘度，使得空穴形態(tài)預(yù)測失真，而在大型復(fù)雜三維流動問題中應(yīng)用直接數(shù)值模擬（DNS）或者大渦模擬（LES）時很難找到一個網(wǎng)格無關(guān)解，且消耗大量的計算資源。因此，為了保證計算精度及節(jié)省計算成本，本文采用一種修正的濾波器湍流模型（MFBM）［11］。該模型通過一個函數(shù)fMFBM聯(lián)立原始濾波器模型（FBM）和密度修正模型（DCM）來共同修正傳統(tǒng)雙方程湍流模型模擬空化兩相流時的不足，修正后的模型兼?zhèn)銯BM和DCM的優(yōu)點。MFBM的湍流黏度系數(shù)為：

式中，n=10、Cμ=0.09、C3=1.0，λ為濾波尺。

2 計算域及網(wǎng)格劃分

研究對象為一比轉(zhuǎn)速ns=692的軸流泵，設(shè)計流量為Q=312 m3/h、設(shè)計揚程H=2.95 m、轉(zhuǎn)速n=1 450 r/min，其過流部件主要包括進出口管路、葉輪和導(dǎo)葉，三維造型過程中保留了葉頂間隙處水體，其厚度為0.5 mm。過流部件如圖2所示，整個計算域網(wǎng)格劃分在ICEM CFD軟件完成，采用結(jié)構(gòu)化網(wǎng)格。

圖2 模型軸流泵計算域

為了盡量消除網(wǎng)格數(shù)量造成的求解誤差，分別對計算域劃分了6種不同數(shù)量的網(wǎng)格作網(wǎng)格獨立性分析，網(wǎng)格節(jié)點數(shù)分別為N1=1 170 764、N2=1 756 146、N3=2 341 528、N4=2 926 910、N5=3 212 292和N6=3 697 674。在最優(yōu)工況下（Q=302 m3/h）［12］，計算基于不同網(wǎng)格模型泵的揚程來確定最終網(wǎng)格，計算結(jié)果如圖3所示，從圖中可以看到，當網(wǎng)格節(jié)點數(shù)大于N5時，揚程的相對變化在1%以內(nèi)，表明繼續(xù)加密網(wǎng)格對數(shù)值計算結(jié)果的影響可以忽略。N5網(wǎng)格計算得到的揚程為3.20 m，對應(yīng)的試驗揚程為3.04 m，二者相對誤差為5.26%，計算精度較高。因此，采用網(wǎng)格節(jié)點數(shù)為N5的網(wǎng)格作為本文最終的計算網(wǎng)格，其過流部件網(wǎng)格總成及葉輪網(wǎng)格示意如圖4所示。數(shù)值計算時，以非空化流場的計算結(jié)果作為空化非定常數(shù)值計算的初始值，在非定常數(shù)值計算中設(shè)置時間步長為Δt=0.11 ms。采用托馬（空化）系數(shù)來判斷空化程度，表達式為：

式中NPSHA——裝置空化余量；

H——揚程。

圖3 不同數(shù)量網(wǎng)格下計算揚程變化

圖4 網(wǎng)格示意

3 數(shù)值方法驗證

該模型軸流泵空化特性測試在江蘇大學(xué)國家水泵及系統(tǒng)工程技術(shù)研究中心完成，按照回轉(zhuǎn)動力泵水力性能驗收試驗標準GB/T3216-2005進行模型泵性能測試，工作介質(zhì)為常溫清水，試驗裝置簡圖如圖5所示。

圖5 軸流泵空化測試系統(tǒng)示意

圖6示出0.8QBEP，1.0QBEP和1.2QBEP3種工況下不同空化模型計算得到的空化性能曲線和臨界空化余量（揚程下降3%），與相應(yīng)試驗測量值的對比，其中計算揚程均為時均化后的結(jié)果，計算總時間為0.33 s。

圖6 空化特性預(yù)測值與試驗測量值的對比

從圖中可以看到，不同工況下，當空化數(shù)足夠大時，泵內(nèi)沒有發(fā)生空化，不同空化模型預(yù)測到的揚程十分接近。隨著空化數(shù)的減小，當空化發(fā)生時，各工況下?lián)P程的測量值相比于計算值提前發(fā)生斷裂，說明數(shù)值計算會一定程度上低估泵內(nèi)空化的初生及發(fā)展過程，但非線性模型預(yù)測的揚程斷裂時的空化數(shù)明顯大于Schnerr-Sauer模型的計算結(jié)果，表明非線性模型會更早預(yù)測到空化的發(fā)生，造成這種現(xiàn)象的原因有兩點：（1）非線性模型中考慮了NCG的影響，該部分氣體充當“空化核”的作用，加速了空化的發(fā)生；（2）非線性模型中考慮了湍流脈動對空化臨界壓力的影響。從圖6（d）中可以看到，非線性模型預(yù)測到的臨界空化余量與試驗測量值吻合度更高。為了進一步對比兩種空化模型捕捉空化特征的能力，截取了空化充分發(fā)展階段（揚程下降10%）不同空化模型預(yù)測的葉輪流道空泡形態(tài)，鄰近比轉(zhuǎn)速（ns=728）相同空化數(shù)條件下某一軸流泵的可視化試驗結(jié)果作為數(shù)值計算結(jié)果的對比如圖7所示?？梢园l(fā)現(xiàn)，不同空化模型的預(yù)測結(jié)果均捕捉到了尾緣處的小體積空泡團向葉輪流道脫落，此時實驗結(jié)果中葉頂間隙內(nèi)附著空穴與延伸至葉輪流道內(nèi)的空泡混合形成的總空泡團外形上呈三角形，非線性模型預(yù)測到的葉頂間隙附近處空泡團的三角形外形明顯與試驗結(jié)果吻合更好，說明非線性模型預(yù)測的葉頂間隙內(nèi)空化發(fā)展更加充分。綜上說明，非線性空化模型相比于傳統(tǒng)的Schnerr-Sauer模型預(yù)測軸流泵空化性能及對空化細節(jié)捕捉的能力更強。

圖7 空化形態(tài)與試驗結(jié)果的對比

4 葉輪內(nèi)2種主要空化特征

4.1 葉頂間隙空化

軸流泵葉輪內(nèi)空化主要包含葉頂間隙空化和葉片吸力面空化兩大類，張德勝等［13］參考軸流式風(fēng)機葉輪內(nèi)流動特點，將軸流泵葉頂間隙處流動進行了歸類和描述，認為葉頂間隙空化包括間隙內(nèi)附著空化、泄漏渦渦帶空化、剪切層分離渦空化及尾緣脫離渦空化。鑒于非線性模型對空化捕捉能力優(yōu)于傳統(tǒng)Schnerr-Sauer模型，截取了1.0QBEP，空化數(shù)σ=0.59工況下非線性模型計算得到的葉輪內(nèi)相對速度流線分布，以及葉頂間隙附近處空泡分布分別如圖8，9所示。從圖中可以看到，非線性空化模型基本上捕捉到了葉頂間隙附近上述4種不同類型的空化。

圖8 葉輪內(nèi)相對速度分布（1.0 QBEP，σ=0.59）

圖9 葉頂間隙附近空化分類

圖10示出1.0QBEP工況不同空化數(shù)下，非線性模型計算得到的葉頂間隙附近處的空泡形態(tài)分布，為了使計算結(jié)果精確度更高，該圖中的所有的空泡分布均采用非定常求解器計算后處理得到，每一工況均運行了2 000迭代步，即每隔t=0.22 s的時間處理一張分圖?？梢钥闯?，當空化數(shù)降低至σ=2.09時，葉頂間隙內(nèi)初生空化形成；當空化數(shù)大幅降低至1.41時，間隙內(nèi)附著型片狀空穴長度出現(xiàn)了較為明顯的增長，如圖10（b）所示。隨著空化數(shù)緩慢降低至σ=1.26時，間隙內(nèi)附著型片狀空穴的長度出現(xiàn)了顯著的增長；隨后，當空化數(shù)降低至σ=0.94時，發(fā)現(xiàn)附著空穴的長度出現(xiàn)了迅速增長。由此可見，葉頂間隙附著型空化生長速度先極為緩慢，當空化數(shù)降低至一定程度時出現(xiàn)快速增長。當空化數(shù)低至σ=0.59時，此時葉頂間隙附近的空化由葉頂泄漏渦渦帶空化、剪切層分離渦空化以及葉頂翼型脫離渦空化構(gòu)成。同時還說明，葉頂泄漏渦渦帶空化、剪切層分離渦空化及葉頂翼型尾緣脫離渦空化的發(fā)生明顯滯后于葉頂間隙內(nèi)的附著型空化。

圖10 非線性模型計算得到的不同工況下葉頂間隙附近空穴形態(tài)（Δ t=0.22 s，Q=1.0 QBEP）

圖11示出1.0QBEP工況不同空化數(shù)下，非線性模型計算得到的50%弦長葉頂間隙內(nèi)的空泡體積分數(shù)及速度矢量分布。從圖中可以看到，當空化數(shù)σ=1.05時，在葉片頂端，由于葉片和泵腔內(nèi)壁面存在間隙，流體的粘性力不足以克服葉頂翼型工作面和吸力面之間的壓差，葉片壓力面后的流體收縮被迫加速流入葉頂間隙，當流體通過葉頂間隙區(qū)域后以射流的形式逆向流入主流，并與葉片吸力面葉頂附近區(qū)域的主流相互作用形成逆時針的泄漏渦，此時因泄漏渦渦心處的壓力值高于該溫度下的飽和蒸汽壓而沒有發(fā)生泄漏渦渦帶空化。當空化數(shù)降低至σ=0.91時，葉頂處吸力面發(fā)生了空化。隨著空化數(shù)降低至σ=0.85時，在射流區(qū)域內(nèi)形成了射流剪切層空化，而且此時吸力面上附著空化厚度相比于前一工況出現(xiàn)了急劇增加。當空化數(shù)降低至σ=0.78時，泄漏渦渦心處發(fā)生空化，泄漏渦渦帶空化范圍隨著空化數(shù)的降低而不斷擴大，其內(nèi)部所包含的空泡的體積分數(shù)也逐漸增大。在主流的作用下泄漏渦渦帶逐漸向吸力面靠近并與葉頂處吸力面空化發(fā)生摻混。葉頂泄漏渦是一種高度的三維流動結(jié)構(gòu)，葉頂泄漏渦與渦帶空化互相作用。隨著空化數(shù)的降低，渦帶空化呈現(xiàn)出不同的形狀，其軸向位置逐漸向吸力面靠近，泄漏渦渦心處的最低壓力大小決定了渦帶空化區(qū)空泡相的體積分。

圖11 非線性模型計算得到的不同工況50%弦長處葉頂間隙內(nèi)流動（Δt=0.22 s，Q=1.0 QBEP）

4.2 葉片吸力面空化

圖12示出1.0QBEP工況不同空化數(shù)下，非線性模型計算得到的葉片吸力面上空泡體積分數(shù)分布。從圖中可以看到，當空化數(shù)為σ=1.47時，初生空化出現(xiàn)在葉片吸力面進口靠近輪緣處。對比圖10（a）所對應(yīng)的空化數(shù)可知，軸流泵葉片吸力面空化明顯滯后于葉頂間隙內(nèi)空化。隨著空化數(shù)的降低，空化區(qū)逐漸向葉片出口及輪轂方向擴展，靠近輪緣處區(qū)域內(nèi)空泡相體積分數(shù)急劇增加。當空化數(shù)降低至σ=0.85時，吸力面空化已經(jīng)發(fā)展至輪轂，并在2/3倍弦長處形成了一條與葉輪軸向垂直的空泡帶狀區(qū)域，此時輪轂表面區(qū)域會因發(fā)生空化而附著有空泡團，即轂渦空化。在臨界空化工況點時，吸力面上空化帶狀區(qū)分別向葉片進出口方向延伸，并逐漸與輪緣處的空化區(qū)摻混一起，靠近輪轂處空化區(qū)空泡相體積分數(shù)相比前一工況急劇增加。隨著空化數(shù)繼續(xù)降低，轂渦空化程度加劇，摻混后的吸力面空穴末端逐漸向葉片出口延伸，如圖12（f）所示。當空化數(shù)降低至σ=0.63時，吸力面上2/3倍弦長處至葉片出口，在整個葉高范圍內(nèi)完全被高體積分數(shù)的空泡覆蓋。對比圖12（g）和（h）可以發(fā)現(xiàn)，當空化數(shù)降低至一定程度時，葉片出口靠近輪轂處發(fā)生了嚴重空化，空化區(qū)內(nèi)空泡體積分數(shù)極高，在整個吸力面空化區(qū)域的后部靠近輪轂的位置，空泡體積分數(shù)發(fā)生了急劇變化，這種劇烈變化是由于空泡運動至葉片出口附近的高壓區(qū)發(fā)生潰滅所造成?？梢酝茢?，當軸流泵葉輪內(nèi)發(fā)生嚴重空化時葉片出口位置且靠近輪轂處最容易發(fā)生空蝕破壞。

圖12 非線性模型計算得到的葉片吸力面空泡體積分數(shù)云圖（Δ t=0.22 s，Q=1.0 QBEP）

5 結(jié)論

（1）由于非線性空化模型考慮了NCG及湍流脈動的影響，使得該模型預(yù)測的軸流泵葉輪內(nèi)空化初生及發(fā)展過程早于傳統(tǒng)的Scherr-Sauer模型的計算結(jié)果，且非線性空化模型預(yù)測的空化性能曲線與試驗測量值吻合度更好，對空化細節(jié)捕捉的能力也更強。

（2）葉頂間隙附著型空化生長速度開始比較緩慢，當空化數(shù)降低至一定程度時出現(xiàn)快速增長；葉頂泄漏渦渦帶空化、剪切層分離渦空化及葉頂翼型尾緣脫離渦空化的發(fā)生明顯滯后于葉頂間隙內(nèi)的附著型空化。

（3）當軸流泵葉片吸力面發(fā)生嚴重空化時，空化區(qū)域的后部靠近輪轂的位置空泡體積分數(shù)發(fā)生了急劇變化，該處最容易發(fā)生空蝕破壞。