液體火箭發(fā)動機誘導輪非定?？栈鲃犹匦詳?shù)值計算

康亞卓，趙興安，吳欽，黃彪，王國玉

(北京理工大學機械與車輛學院，北京 100081)

隨著航天技術(shù)的進一步發(fā)展，要求運載火箭具備更大的負載能力，而渦輪泵是液體火箭發(fā)動機系統(tǒng)的重要組件之一，其功能是由渦輪帶動離心泵高速運轉(zhuǎn)以提高火箭推進劑的壓力，然后再將推進劑送入火箭發(fā)動機的推力室進行燃燒，從而產(chǎn)生巨大的推力，為火箭飛行提供所需的動力[1].近年來，在高速運轉(zhuǎn)的火箭發(fā)動機渦輪泵前加裝誘導輪已成為保證渦輪泵獲取優(yōu)越空化性能的關(guān)鍵技術(shù).誘導輪的主要作用是對進入渦輪泵的推進劑進行加壓產(chǎn)生一定揚程以提高渦輪泵葉輪入口壓力，避免離心泵發(fā)生空化，從而提高渦輪泵空化性能并降低推進劑貯箱壓力和重量.然而在誘導輪的高速運轉(zhuǎn)過程中，其運行工況常常會發(fā)生變化，尤其是當其內(nèi)部局部壓力降低至飽和蒸氣壓以下發(fā)生空化時，將導致誘導輪內(nèi)部強烈的空化現(xiàn)象，引起性能下降、系統(tǒng)結(jié)構(gòu)振動加強，從而影響火箭發(fā)動機系統(tǒng)可靠性和渦輪泵誘導輪使用壽命[2-4].因此，渦輪泵誘導輪空化流動特性研究已成為現(xiàn)代高能液體火箭發(fā)動機研制中的關(guān)鍵技術(shù)問題之一，受到各航天大國的廣泛重視.

目前，國內(nèi)外關(guān)于誘導輪空化研究主要以試驗和數(shù)值計算研究為主.COUTIER等[5]采用壓力邊界條件、基于滑移網(wǎng)格技術(shù)對一單級單吸蝸殼離心泵進行了三維空化非定常相變流動的數(shù)值計算,并對渦輪泵誘導輪空化流動進行了試驗研究，分析了空化發(fā)展過程中空泡在誘導輪內(nèi)部的空間分布和發(fā)展規(guī)律.李森虎等[6]通過試驗研究了離心泵內(nèi)部空化流動，分析了有無誘導輪2種情況下的離心泵及誘導輪內(nèi)部流場分布，表明加誘導輪之后有效減弱了離心泵內(nèi)部空化.TANI等[7]探究了入口流量和流場結(jié)構(gòu)對誘導輪空化不穩(wěn)定性的影響，表明逆壓梯度造成的回流是引起誘導輪非定常特性的主要因素，當流量增加時，回流減弱，但是其與空泡潰滅之間的交互作用增強.劉春哲等[8]對單吸蝸殼離心泵進行了非定常流動的數(shù)值計算，得到了葉片旋轉(zhuǎn)空化的周期.

文中采用基于旋轉(zhuǎn)曲率修正的湍流模型對誘導輪空化特性進行數(shù)值計算，并通過試驗進行驗證，分析誘導輪非定?？栈鲃犹匦砸约傲鲌鰤毫γ}動特性.

1 數(shù)值計算方法

1.1 基本控制方程

采用均相流模型[9]進行空化流動計算，模型基本控制方程為

(1)

(2)

其中，

μm=μvαv+μlαl，

(3)

ρm=ρvαv+ρlαl，

(4)

以上式中：U為流場速度；p為壓力；μ，μt和μm分別為流體的黏性系數(shù)、湍流黏性系數(shù)和混合相動力黏性系數(shù)；ρv，ρl和ρm分別為氣相密度、液相密度和混合相密度；αv，αl分別為氣相和液相體積分數(shù)；μv，μl分別為氣相和液相黏性系數(shù).

1.2 考慮旋轉(zhuǎn)修正的湍流模型

為了更好地模擬流場中的旋轉(zhuǎn)效應(yīng)，采用標準k-ε模型[10-11]，并對其中湍動能生成項進行旋轉(zhuǎn)曲率修正.

標準k-ε模型的控制方程為

(5)

(6)

式中：Cε1=1.44，Cε2=1.92，σε=1.3，σk=1.0.

采用LAUNDER等[12-13]提出的滿足伽利略不變性修正方法，對湍動能生成項Pt進行旋轉(zhuǎn)-曲率修正，乘以修正系數(shù)fr，即

Pt→Ptfr，

(7)

(8)

(9)

(10)

