李龍賢，丁振曉，吳玉珍

(北京航天動力研究所，北京 100076)

0 引言

誘導(dǎo)輪是液體火箭發(fā)動機(jī)渦輪泵的重要組成部件，可在較低的入口壓力或部分空化條件下工作，同時為主葉輪提供必需的入口壓力防止汽蝕破壞。在空化過程中，工質(zhì)汽化需要吸收汽化潛熱，而這部分熱量只能從工質(zhì)本身獲得，這會導(dǎo)致發(fā)生空化的區(qū)域的溫度要低于未發(fā)生空化的單質(zhì)液相的區(qū)域，這種效應(yīng)叫做工質(zhì)的熱力學(xué)效應(yīng)。熱力學(xué)效應(yīng)對氣泡的形成和發(fā)展的過程會產(chǎn)生抑制作用，進(jìn)而改善誘導(dǎo)輪空化下的性能，低溫工質(zhì)的熱力學(xué)效應(yīng)對誘導(dǎo)輪空化條件下的影響尤為顯著。但由于計(jì)算模型和試驗(yàn)條件的限制，目前國內(nèi)對此問題的認(rèn)識停留在定性的層面。為了更準(zhǔn)確地了解誘導(dǎo)輪在低溫工質(zhì)中的空化特性，必須定量地分析熱力學(xué)效應(yīng)所產(chǎn)生的空化抑制作用，這對于研制抗空化高性能低溫誘導(dǎo)輪、提高低溫液體火箭發(fā)動機(jī)整體性能及火箭有效載荷具有重要意義。日本和歐洲針對低溫工質(zhì)的熱力學(xué)特性在理論、數(shù)值計(jì)算和實(shí)驗(yàn)等方面開展了大量的研究。Fruman等針對低溫?zé)崦袅黧w提出了一種空化修正模型[1]，并以R114為研究對象，進(jìn)行數(shù)值模擬和實(shí)驗(yàn)研究。Franc等基于Rayleigh方法提出了一種考慮熱力學(xué)效應(yīng)的空化模型修正方法[2]，對汽液相間熱交換提出了兩種不同的設(shè)想，并分別在不同溫度下對兩種介質(zhì)R114和水進(jìn)行了試驗(yàn)驗(yàn)證。日本學(xué)者Yoshiki Yoshida等采用液氮為工質(zhì)對熱力學(xué)特性進(jìn)行了研究[3]。國內(nèi)在進(jìn)行低溫工質(zhì)渦輪泵特性試驗(yàn)時大多采用常溫水，忽略了低溫工質(zhì)的熱力學(xué)特性影響。本文結(jié)合國外的研究成果，以液氧為例在低溫?zé)崃W(xué)特性的理論和數(shù)值計(jì)算方面進(jìn)行研究探索。

通常情況下，低溫液體與常溫液體(比如水)的物理特性有很大的區(qū)別，相對于常溫液體，低溫液體可壓縮性較大、氣相與液相的密度差別較小、汽化潛熱較小等?？栈臒崃W(xué)特性影響解釋如下：汽化過程需要吸收汽化潛熱，這部分熱量由液相傳遞給氣液交界面。而只有氣泡內(nèi)的溫度低于液體的溫度形成溫度差時，熱量的傳遞才能進(jìn)行。而工質(zhì)處于飽和狀態(tài)下，溫度和壓力是一一對應(yīng)的。因此，氣泡內(nèi)氣體的壓力(pv(Tc))也要低于液體內(nèi)的壓力(pv(T))。這樣氣泡與參考點(diǎn)處的壓力不平衡值(pv(Tc)-pref)就會降低。所以，存在熱力學(xué)特性時氣泡的增長要慢于不存在熱力學(xué)特性的工質(zhì)。

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

1.1 研究對象

