誘導(dǎo)輪非對稱空化流動數(shù)值計算和試驗研究

李龍賢，溫永鑫

（1.北京航天動力研究所；北京 100076；2.北京電子工程總體研究所；北京 100854）

0 引言

在航天推進(jìn)系統(tǒng)中，對作為中大推力液體火箭發(fā)動機(jī)關(guān)鍵增壓裝置的渦輪泵的工作穩(wěn)定性和可靠性等有嚴(yán)格的要求。誘導(dǎo)輪作為渦輪泵的前置增壓部件，對改善渦輪泵的抽吸性能、降低推進(jìn)劑貯箱內(nèi)部壓力及質(zhì)量具有重要作用。高性能誘導(dǎo)輪轉(zhuǎn)速較高，通常允許流場內(nèi)存在一定程度的空化，但空化的產(chǎn)生會使流場趨于不穩(wěn)定，甚至對發(fā)動機(jī)的工作穩(wěn)定性和可靠性產(chǎn)生影響。

誘導(dǎo)輪內(nèi)的空化存在多種形式，非對稱空化是已知的一種存在于渦輪泵揚(yáng)程下降過程中的復(fù)雜流動現(xiàn)象。對于附著型非對稱空化，誘導(dǎo)輪內(nèi)各通道的空化區(qū)長度表現(xiàn)出明顯的不同。在此工況下，須單獨(dú)考慮和分析每個葉片，但這將使得空化模型變得非常復(fù)雜。Horiguchi等[1]利用奇點(diǎn)方法給出了預(yù)測交替型空化的近似模型，為非對稱型空化的預(yù)測提供了重要的借鑒。本文基于奇點(diǎn)法理論，并結(jié)合空化波動理論對某低溫介質(zhì)三葉片誘導(dǎo)輪的非對稱空化特征進(jìn)行研究。

1 研究方法

1.1 研究對象

以某低溫介質(zhì)三葉片變螺距誘導(dǎo)輪為研究對象。為了計算及試驗方便，將離心輪、擴(kuò)壓器、誘導(dǎo)輪等組成渦輪泵整體，模型如圖1所示，誘導(dǎo)輪實物如圖2所示。誘導(dǎo)輪的實際工作介質(zhì)為低溫液氧，但在本文研究中由于采用低溫介質(zhì)進(jìn)行試驗難度和風(fēng)險較大，故采用密度和黏度相近的常溫水代替。

圖1 液氧渦輪泵幾何模型Fig.1 Liquid oxygen turbopump geometry model

圖2 誘導(dǎo)輪實物圖Fig.2 Physical view of inducer of liquid oxygen turbopump

1.2 控制方程

控制方程主要由連續(xù)方程、動量方程、能量方程構(gòu)成。其中連續(xù)方程為混合相質(zhì)量守恒方程，在連續(xù)方程中增加一個氣態(tài)相或液態(tài)相質(zhì)量守恒方程并帶入質(zhì)量輸運(yùn)源相，便可導(dǎo)入空化模型[2]。

（1）連續(xù)方程

式中：t為時間；ρ為氣液兩相流體密度，ρ=αlρl+αvρv；U為速度向量；?(ρU)為速度散度。αl為液體體積分?jǐn)?shù)；ρl為液體密度；αv為氣體體積分?jǐn)?shù)；ρv為氣體密度。

（2）動量方程

式中：U×U為向量積；τ為表面力，?τ為表面力τ的散度；p為流場壓力；SM為由體積力牽引的動量源項。

（3）能量方程

式中：e為工質(zhì)內(nèi)能；U為速度標(biāo)量值；T為工質(zhì)溫度；q?為工質(zhì)與外界的熱交換量；SΦ為耗散函數(shù)。

（4）空化模型及空泡動力學(xué)方程

空化模型為基于均質(zhì)多相質(zhì)量輸運(yùn)模型的空化模型（又稱1方程模型），并假設(shè)相間無滑移。數(shù)值計算過程中，包含多相流及空化的計算模型適用于常規(guī)的黏性流體力學(xué)計算方法，常規(guī)計算單相流的湍流模型對于多相流流場計算依然適用[3]。

式中：Re、Rc分別表示氣泡產(chǎn)生和氣泡潰滅的質(zhì)量輸運(yùn)源項。

空泡動力學(xué)方程采用Rayleigh-Plesset方程[4]，即：