以上式中：常數(shù)Cr1，Cr2和Cr3取值分別為1.0，2.0，1.0；Sij和ωij分別為應(yīng)變S和渦量ω的張量分量，S2=2SijSij，ω2=ωijωij，參數(shù)D2=max(S2,0.09).

1.3 空化模型

采用Zwart空化模型[14]模擬非定?？栈匦砸约皦毫γ}動，該模型的蒸發(fā)項和凝結(jié)項分別為

(11)

(12)

式中：αnuc為不溶性氣體體積分數(shù)；RB為空泡半徑；p和pv分別表示當?shù)貕簭姾蜌饣瘔簭?；Cprod和Cdest分別為凝結(jié)系數(shù)和蒸發(fā)系數(shù).相關(guān)經(jīng)驗系數(shù)在計算中取值分別為αnuc=5.0×10-4，RB=1.0×10-6m，Cprod=0.01，Cdest=50.

1.4 數(shù)值計算模型

1.4.1 誘導輪幾何建模

選取高揚程渦輪泵誘導輪DAPAMITO4[15-16]為研究對象，該誘導輪由意大利比薩Alta S.P.A.航天航空公司基于一種降階模型設(shè)計而成.DAPAMITO4誘導輪的主要設(shè)計參數(shù)分別為轉(zhuǎn)速Ω=3 000 r/min，流量系數(shù)Φ=0.070，葉片數(shù)Z=4，輪緣半徑rT=81.0 mm，進口葉頂安放角γTle=81.10°，進口輪轂半徑(葉片充分發(fā)展段)rHle=48.0 mm，出口輪轂半徑rHte=58.5 mm，葉片平均高度hm=27.75 mm，軸向長度(葉片充分發(fā)展段)ca=63.5 mm，進口輪轂半徑rH1=35.0 mm，軸向長度L=90.0 mm.定義流量系數(shù)Φ為

(13)

計算流體域包括進口段、誘導輪和出口段3部分，其中進口段和出口段分別為5倍以及8倍誘導輪輪緣直徑的管道.計算時，設(shè)置誘導輪的葉頂間隙為1 mm，總壓進口、質(zhì)量流量出口、流域所有壁面均采用無滑移壁面.動靜域之間以及葉頂間隙連接面采用Interface連接，進口段-誘導輪與誘導輪-出口段交界面類型采用Rotor Stator界面.建模后的計算域、邊界條件以及部分尺寸如圖1所示.

圖1 流體計算域以及邊界條件

Fig.1 Computational fluid domains and boundary conditions

1.4.2 網(wǎng)格劃分

整個誘導輪計算域均采用六面體結(jié)構(gòu)化網(wǎng)格劃分，如圖2所示，其中進口段與出口段網(wǎng)格利用ICEM CFD軟件進行劃分，網(wǎng)格數(shù)分別約為34萬、20萬，網(wǎng)格質(zhì)量均在0.7以上；誘導輪區(qū)域網(wǎng)格利用Turbo Grid軟件進行劃分，并對誘導輪葉片以及輪轂等壁面處網(wǎng)格進行加密處理，網(wǎng)格數(shù)約為100萬.

Fig.2 Meshes in inlet, outlet and inducer fluid domains

2 結(jié)果與討論

2.1 試驗驗證

定義空化數(shù)σ和揚程系數(shù)Ψ分別為

(14)

(15)

式中：pin為誘導輪進口壓力；pv為水的飽和蒸氣壓；ρ為水的密度；Δp為誘導輪出口壓增.

為驗證數(shù)值計算方法的正確性，對DAPAMIT O4誘導輪在轉(zhuǎn)速Ω=3 000 r/min、流量系數(shù)Φ=0.053工況下的空化特性進行計算，并與試驗數(shù)據(jù)[15]進行對比，如圖3所示.

圖3 誘導輪揚程系數(shù)曲線

由圖3可以看出：誘導輪揚程系數(shù)數(shù)值計算結(jié)果與試驗吻合較好，最大誤差為3.6%，發(fā)生在空化數(shù)σ=0.11處；在空化數(shù)較大時誘導輪揚程系數(shù)較大，空化數(shù)較小時揚程系數(shù)隨空化數(shù)的降低逐漸減小.

圖4進一步給出了不同空化數(shù)下數(shù)值計算的誘導輪空泡形態(tài)與試驗結(jié)果對比，可以看出：數(shù)值計算的空泡形態(tài)與試驗觀測的結(jié)果吻合較好；當σ=0.288時，誘導輪在高轉(zhuǎn)速下形成較小的葉片前緣葉尖空化；隨著空化數(shù)不斷減小，誘導輪前緣空化不斷發(fā)展，當σ=0.092時誘導輪前緣已經(jīng)呈現(xiàn)明顯的云狀空化；當σ=0.060時，誘導輪前緣空化較為嚴重并伴隨著回流渦空化，同時揚程系數(shù)也有明顯的下降.