本文的研究對象是某型液體火箭發(fā)動機(jī)氧泵，如圖1所示，氧泵是由進(jìn)口管、誘導(dǎo)輪、導(dǎo)流支座、離心輪、擴(kuò)壓器、蝸殼和出口管組成。誘導(dǎo)輪是泵的重要部件，為三葉片變螺距誘導(dǎo)輪，其實(shí)物照片如圖2所示。

圖1 氧泵幾何模型Fig.1 Oxygen pump geometry model

圖2 氧泵誘導(dǎo)輪Fig.2 Oxygen pump inducer

泵的設(shè)計(jì)參數(shù)：流量124 L/s，轉(zhuǎn)速為18 000 rpm，揚(yáng)程1 135 m，工作介質(zhì)為低溫工質(zhì)液氧。

1.2 數(shù)值計(jì)算條件

1.2.1 網(wǎng)格劃分

為了提高計(jì)算的收斂性，將入口向前延伸500 mm(約為入口直徑的3.5倍)，出口向后延伸200 mm(約為出口直徑的3倍)，在Ansys Workbench下劃分網(wǎng)格，復(fù)雜幾何結(jié)構(gòu)區(qū)采用非結(jié)構(gòu)化四面體網(wǎng)格，進(jìn)/出口管延長區(qū)采用六面體網(wǎng)格。總網(wǎng)格數(shù)約為420萬左右，計(jì)算域網(wǎng)格如圖3所示。

圖3 計(jì)算域網(wǎng)格Fig.3 Calculation field mesh

1.2.2 邊界條件

進(jìn)口設(shè)置為壓力進(jìn)口，出口設(shè)置為質(zhì)量出口，各固體壁面都采用絕熱無滑移壁面邊界條件，誘導(dǎo)輪外殼、誘導(dǎo)輪輪轂、誘導(dǎo)輪葉片、離心輪前后蓋板和離心輪葉片設(shè)置為移動壁面，其他固體壁面都設(shè)置為靜止壁面。

1.2.3 計(jì)算域

整個計(jì)算域由進(jìn)口管延伸段、誘導(dǎo)輪、導(dǎo)流支座、離心輪、蝸殼和出口管延伸段6個子域構(gòu)成，如圖3所示，其中誘導(dǎo)輪和離心輪為轉(zhuǎn)子域，其他為靜子域。

數(shù)值計(jì)算利用CFX計(jì)算軟件，采用有限體積法對計(jì)算區(qū)域進(jìn)行離散，湍流模型選用k-ε雙方程模型，所有控制方程計(jì)算采用基于SIMPLE的標(biāo)準(zhǔn)壓力修正算法。

1.2.4 控制方程

計(jì)算采用基于Rayleigh-Plesset空泡動力學(xué)方程的推導(dǎo)出的Singhal空化模型，基本相為液態(tài)工質(zhì)，第二相為氣態(tài)工質(zhì)。采用雷諾時均方法(RANS)，通過求解流體混合相的連續(xù)性方程、動量方程、能量方程以及第二相(氣相)的體積份額方程和相對速度的代數(shù)表達(dá)式，模擬泵內(nèi)空化流場[4]。

混合模型允許各相之間相互遷移，所以對控制容積的體積分?jǐn)?shù)α可以是0和1之間的任意值，取決于該相所占有的空間。

1.2.4.1 連續(xù)性方程

(1)

其中

ρ=α1ρ1+αvρv

式中：t為時間；ρ為汽液混合流體的密度；α為體積分?jǐn)?shù)；下標(biāo)l為液相；下標(biāo)v為汽相；ρv為真實(shí)氣體密度；V為速度矢量；·(ρV)為對速度的散度。

1.2.4.2 動量方程

(2)

式中：τ為表面力；SM為由體積力引起的動量源項(xiàng)。V·V表示并向量積

(3)

1.2.4.3 能量方程

(4)

式中：e為內(nèi)能；T為溫度；q為與外界的熱交換率；SΦ為耗散函數(shù)，表示流場中粘性切應(yīng)力的所有作用，表示由于變形對流體質(zhì)點(diǎn)做功得到的能量源項(xiàng)，這些功由于機(jī)械作用產(chǎn)生，使流體運(yùn)動并轉(zhuǎn)變成內(nèi)能或熱；SM為由體積力引起的動量源項(xiàng)；V·SM為體積力做功。