式中：RB為空泡半徑；pv為飽和壓力；p∞為遠(yuǎn)場壓力；S為液體表面張力系數(shù)；νl為液體運(yùn)動黏度。

采用上述組合模型計算出的結(jié)果能較好地反映出空化區(qū)隨流場壓力的變化情況?；赗ayleigh-Plesset空泡動力學(xué)理論推導(dǎo)出的空化模型為均勻化模型，目前很難用質(zhì)量輸運(yùn)空化模型準(zhǔn)確預(yù)測附著型非對稱流動[5]流場內(nèi)空化區(qū)波動以及不對稱特性。因此，本文引入奇點(diǎn)法表示空化波動方程，即：

式中：Ue-iαz表示主流流動勢；N表示葉片數(shù)；wn為流場受第n個葉片的擾動，其物理意義是受源項及葉片上漩渦影響的流動勢，單位為J/g。

采用空化數(shù)σ定量描述空化的劇烈程度[6]，空化數(shù)為無量綱參數(shù)，其數(shù)學(xué)表達(dá)式為：

式中：p∞和u∞分別為液體流動參考壓力和參考速度，在葉輪機(jī)械流場中，p∞通常取入口總壓pin；u∞取葉尖圓周速度uT，因此空化數(shù)可以表示為[7]：

式中：Ω為角速度；rT為葉片圓周半徑。若不計比尺效應(yīng)，理論上如果空化數(shù)σ相同，空化程度也相同[8]。

1.3 試驗研究

圖3為誘導(dǎo)輪空化試驗及測試系統(tǒng)示意圖。系統(tǒng)主要由三個分系統(tǒng)組成：泵組合體外特性參數(shù)試驗系統(tǒng)、流場脈動壓力測量系統(tǒng)和空化流場可視化試驗系統(tǒng)。由于直接用低溫工質(zhì)進(jìn)行試驗的難度大，成本較高，試驗風(fēng)險較大，因此，采用與實際低溫工質(zhì)液態(tài)氧黏度和密度等物性參數(shù)相近的常溫水進(jìn)行試驗。

圖3 誘導(dǎo)輪空化試驗及測試系統(tǒng)Fig.3 Diagram of inducer cavitation test system

圖4為試驗系統(tǒng)現(xiàn)場圖。試驗過程中采用高速相機(jī)通過觀察窗口連續(xù)高速拍攝流場的瞬態(tài)圖像。

圖4 可視化試驗系統(tǒng)現(xiàn)場圖Fig.4 Visualization test system

誘導(dǎo)輪空化流場觀測窗口采用透明有機(jī)玻璃制成，圖5為觀測窗口實物圖。

圖5 誘導(dǎo)輪空化觀測窗口實物圖Fig.5 Inducer cavitation observation window

2 數(shù)值計算和試驗結(jié)果分析

2.1 數(shù)值計算結(jié)果分析

圖6為試驗和數(shù)值仿真計算結(jié)果的對比。

圖6 空化數(shù)-泵出口壓力系數(shù)曲線Fig.6 Cavitation number-pump outlet pressure coefficient curve

由圖6可以看出，由數(shù)值計算得出的泵出口壓力系數(shù)ψ隨空化數(shù)σ的變化關(guān)系曲線與試驗結(jié)果吻合較好，說明由數(shù)值計算得出的泵組合體外特性參數(shù)比較準(zhǔn)確，最大誤差在5%以內(nèi)。壓力系數(shù)：

式中：pout為泵出口壓力，壓力系數(shù)和空化數(shù)均為無量綱參數(shù)。

通過數(shù)值計算和試驗可以發(fā)現(xiàn)，本文研究的誘導(dǎo)輪，非對稱空化區(qū)開始出現(xiàn)的空化數(shù)σ≈0.03，空化數(shù)區(qū)間為0.01≤σ≤0.03，非對稱空化出現(xiàn)的空化數(shù)與斷裂揚(yáng)程發(fā)生時空化數(shù)很接近（斷裂揚(yáng)程指泵出口壓力系數(shù)急劇下降階段所對應(yīng)的泵揚(yáng)程）。