圖4 不同空化數(shù)下誘導輪空化形態(tài)對比

綜上所述，文中采用的數(shù)值計算方法能夠較為準確地預測誘導輪內(nèi)部的空化流動.

2.2 誘導輪非定常空化流動特性

為了進一步分析誘導輪內(nèi)部非定?？栈鲃犹匦裕槍Ζ?0.053，σ=0.077典型工況下誘導輪空化流動進行非定常數(shù)值計算，并對誘導輪內(nèi)部非定常流動特性進行研究.

圖5為誘導輪流道內(nèi)空化體積脈動的功率譜密度PSD分布，可以看出，誘導輪流道內(nèi)部的空泡體積脈動頻率以低頻為主，在f=14.6，53.6 Hz處存在明顯峰值.

圖5 誘導輪流道空化體積脈動的功率譜密度分布

Fig.5 Power spectral density distribution of cavity volume pulsation in inducer flow passage

圖6為1個旋轉(zhuǎn)周期誘導輪流道內(nèi)空泡分布，可以看出：誘導輪在該工況下空化現(xiàn)象比較明顯，各葉片均有大面積的附著型空化，空化主要發(fā)生在葉片前緣和葉片進口輪緣處，每個葉片上的空泡形態(tài)大小不一，并且隨著時間的推移各葉片空泡形態(tài)也在不斷變化；在t=t0時刻，誘導輪的進口產(chǎn)生回流渦空化；隨著誘導輪的旋轉(zhuǎn)，在t=t0+1T/10時刻，一個回流渦空化潰滅，而另一個回流渦空化卻在緩慢增大；從t=t0+1T/10到t=t0+6T/10時間內(nèi)，該回流渦空泡一直呈現(xiàn)增長的趨勢；從t=t0+7T/10到t=t0+T時間段內(nèi)，回流渦空泡大小又逐漸減小.在整個旋轉(zhuǎn)周期中，回流渦空化旋轉(zhuǎn)方向與誘導輪旋轉(zhuǎn)方向一致，但旋轉(zhuǎn)速度遠小于誘導輪轉(zhuǎn)速.

由于受旋轉(zhuǎn)、湍流、空化影響，誘導輪的內(nèi)部流場具有顯著的壓力脈動特性.為了研究誘導輪內(nèi)部流場的壓力脈動和空泡體積脈動規(guī)律，在誘導輪流域均勻布置24個壓力測點，各個壓力測點的相對位置如圖7所示，其中監(jiān)測點P1,P2,…,P8位于誘導輪進口截面，監(jiān)測點P9,P10,…,P16位于距離誘導輪進口35 mm的葉頂間隙處，監(jiān)測點P17,P18,…,P24則位于距離誘導輪出口截面下游8 mm處.在誘導輪空化流動非定常數(shù)值計算過程中，對24個監(jiān)測點處的壓力進行實時監(jiān)測.

圖6 1個完整旋轉(zhuǎn)周期內(nèi)誘導輪空泡形態(tài)變化

圖7 誘導輪流體域監(jiān)測點布置示意圖

圖8為1個周期內(nèi)誘導輪進口截面相鄰監(jiān)測點P1和P2壓力隨時間變化關(guān)系，可以看出，受誘導輪旋轉(zhuǎn)的影響，2個點的壓力呈現(xiàn)明顯的周期性，壓力幅值在5～50 kPa內(nèi)波動，且波峰和波谷對應(yīng)的壓力值也呈現(xiàn)周期性變化.

圖8 監(jiān)測點P1和P2壓力

Fig.8 Pressure curves in time-domain atP1andP2

進一步研究進口壓力脈動，圖9給出了誘導輪進口截面壓力測點P1,P2,…,P8的功率譜密度分布，可以看出：監(jiān)測點P1,P2,…,P8的功率譜密度分布基本一致，監(jiān)測點的壓力脈動主頻與誘導輪葉片通過頻率相同，為200 Hz，即誘導輪進口的壓力脈動主要受葉片旋轉(zhuǎn)影響；同時在頻率為100，300，400，500 Hz的葉輪轉(zhuǎn)頻倍頻處，功率譜密度分布均有明顯的峰值；在f=50 Hz(誘導輪旋轉(zhuǎn)頻率)處，壓力脈動也有1個較小的尖峰，但明顯小于其他峰值.

圖9 誘導輪進口壓力測點P1—P8的功率譜密度隨時間變化度分布圖

Fig.9 Power spectral density of inlet pressure at monitoring pointsP1toP8