1.2.4.4 空泡動力學(xué)方程

本文采用的空化模型為基于均質(zhì)多相質(zhì)量輸運(yùn)方程的空化模型，除混合物質(zhì)量守恒外，只需要增加一個液相或氣相的質(zhì)量守恒方程

(5)

式中Re，Rc分別為汽泡產(chǎn)生和潰滅的質(zhì)量輸運(yùn)源項(xiàng)。

空化動力學(xué)方程采用Rayleigh-Plesset方程

pv(T)-p∞=

(6)

式中：RB為氣泡半徑；pv為飽和壓力；p1為遠(yuǎn)場壓力；ρ1為液體密度；S為氣泡和周圍液體之間的表面張力系數(shù)。

1.2.4.5 Singhal空化模型

在大多數(shù)工程應(yīng)用環(huán)境中，可以認(rèn)為在空化初生階段有足夠的空泡核。這樣就可以把關(guān)注的焦點(diǎn)聚集到空泡的生長和潰滅上。不考慮液相與汽相之間的滑移(即空泡隨液體的流動而運(yùn)動)，在CFD程序中空化模型的實(shí)際形式忽略二階項(xiàng)和表面張力的影響，Rayleigh-Plesset方程可以簡化為

(7)

得空泡體積變化率

(8)

空泡質(zhì)量變化率

(9)

假設(shè)每單位體積有NB個空泡，則

(10)

式中αv為汽相的體積分?jǐn)?shù)。

結(jié)合式(9)和式(10)得單位質(zhì)量的工質(zhì)相間質(zhì)量輸運(yùn)率

(11)

不考慮熱力學(xué)效應(yīng)的影響，則汽化和凝結(jié)速率可表示為如下形式(Singhal空化模型)：

當(dāng)p

(12)

當(dāng)p>pv時，汽泡凝結(jié)為液體

(13)

式中：αnuc為空化核的體積分?jǐn)?shù)；Cvap，Ccon分別為汽化和凝結(jié)源項(xiàng)的修正系數(shù)，通常汽化過程比凝結(jié)過程快得多[5]。經(jīng)驗(yàn)系數(shù)分別取為αnuc=5×104，Cvap=50，Ccon=0.01,RB=1.0×10-6m。

1.2.4.6 基于熱力學(xué)效應(yīng)對空化模型修正

用于描述空化發(fā)生可能性的無量綱參數(shù)空化數(shù)的定義式

(14)

在實(shí)際流體中，空化發(fā)生時的臨界壓力是跟當(dāng)?shù)販囟认嚓P(guān)的，對于當(dāng)?shù)販囟攘黧w的空化數(shù)

(15)

式中Tc為空泡內(nèi)的當(dāng)?shù)販囟?，由?14)和式(15)可得

(16)

將式(15)代入式(5)得

(17)

下面的問題是如何求解ΔT，氣液兩相中的熱平衡關(guān)系可以表示為

ρvvvL=ρ1v1Cp1ΔT

(18)

式中：ρ1和ρv分別為液相工質(zhì)和氣相工質(zhì)的密度；v1和vv分別為液相工質(zhì)和氣相工質(zhì)的體積流量；L為飽和壓力對應(yīng)飽和溫度下的汽化潛熱；Cp1為液相工質(zhì)的定壓比熱容。

在穩(wěn)態(tài)工況的假設(shè)下，溫度的下降主要由相與相間的體積流量之比來評估，也就是所謂的B因子法，這是由Ruggeri和Moore在1969年首先提出的，其表達(dá)式如下

(19)

關(guān)于如何求解B因子，很多學(xué)者進(jìn)行了相關(guān)的研究[6-8]，本文采用Franc提出的方法

vv≈(1-α1)Vδc