以空化數(shù)σ=0.022為例，介紹試驗過程中觀測到的非對稱空化區(qū)的特征。圖7為σ=0.022情況下誘導(dǎo)輪葉片和導(dǎo)葉葉片表面動態(tài)壓力分布云圖。可以看到，此時低壓區(qū)已覆蓋葉片前緣部分，壓力最低的區(qū)域分布在吸力面修圓的表面和葉尖處。

圖7 葉片表面壓力動態(tài)分布試驗云圖Fig.7 Test image of pressure dynamic distribution on blade surface

圖8為σ=0.022工況下誘導(dǎo)輪葉片通道內(nèi)空化區(qū)隨時間變化的動態(tài)云圖，（a）中的0°（0 s）為初始角度（初始時刻）流場試驗云圖，（b）（c）（d）依次為誘導(dǎo)輪轉(zhuǎn)過90°、180°、270°時在同一觀測位置觀察到的誘導(dǎo)輪葉片流場云圖?？梢钥闯觯藭r空化區(qū)主要分布在葉片前緣附近吸力面表面及葉片修圓吸力面上，葉片通道間的空化區(qū)表現(xiàn)出明顯的非對稱性，同時每個葉片上的附著空化區(qū)也在不斷地伸縮變化。

圖8 葉片通道內(nèi)空化區(qū)動態(tài)分布試驗云圖Fig.8 Cavitation area dynamic distribution on stator and rotor blade surface

圖9為采用數(shù)值計算得到的葉片通道內(nèi)空化區(qū) 的動態(tài)分布圖。除了附著在葉片修圓段外緣以外的空化區(qū)，靠近葉根的空化區(qū)以云狀氣團(tuán)的形式 分布在葉片通道入口，對通道造成阻塞。

圖9 葉片通道內(nèi)非對稱空化區(qū)數(shù)值計算云圖Fig.9 Nonsymmetrical cavitation area in flow passge obtained by numerical simulation

非對稱空化出現(xiàn)時伴隨著大量的空化區(qū)氣泡的脫落和破碎，按照空化的類別，此時的空化屬于片空化和云空化的混合態(tài)[9]。云空化是由片空化發(fā)展到一定程度氣泡在下游脫落破碎產(chǎn)生的，是由大量微泡群和液態(tài)工質(zhì)混合而成的兩相流。云空化的形成與片空化的尾部回射流以及片空化的整體失穩(wěn)密切相關(guān)。云空化的出現(xiàn)總是與流動失穩(wěn)、旋渦和湍流脈動等復(fù)雜流動現(xiàn)象聯(lián)系在一起，其發(fā)生、發(fā)展和潰滅會帶來結(jié)構(gòu)的震動、噪聲和空蝕[10]。此時泵在外特性上表現(xiàn)為揚(yáng)程降低3%，工程上認(rèn)為泵內(nèi)開始發(fā)生空蝕。云空化的出現(xiàn)從微觀的角度解釋了空蝕出現(xiàn)的原因。

圖10為空化數(shù)σ=0.022時計算得到的誘導(dǎo)輪葉片修圓末端橫截面上的壓力動態(tài)分布?？梢钥闯?，此時截面大部分被低壓區(qū)覆蓋，壓力較高的區(qū)域只出現(xiàn)在葉片壓力面葉尖處及葉片表面小部分區(qū)域，壓力的等壓線隨時間變化非常明顯。

圖10 葉片修圓末端橫截面壓力動態(tài)分布數(shù)值計算云圖Fig.10 Pressure dynamic distribution on blade tip relief cross section by numerical calculation

圖11為σ=0.022時計算得到的該截面空化區(qū)隨時間的動態(tài)分布?？梢钥吹剑栈瘏^(qū)主要分布在吸力面葉尖及葉中表面上。從圖11可以清晰地看到空化區(qū)在三個葉片的非對稱性分布。為便于分析，將誘導(dǎo)輪葉片編號，圓心處的箭頭方向為誘導(dǎo)輪旋轉(zhuǎn)方向，附著空化區(qū)面積最大的為葉片1，沿著與誘導(dǎo)輪旋轉(zhuǎn)方向相反的方向依次為葉片2和葉片3；葉片1、2、3下游的葉片通道分別為通道1、通道2和通道3。葉片1表面的空化區(qū)體積分?jǐn)?shù)明顯高于葉片2和葉片3表面；從圖11可以較明顯地看到空化區(qū)氣泡脫落過程?？栈瘏^(qū)氣泡脫落發(fā)生在葉片吸力面，部分空化區(qū)氣泡脫落并逐漸消失（圖11紅色虛線圈內(nèi)為空化區(qū)變化過程）?？栈瘏^(qū)氣泡脫落是云空化形成的誘導(dǎo)因素。