圖10為1個周期內(nèi)誘導輪相鄰監(jiān)測點P9和P10壓力隨時間變化關(guān)系.可以看出，在1個旋轉(zhuǎn)周期內(nèi)，由于誘導輪旋轉(zhuǎn)的影響，2個監(jiān)測點壓力均呈現(xiàn)周期性變化規(guī)律.監(jiān)測點P9,P10,…,P16均會周期性地掠過空化區(qū)域，因此壓力波動特別明顯，最小壓力為飽和蒸氣壓3 169 Pa.監(jiān)測點轉(zhuǎn)過空化區(qū)域后，壓力增大過程中存在1個小的壓力轉(zhuǎn)折區(qū)域，這是由于空泡尾部的局部高壓而產(chǎn)生的.整個旋轉(zhuǎn)周期中，監(jiān)測點P9和P10的壓力波動為3 169～141 000 Pa.

圖10 監(jiān)測點P9和P10壓力隨時間變化

Fig.10 Pressure curves in time-domain atP9andP10

為了進一步研究監(jiān)測點P9,P10,…,P16的壓力脈動情況，圖11給出了誘導輪葉頂間隙截面壓力測點P9,P10,…,P16的功率譜密度分布.監(jiān)測點P9,P10,…,P16位于誘導輪葉輪內(nèi)部，受誘導輪旋轉(zhuǎn)影響更大，因而主頻為葉頻200 Hz.在頻率為100，300，400，500，600，700，800 Hz等整數(shù)倍誘導輪轉(zhuǎn)頻處均有不同程度的壓力脈動峰值，并且400，600，800 Hz等整數(shù)倍葉頻處的峰值，相較于300，500，700 Hz明顯要高.

圖11 導輪葉頂間隙壓力測點P9—P16的功率譜密度分布圖

Fig.11 Power spectral density of inlet pressure at monitoring pointsP9toP16

圖12給出了1個周期內(nèi)誘導輪出口截面相鄰監(jiān)測點P17和P18壓力變化.可以看出，2個監(jiān)測點的壓力均呈現(xiàn)準周期性，壓力幅值為146～167 kPa.由于誘導輪出口截面壓力較高，壓力波動相對較小.

進一步對誘導輪出口壓力脈動進行研究，對監(jiān)測點P17,P18,…,P24壓力進行FFT變換，出口截面監(jiān)測點壓力的功率譜密度分布圖如圖13所示.由圖13a可以看出，壓力脈動主頻為200 Hz，在400，600，800，1 000 Hz等整數(shù)倍葉頻處有明顯的幅值尖峰，此外在14.6，26.8，136.6，148.8 Hz頻率處均有壓力脈動峰值.由圖13b可以看出，各個監(jiān)測點的功率譜密度分布基本一致，在低頻處對應(yīng)幅值大小有差別.

圖12 監(jiān)測點P17和P18壓力隨時間變化

圖13 誘導輪出口截面壓力測點P17—P24的功率譜密度分布

Fig.13 Power spectral density of inlet pressure at pointsP17toP24

3 結(jié) 論

采用基于旋轉(zhuǎn)曲率修正的湍流模型對誘導輪空化特性進行數(shù)值計算,針對在Ω=3 000 r/min，Φ=0.053工況下對DAPAMITO4誘導輪的空化特性進行預測，并與試驗結(jié)果進行對比，分析了誘導輪非定?？栈鲃犹匦砸约傲鲌鰤毫γ}動特性，得到結(jié)論如下：

1) 數(shù)值計算結(jié)果與試驗數(shù)據(jù)吻合較好，誘導輪揚程系數(shù)的最大誤差為3.6%，發(fā)生在空化數(shù)σ=0.11處.空化數(shù)較小時揚程系數(shù)隨空化數(shù)的降低而減小.

2) 針對Φ=0.053，σ=0.077工況下的誘導輪非定?？栈鲃訑?shù)值計算結(jié)果表明，空泡主要分布在葉片前緣和葉片進口輪緣處，每個葉片上的空泡形態(tài)大小不一，并且隨著時間的推移各葉片空泡形態(tài)不斷波動.

3) 誘導輪的內(nèi)部流場具有顯著的壓力脈動特性，進口截面壓力脈動各點的功率譜密度分布基本一致，壓力脈動主頻為200 Hz，與誘導輪葉片通過頻率相同，說明誘導輪的壓力脈動主要受葉片旋轉(zhuǎn)影響.葉頂間隙處壓力脈動受誘導輪旋轉(zhuǎn)影響更大，主頻為200 Hz,在整數(shù)倍誘導輪轉(zhuǎn)頻處均有不同程度的壓力脈動峰值.誘導輪出口壓力脈動主頻仍然為200 Hz，在整數(shù)倍葉頻處有明顯的幅值尖峰.各個壓力測點的功率譜密度分布基本一致，在低頻處對應(yīng)幅值大小有差別.