(20)

v1≈α1Vδc

(21)

式中δc為葉片上空穴的厚度。B因子可以表示為

(22)

式(20)與式(17)聯(lián)立可得

(23)

將式(21)代入式(16)，可得

(24)

即

(25)

令

(26)

由此可見熱力學(xué)效應(yīng)對空化的影響主要體現(xiàn)在參數(shù)∑，該參數(shù)作為一個判據(jù)，以此來判斷空化過程受熱力學(xué)影響的程度。

采用式(13)對氣液兩相間的質(zhì)量傳輸率進(jìn)行修正，可得

(p

(27)

(p>pv，液化過程)

(28)

1.2.5 工質(zhì)飽和壓力與飽和溫度關(guān)系的擬合

式(25)和式(26)中臨界壓力pv是溫度T的函數(shù)。當(dāng)物質(zhì)達(dá)到飽和狀態(tài)時，飽和壓力pv和溫度T之間存在單值函數(shù)關(guān)系

pv=f(tc)

(29)

這一關(guān)系式就是克勞修斯-克拉貝隆方程[9]

(30)

式中R為氣體常數(shù)。

本文的研究對象氧泵誘導(dǎo)輪采用的真實(shí)工質(zhì)為液態(tài)氧，將不同溫度下液氧的臨界壓力數(shù)值進(jìn)行函數(shù)擬合，得到飽和壓力pv與溫度T的關(guān)系[10]。

設(shè)隨機(jī)向量x對隨機(jī)向量y的多元線性回歸模型為[11]

y=Bx+ζ

(31)

式中：B為未知參數(shù)矩陣；ζ為零均值正態(tài)分布，x與y采用列向量形式表達(dá)。

采用多元線性回歸模型，設(shè)液氧飽和壓力pv與溫度T的向量[1,T,T2,T3,T4,T5,T6]τ，有下列關(guān)系

pv=BOv·[1,T,T2,T3,T4,T5,T6]τ

(32)

將不同溫度下對應(yīng)的飽和壓力作為液氧飽和壓力和溫度對應(yīng)樣本，求得最小二乘估計(jì)

BOv=[b0,b1,b2,b3,b4,b5,b6]=

[4.473×106，-3.087×105，8563，-121.1，

0.91，-0.003 438，5.475×10-6]

得到液氧飽和壓力pv與溫度T的擬合曲線如圖4所示。其關(guān)系式如下

pv=4.473×106-3.087×105T+ 8 563T2-121.1T3+0.91T4- 0.003 438T5+5.475×10-6T6

(33)

圖4 液態(tài)氧飽和壓力與溫度的擬合曲線Fig.4 Fitting curve for liquid oxygen saturation pressure & temperature

2 計(jì)算結(jié)果分析

2.1 Singhal空化模型應(yīng)用于常溫水的數(shù)值計(jì)算

空化模型的修正基于Singhal模型進(jìn)行，首先對Singhal模型對空化流場計(jì)算的準(zhǔn)確性進(jìn)行驗(yàn)證。圖5所示為采用Singhal空化模型通過數(shù)值計(jì)算所得的泵的汽蝕余量-揚(yáng)程曲線，可以看出，由數(shù)值計(jì)算結(jié)果與水力試驗(yàn)數(shù)據(jù)基本吻合。

圖5 揚(yáng)程-NPSH曲線Fig.5 Hydraulic head-NPSH curve

圖6所示為采用Singhal空化模型以水為工質(zhì)計(jì)算得到的不同空化數(shù)下的誘導(dǎo)輪葉片空化區(qū)分布與可視化試驗(yàn)結(jié)果的對比?？栈囼?yàn)為本研究團(tuán)隊(duì)在實(shí)驗(yàn)室親自進(jìn)行并采集到的結(jié)果?？梢钥吹讲捎肧inghal空化模型的空化區(qū)大小、形狀以及在葉片上的分布與同等工況下試驗(yàn)測得結(jié)果非常相似。計(jì)算和試驗(yàn)表明空化區(qū)的空化初生位置位于誘導(dǎo)輪葉片修緣段末端，隨空化數(shù)的減小，空化區(qū)逐漸向葉片根部和后緣處蔓延。采用Singhal空化模型對泵內(nèi)空化流場的預(yù)測比較準(zhǔn)確。