圖11 葉片修圓末端橫截面空化區(qū)動態(tài)分布數(shù)值計算云圖Fig.11 Cavitation area dynamic distribution on blade tip relief cross section by numerical calculation

圖12為空化數(shù)σ=0.022時誘導(dǎo)輪三個葉片表面平均空化區(qū)體積分?jǐn)?shù)?？梢钥吹剑藭r葉片上的空化區(qū)以大致相同的周期隨時間不斷波動，空化區(qū)波動周期約為葉輪旋轉(zhuǎn)周期的3倍；葉片3表面的空化體積分?jǐn)?shù)平均值比葉片2高出約25%，三個葉片上的空化區(qū)呈現(xiàn)出明顯的非對稱性。

圖12 誘導(dǎo)輪葉片表面空化區(qū)體積分?jǐn)?shù)波動曲線Fig.12 Fluctuation curve of inducer blade surface cavitation area volume fraction（σ=0.022）

2.2 試驗結(jié)果分析

圖13～15分別為誘導(dǎo)輪3個葉片依次經(jīng)過觀察窗口的連續(xù)成像，試驗工況為轉(zhuǎn)速4 000 r/min，空化數(shù)σ=0.022，流量系數(shù)φ=0.059，在此階段發(fā)現(xiàn)了比較強(qiáng)烈的非對稱空化，通道1的空化區(qū)比其他兩個要劇烈得多。Tsuijimoto[9]也在試驗中發(fā)現(xiàn)了這種非對稱空化，并證實這種空化不穩(wěn)定形式出現(xiàn)的空化數(shù)要小于旋轉(zhuǎn)空化。這種形式的非對稱流動流道中的空化區(qū)較穩(wěn)定，空化區(qū)隨葉片旋轉(zhuǎn)而同步旋轉(zhuǎn)，這與旋轉(zhuǎn)空化截然不同，將這種非對稱空化稱之為附著型非對稱空化。這種流道中的非對稱空化區(qū)為云空化的一種，誘導(dǎo)輪中云空化的出現(xiàn)是揚(yáng)程開始顯著下降和葉片空化破壞的誘導(dǎo)因素，云空化出現(xiàn)時對應(yīng)的空化數(shù)與斷裂揚(yáng)程對應(yīng)的空化數(shù)非常接近。

圖13 通道1空化區(qū)連續(xù)成像Fig.13 Continuous imaging for the first blade channel cavitation area

圖16為從試驗結(jié)果中提取的空化區(qū)光學(xué)圖像，從該組圖像中可以更清楚地看到誘導(dǎo)輪3個葉片通道中的空化區(qū)差異。

圖14 通道2空化區(qū)連續(xù)成像Fig.14 Continuous imaging for the second blade channel cavitation area

圖16 空化區(qū)光學(xué)圖像Fig.16 Cavitatoin area optical image

圖15 通道3空化區(qū)連續(xù)成像Fig.15 Continuous imaging for the third blade channel cavitation area

圖17為經(jīng)過圖像處理后得到的空化區(qū)等高線圖像。

圖17 不同通道的空化區(qū)等高線圖像Fig.17 Cavitation area contour line image

對圖17的圖像進(jìn)行光譜亮度分析，得到空化區(qū)頻譜能量概率統(tǒng)計圖，如圖18所示，這是對空化圖像特征的定量分析[11]?？梢钥吹剑ǖ?的頻譜能量峰值均明顯大于其他兩個通道，表明通道1內(nèi)的空泡數(shù)量、空泡密度均大于另外兩個通道，通道2和通道3內(nèi)的空化區(qū)分布規(guī)律比較接近。

圖18 頻譜能量概率統(tǒng)計Fig.18 Spectrum energy probability statistics

利用圖像的灰度直方圖提取諸如均值、方差、能量及熵等來描述圖像特征[13]。設(shè)（x，y）為圖像中的一點(diǎn)，該點(diǎn)和離它只有微小距離的點(diǎn)（x+Δx，y+Δy）的灰度差值gΔ為：

