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

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

0 引言

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

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

1 數(shù)值計(jì)算模型

1.1 非線性空泡模型

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

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

p∞（t）——無(wú)窮遠(yuǎn)處壓力；

t——時(shí)間；

ρl——液相密度；

R——空泡半徑；

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

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

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

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

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

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

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

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

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

Sc——凝結(jié)率；

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

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

k——湍動(dòng)能。

1.2 湍流模型

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

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

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

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

圖2 模型軸流泵計(jì)算域

為了盡量消除網(wǎng)格數(shù)量造成的求解誤差，分別對(duì)計(jì)算域劃分了6種不同數(shù)量的網(wǎng)格作網(wǎng)格獨(dú)立性分析，網(wǎng)格節(jié)點(diǎn)數(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］，計(jì)算基于不同網(wǎng)格模型泵的揚(yáng)程來(lái)確定最終網(wǎng)格，計(jì)算結(jié)果如圖3所示，從圖中可以看到，當(dāng)網(wǎng)格節(jié)點(diǎn)數(shù)大于N5時(shí)，揚(yáng)程的相對(duì)變化在1%以內(nèi)，表明繼續(xù)加密網(wǎng)格對(duì)數(shù)值計(jì)算結(jié)果的影響可以忽略。N5網(wǎng)格計(jì)算得到的揚(yáng)程為3.20 m，對(duì)應(yīng)的試驗(yàn)揚(yáng)程為3.04 m，二者相對(duì)誤差為5.26%，計(jì)算精度較高。因此，采用網(wǎng)格節(jié)點(diǎn)數(shù)為N5的網(wǎng)格作為本文最終的計(jì)算網(wǎng)格，其過(guò)流部件網(wǎng)格總成及葉輪網(wǎng)格示意如圖4所示。數(shù)值計(jì)算時(shí)，以非空化流場(chǎng)的計(jì)算結(jié)果作為空化非定常數(shù)值計(jì)算的初始值，在非定常數(shù)值計(jì)算中設(shè)置時(shí)間步長(zhǎng)為Δt=0.11 ms。采用托馬（空化）系數(shù)來(lái)判斷空化程度，表達(dá)式為：

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

H——揚(yáng)程。

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

圖4 網(wǎng)格示意

3 數(shù)值方法驗(yàn)證

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

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

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

圖6 空化特性預(yù)測(cè)值與試驗(yàn)測(cè)量值的對(duì)比

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

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

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

4.1 葉頂間隙空化

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

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

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

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

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

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

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

4.2 葉片吸力面空化

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

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

5 結(jié)論

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

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

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