圖6 以水為工質(zhì)數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果的對比Fig.6 Comparison between calculation and test result with water as working medium

2.2 采用熱力學(xué)效應(yīng)修正后的模型計(jì)算結(jié)果

采用Singhal空化模型和熱力學(xué)效應(yīng)修正后的模型對誘導(dǎo)輪流場進(jìn)行數(shù)值計(jì)算，計(jì)算采用的工質(zhì)為液氧。圖7～圖10為不同空化數(shù)下模型修正前后葉片上的空化區(qū)分布對比，從數(shù)值計(jì)算的結(jié)果可以看出，添加熱力學(xué)修正項(xiàng)以后，各個空化數(shù)下計(jì)算的空化區(qū)較修正前的空化模型均有所減小，這與理論分析一致。圖7和圖8是空化數(shù)分別為σ=0.06和σ=0.03時模型修正前后的計(jì)算結(jié)果對比，計(jì)算結(jié)果表明此階段空化區(qū)呈細(xì)長條狀附著在葉片修緣處，添加熱力學(xué)效應(yīng)修正項(xiàng)以后空化區(qū)長度明顯縮短。當(dāng)空化數(shù)σ≤0.022以后(圖9和圖10所示)，空化區(qū)逐漸蔓延到葉片流道中，模型修正前后流場中的空化區(qū)分布差別更為明顯，可見低溫工質(zhì)的熱力學(xué)特性對誘導(dǎo)輪進(jìn)入深度空化的進(jìn)程有很大影響。

圖7 模型修正前后的計(jì)算結(jié)果空化區(qū)分布對比(σ=0.06)Fig.7 Cavitation distribution comparison between non-updated and updated model(σ=0.06)

圖8 模型修正前后的計(jì)算結(jié)果空化區(qū)分布對比(σ=0.03)Fig.8 Cavitation distribution comparison between non-updated and updated model(σ=0.03)

圖9 模型修正前后的計(jì)算結(jié)果空化區(qū)分布對比(σ=0.022)Fig.9 Cavitation distribution comparison between non-updated and updated model(σ=0.022)

圖10 模型修正前后的計(jì)算結(jié)果空化區(qū)分布對比(σ=0.02)Fig.10 Cavitation distribution comparison between non-updated and updated model(σ=0.02)

2.3 數(shù)值計(jì)算結(jié)果與試車數(shù)據(jù)的對比

數(shù)值計(jì)算過程中由于計(jì)算資源的限制，只能將控制參數(shù)進(jìn)行平均化處理，在特定時間點(diǎn)采用定常計(jì)算，通過邊界條件的設(shè)定控制轉(zhuǎn)速、流量和進(jìn)口壓力，對計(jì)算出的泵出口壓力進(jìn)行比較。

圖11所示為發(fā)動機(jī)試車過程中空化開始后泵出口壓力曲線試車結(jié)果、修正前空化模型的計(jì)算結(jié)果和修正后空化模型的計(jì)算結(jié)果對比。試車過程中，在21 s開始程序調(diào)節(jié)降低發(fā)動機(jī)入口壓力考察誘導(dǎo)輪及發(fā)動機(jī)的抗空化性能。可以看到在21 s以前，此時空化數(shù)σ>0.026，誘導(dǎo)輪內(nèi)的空化尚未充分發(fā)展，試驗(yàn)結(jié)果和兩種模型的數(shù)值計(jì)算結(jié)果非常接近。21 s以后，σ>0.025，空化區(qū)加速發(fā)展，修正前的空化模型計(jì)算的泵出口壓力明顯小于試驗(yàn)結(jié)果，說明不考慮熱力學(xué)效應(yīng)的空化模型計(jì)算的泵內(nèi)的空化程度大于實(shí)際情況。采用修正后的空化模型得出的計(jì)算結(jié)果更逼近試驗(yàn)結(jié)果，從圖中曲線可以看出，空化區(qū)開始發(fā)展以后，模型修正后揚(yáng)程下降幅度明顯減小，但計(jì)算結(jié)果要略大于試驗(yàn)結(jié)果，且在工況波動較大的位置(24 s處)與試驗(yàn)數(shù)據(jù)有較大差異，說明空化模型在反映工況的波動時不夠精細(xì)，計(jì)算空化模型仍有很大的改進(jìn)空間。