式中：設(shè)灰度差值的所有可能取值有m級，點(diǎn)（x，y）在整個圖像上移動，統(tǒng)計出gΔ取各個數(shù)值的次數(shù)，由此做出gΔ的直方圖。由直方圖可以得到gΔ取值的概率p（i），p（i）表示第i級灰度出現(xiàn)的概率，相關(guān)的特征參數(shù)有[14]：灰度平均值灰度對比度，其中灰度平均值表示圖像的亮度，平均值越大，圖像平均亮度越亮，圖像中的空泡越多。對比度表示圖像的亮度差值，對比度越大，圖像中的空泡越密集，熵表示圖像的平整程度，熵越大，空化區(qū)分布越均勻。

對表1空化區(qū)圖像特征參數(shù)進(jìn)行分析，可以看到，通道1的灰度平均值、對比度及熵均明顯大于其他兩個通道的空化區(qū)，表明通道1內(nèi)的空泡數(shù)量、空泡密度均大于另外兩個通道。通道2和通道3內(nèi)的空化區(qū)分布規(guī)律非常接近。

表1 空化圖像特征分析（σ=0.022）Tab.1 Characteristic analysis for cavitation image（σ=0.022）

圖19為經(jīng)過平均化處理的不同葉片背面上的附著空化區(qū)長度及不同通道的流量隨空化數(shù)的變化關(guān)系。

圖19 通道空化區(qū)長度及流量的歸一化處理Fig.19 Averaging treatment for cavitation area length and flux in different channel

圖19曲線上標(biāo)注的數(shù)字對應(yīng)葉片編號，lc表示空化區(qū)長度，lm表示空化混合區(qū)長度，Q表示葉片流道內(nèi)流量?？梢钥吹?，葉片上出現(xiàn)不對稱附著空化區(qū)的空化數(shù)范圍為0.02≤σ≤0.035，非對稱空化出現(xiàn)時通道內(nèi)的流量也隨之不對稱，空化最嚴(yán)重的通道1流量最小，旋轉(zhuǎn)阻塞只是出現(xiàn)在空化程度最嚴(yán)重的通道1中[5]。隨著空化數(shù)繼續(xù)降低，空化區(qū)的不對稱性也逐漸降低，當(dāng)空化數(shù)進(jìn)一步降低，誘導(dǎo)輪揚(yáng)程接近斷裂揚(yáng)程，此時空化區(qū)已蔓延至葉片后緣，通道內(nèi)的空化區(qū)不對稱性消失。

3 結(jié)論

誘導(dǎo)輪內(nèi)的空化存在多種形式，非對稱空化是已知的存在于渦輪泵揚(yáng)程下降過程中的復(fù)雜流動現(xiàn)象的一種。本文采用數(shù)值計算和試驗方法對某低溫介質(zhì)誘導(dǎo)輪流場中的非對稱空化進(jìn)行了研究，得出以下結(jié)論：

（1）數(shù)值計算和可視化試驗證明了非對稱空化的存在，數(shù)值計算結(jié)果與可視化試驗結(jié)果一致性較好，說明本文中所采用的奇點(diǎn)法空化波動方程數(shù)值計算方法預(yù)測的非對稱空化具有較好的準(zhǔn)確性；

（2）揚(yáng)程明顯下降時葉片通道出現(xiàn)明顯的非對稱性，對應(yīng)誘導(dǎo)輪的空化數(shù)為σ=0.03；

（3）三葉片誘導(dǎo)輪內(nèi)發(fā)生非對稱空化時，其中一個通道內(nèi)的空化區(qū)明顯大于其他兩個通道，且每個通道的空化區(qū)大小以一個特定值為中心隨時間成周期性波動；

（4）非對稱空化出現(xiàn)后，通道內(nèi)的流量也隨之變得不對稱，空化最嚴(yán)重的通道對應(yīng)的流量最小。旋轉(zhuǎn)阻塞只是出現(xiàn)在空化程度最嚴(yán)重的通道中；

（5）非對稱空化的出現(xiàn)只發(fā)生在特定空化區(qū)域，對于本文研究的誘導(dǎo)輪0.02≤σ≤0.035?？栈瘮?shù)進(jìn)一步降低，誘導(dǎo)輪揚(yáng)程接近斷裂揚(yáng)程，此時空化區(qū)已蔓延至葉片后緣，葉片通道內(nèi)的空化區(qū)不對稱性消失。