圖11 數(shù)值計(jì)算結(jié)果與試車數(shù)據(jù)的比較Fig.11 Comparison between numerical calculation and trial run result

圖12所示為采用相似計(jì)算得到的實(shí)際工質(zhì)泵與水試泵汽蝕性能的比較，相似計(jì)算是采用相似定律將試車數(shù)據(jù)和計(jì)算數(shù)據(jù)折算到與水試相同的工況，并與水試泵汽蝕性能曲線進(jìn)行對比?？梢钥吹讲捎靡貉鯙楣べ|(zhì)的試車數(shù)據(jù)較水試數(shù)據(jù)汽蝕性能稍有改善，進(jìn)口所需的NPSHr(必需汽蝕余量)從26.05 m降到約22 m。這是因?yàn)榈蜏毓べ|(zhì)相對于常溫水對空化有一定的抑制作用，同一臺泵驅(qū)動低溫工質(zhì)表現(xiàn)出的抗空化性能要優(yōu)于驅(qū)動常溫水的性能。修正前的模型計(jì)算所得數(shù)據(jù)(近似計(jì)算處理以后)與試車數(shù)據(jù)(近似計(jì)算處理以后)偏差較大，與水試數(shù)據(jù)比較接近；修正后的模型與試車數(shù)據(jù)吻合性較好。

圖12 相似計(jì)算后實(shí)際工質(zhì)泵與水試泵汽蝕性能對比Fig.12 Similarity computation performance comparison between real working medium and water trial

3 結(jié)論

本文考慮熱力學(xué)效應(yīng)對低溫工質(zhì)液氧空化性能的影響，對Singhal空化模型進(jìn)行修正，將修正后的模型編譯成程序模塊應(yīng)用于某型號氫氧火箭發(fā)動機(jī)氧泵的數(shù)值計(jì)算，并將計(jì)算結(jié)果與發(fā)動機(jī)空化性能試車數(shù)據(jù)進(jìn)行對比分析，可以得到以下結(jié)論

1)當(dāng)工質(zhì)的熱力學(xué)效應(yīng)可以忽略不計(jì)時(如常溫水)，應(yīng)用Singhal空化模型計(jì)算誘導(dǎo)輪內(nèi)的空化區(qū)分布與試驗(yàn)觀測到的空化區(qū)分布非常接近，計(jì)算結(jié)果比較準(zhǔn)確。

2)直接應(yīng)用Singhal空化模型計(jì)算低溫工質(zhì)液氧的空化流場，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)偏差較大。

3)采用熱力學(xué)效應(yīng)修正后的空化模型應(yīng)用于氧誘導(dǎo)輪的數(shù)值計(jì)算，能反映出熱力學(xué)效應(yīng)對空化的抑制作用的趨勢以及空化區(qū)與液態(tài)流體之間的溫度差異。

4)數(shù)值計(jì)算結(jié)果與試車結(jié)果對比，應(yīng)用修正后的空化模型對參數(shù)穩(wěn)定工況下的計(jì)算結(jié)果與試車數(shù)據(jù)吻合度較好。相似計(jì)算處理以后，發(fā)現(xiàn)采用液氧為工質(zhì)的泵汽蝕性能較水試泵汽蝕性能略有改善。