黃彪， 樊亞丁， 梁文棟， 吳欽， 王國(guó)玉

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

液氧、液氫等低溫液體燃料是大推力運(yùn)載火箭發(fā)動(dòng)機(jī)的動(dòng)力選擇之一. 我國(guó)新一代液體燃料火箭發(fā)動(dòng)機(jī)的研制要求其核心部件渦輪泵有更高的轉(zhuǎn)速和更低的入口壓力從而提升渦輪泵功率，進(jìn)而提升整機(jī)性能[1]. 為了減小渦輪泵內(nèi)的空化，常在渦輪泵上加裝誘導(dǎo)輪進(jìn)行加壓以提高發(fā)動(dòng)機(jī)性能，因此渦輪泵內(nèi)空化帶來(lái)的危害也主要集中在誘導(dǎo)輪上. 而液氧等低溫介質(zhì)發(fā)生空化時(shí)受熱力學(xué)效應(yīng)影響，其誘導(dǎo)輪空化特性也與常溫水有所不同，因此對(duì)于低溫介質(zhì)誘導(dǎo)輪空化性能的研究有著重要意義. 國(guó)內(nèi)外對(duì)于誘導(dǎo)輪空化做了大量的實(shí)驗(yàn)及數(shù)值計(jì)算研究. 王小波等[2-4]針對(duì)采用常溫水和液氫作為介質(zhì)的誘導(dǎo)輪進(jìn)行流場(chǎng)分析，結(jié)果表明液氫作為介質(zhì)的誘導(dǎo)輪初生空化數(shù)較低.

本文采用經(jīng)過(guò)熱力學(xué)修正的Kubota空化模型[5]和FBM湍流模型[6]對(duì)常溫水和液氧作為介質(zhì)時(shí)的誘導(dǎo)輪空化特性進(jìn)行數(shù)值模擬，分析了兩種介質(zhì)下的誘導(dǎo)輪空化流場(chǎng)特性. 通過(guò)對(duì)比兩種介質(zhì)誘導(dǎo)輪的空化特性曲線，獲得了液氧熱力學(xué)效應(yīng)對(duì)誘導(dǎo)輪空化特性的影響規(guī)律. 針對(duì)熱力學(xué)效應(yīng)隨溫度的變化規(guī)律，通過(guò)引入空化相似參數(shù)，建立了一種誘導(dǎo)輪空化特性理論預(yù)測(cè)方法，通過(guò)已知溫度下的空化特性曲線對(duì)其他溫度工況空化特性進(jìn)行預(yù)測(cè).

1 數(shù)值計(jì)算方法

1.1 基本控制方程

本文采用均相流模型計(jì)算氣液兩相流動(dòng)，則基本控制方程可表述為

(1)

(2)

(3)

(4)

式中:ρm為混合介質(zhì)密度，ρm=ρlαl+ρvαv，其中ρl和ρv為液相和汽相的密度；αl和αv為液相和汽相體積分?jǐn)?shù)；U為速度；x為坐標(biāo)軸；下標(biāo)i,j,k代表坐標(biāo)軸方向；p為壓強(qiáng)；μm和μtur為混合介質(zhì)的層流和湍流黏性系數(shù)；hm為混合介質(zhì)的焓；fv為水蒸氣的質(zhì)量分?jǐn)?shù)；Lev為汽化潛熱；Prtur和Prlam為層流和湍流普朗特?cái)?shù)；m+，m-分別為凝結(jié)項(xiàng)和蒸發(fā)項(xiàng).

1.2 FBM湍流模型

為更好地模擬誘導(dǎo)輪內(nèi)部的流動(dòng)現(xiàn)象，數(shù)值計(jì)算采用濾波器湍流模型(FBM)[7]：

Pt-ρε+Pk b

(5)

(6)

(7)

式中:k、ε分別為湍動(dòng)能和湍流耗散率；Pt為湍動(dòng)能生成項(xiàng)；μt為湍流黏性系數(shù)；Pkb、Pεb分別為k、ε的浮力作用項(xiàng). 模型常數(shù)分別為σk=1，σε=1，Cε1=1.44，Cε2=1.92，Cμ=0.09，λ=2mm[5].

1.3 空化模型

本文采用修正的Kubota模型以考慮溫度變化對(duì)質(zhì)量傳輸過(guò)程的影響，即在Kubota空化模型的蒸發(fā)和凝結(jié)項(xiàng)中添加了熱力學(xué)效應(yīng)項(xiàng)以考慮溫度變化對(duì)質(zhì)量傳輸過(guò)程的影響，修正后的蒸發(fā)項(xiàng)和凝結(jié)項(xiàng)分別為

m-=

(1-αv-αf)，

(8)

m+=

(9)

式中:Lev為汽化潛熱;α為熱擴(kuò)散率;cp為比熱;在水的計(jì)算中，Ce，Cv分別取50，0.01；在液氧計(jì)算中，Cε，Cv分別取1，0.000 2.

研究表明湍動(dòng)能對(duì)空化產(chǎn)生重要的影響[8]，在上述空化模型中，采用文獻(xiàn)[8]中提出的方法來(lái)計(jì)算湍動(dòng)能k對(duì)當(dāng)?shù)仄瘔簭?qiáng)的影響：

pv=[pv(Tl)+pturb/2]=pv(Tl)+0.39ρmk

(10)

式中pv(Tl)為當(dāng)?shù)仫柡驼羝麎?

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

計(jì)算采用三葉片誘導(dǎo)輪模型，相關(guān)結(jié)構(gòu)參數(shù)如表1所示. 數(shù)值計(jì)算區(qū)域由入口管道，誘導(dǎo)輪以及出口管道三部分構(gòu)成，如圖1所示. 為了容納誘導(dǎo)輪葉尖處的回流和減小誘導(dǎo)輪出口流動(dòng)的不穩(wěn)定性，入口管道和出口管道分別設(shè)為8D和5D. 給定壓力入口條件，根據(jù)空化數(shù)σ(σ=(Pin-Pv)/(0.5ρlu2)；Pin為入口壓力；ρl為介質(zhì)密度；u為誘導(dǎo)輪葉尖速度)換算得到；出口邊界為質(zhì)量流量條件，根據(jù)流量系數(shù)換算得到；誘導(dǎo)輪葉片表面均設(shè)置為無(wú)滑移壁面；誘導(dǎo)輪與入口管道和出口管道的動(dòng)靜交界面均設(shè)置為interface面.

表1 誘導(dǎo)輪幾何參數(shù)Tab.1 Geometric parameters of inducer

圖1 數(shù)值計(jì)算區(qū)域Fig.1 Computational domain

采用6面體結(jié)構(gòu)化網(wǎng)格對(duì)入口和出口管道進(jìn)行網(wǎng)格劃分，采用4面體網(wǎng)格對(duì)誘導(dǎo)輪進(jìn)行網(wǎng)格劃分，如圖2所示. 為了分析網(wǎng)格數(shù)對(duì)計(jì)算結(jié)果的影響，圖3給出了轉(zhuǎn)速為9 000 r/min，流量系數(shù)φ為0.188，空化數(shù)為0.067,工況揚(yáng)程系數(shù)ψ隨網(wǎng)格數(shù)N的變化. 可以看出揚(yáng)程系數(shù)在網(wǎng)格數(shù)達(dá)到158萬(wàn)后穩(wěn)定. 考慮計(jì)算經(jīng)濟(jì)性，本文選取的計(jì)算域網(wǎng)格數(shù)為158萬(wàn).

圖2 網(wǎng)格劃分Fig.2 Mesh distribution

圖3 網(wǎng)格無(wú)關(guān)性驗(yàn)證Fig.3 Grid independence validation

2 結(jié)果分析

2.1 液氧誘導(dǎo)輪空化特性研究

圖4給出了n=9 000 r/min，φ=0.188時(shí)不同流體介質(zhì)誘導(dǎo)輪揚(yáng)程系數(shù)隨空化數(shù)的演變規(guī)律，其中橫坐標(biāo)為空化數(shù)，縱坐標(biāo)為量綱一的揚(yáng)程系數(shù)ψ/ψ0，ψ0為相同工況下單相流動(dòng)誘導(dǎo)輪的揚(yáng)程系數(shù). 對(duì)比流體介質(zhì)為298 K水時(shí)的實(shí)驗(yàn)和數(shù)值計(jì)算結(jié)果可以發(fā)現(xiàn)，數(shù)值預(yù)測(cè)揚(yáng)程系數(shù)隨空化數(shù)變化趨勢(shì)與實(shí)驗(yàn)結(jié)果基本一致，分別選取揚(yáng)程下降10%、20%、30%的工況，則實(shí)驗(yàn)和數(shù)值計(jì)算所對(duì)應(yīng)的空化數(shù)分別為0.032、0.026、0.024以及0.036、0.031、0.027，誤差分別為12.5%，19.2%以及12.5%，說(shuō)明本文數(shù)值計(jì)算方法可以較好地模擬誘導(dǎo)輪空化發(fā)展情況.

圖4 誘導(dǎo)輪揚(yáng)程系數(shù)曲線Fig.4 Head coefficients curve of inducer

對(duì)比流體介質(zhì)為298 K水和100.3 K液氧的誘導(dǎo)輪揚(yáng)程系數(shù)曲線可以看出，當(dāng)介質(zhì)為水時(shí)，臨界空化數(shù)σc為0.048，即空化數(shù)下降至0.048后，揚(yáng)程開始下降，而當(dāng)介質(zhì)為液氧時(shí)，臨界空化數(shù)為0.036，液氧的臨界空化數(shù)較水下降了25%. 當(dāng)空化數(shù)繼續(xù)減小時(shí)，揚(yáng)程系數(shù)劇烈下降，當(dāng)揚(yáng)程下降達(dá)30%時(shí)，水和液氧的空化數(shù)分別為0.027和0.017，液氧的空化數(shù)較水減小37%. 因此，相同流動(dòng)工況下液氧有效抑制了誘導(dǎo)輪空化的發(fā)生和發(fā)展.

圖5給出了兩種流體介質(zhì)誘導(dǎo)輪葉片空化面積隨空化數(shù)變化規(guī)律，其中橫坐標(biāo)為空化數(shù)σ，縱坐標(biāo)為葉片空化區(qū)域面積A與葉片面積A0的比值. 空化數(shù)較大時(shí)，空化較弱，兩種介質(zhì)誘導(dǎo)輪葉片空化面積相近；隨著空化數(shù)的減小，誘導(dǎo)輪葉片空化面積逐漸增加，其中，介質(zhì)為水的誘導(dǎo)輪葉片空化面積增長(zhǎng)速度較快，液氧葉片空化面積增長(zhǎng)較為緩慢，這是因?yàn)榻橘|(zhì)為液氧的誘導(dǎo)輪空化受到抑制.

圖5 兩種介質(zhì)空化面積隨空化數(shù)變化規(guī)律Fig.5 Cavitation area at different cavitation numbers

為了進(jìn)一步分析空化對(duì)誘導(dǎo)輪的影響，圖6對(duì)比了不同空化數(shù)下0.9倍葉高處兩種介質(zhì)誘導(dǎo)輪葉片壓力分布，其中橫坐標(biāo)為誘導(dǎo)輪內(nèi)部流線位置，0為葉片前緣，1為葉片尾緣；縱坐標(biāo)Cp為壓力系數(shù)Cp=(P-Pin)/(0.5ρlu2)，其中P為當(dāng)?shù)貕毫?σ=0.067，0.048時(shí)，吸力面部分區(qū)域壓力下降至飽和蒸氣壓，產(chǎn)生低壓區(qū)，介質(zhì)為水的誘導(dǎo)輪葉片吸力面低壓區(qū)相比介質(zhì)為液氧的誘導(dǎo)輪更大. 當(dāng)σ=0.027時(shí)，對(duì)于介質(zhì)為水的誘導(dǎo)輪空化較為嚴(yán)重，使得葉片壓力面也產(chǎn)生較大低壓區(qū)；對(duì)于介質(zhì)為液氧的誘導(dǎo)輪，此時(shí)其壓力面低壓區(qū)域較小，說(shuō)明以液氧作為流體介質(zhì)時(shí)誘導(dǎo)輪空化發(fā)展更為緩慢.

圖6 0.9倍葉高處葉片載荷分布對(duì)比Fig.6 Comparison of pressure loading at 90% radial span

為了進(jìn)一步分析兩種介質(zhì)誘導(dǎo)輪空化特性，圖7給出了σ=0.027時(shí)不同介質(zhì)誘導(dǎo)輪葉片的溫度分布云圖. 可以看出，當(dāng)流體介質(zhì)為水時(shí)，誘導(dǎo)輪葉片溫度并沒有明顯變化；當(dāng)流體介質(zhì)為液氧時(shí)，誘導(dǎo)輪葉片在虛線圍成的空化區(qū)域內(nèi)產(chǎn)生了明顯的溫降，最大幅度達(dá)2 K，從而降低了空化區(qū)的飽和蒸氣壓，抑制了空化的發(fā)展.

圖7 兩種介質(zhì)誘導(dǎo)輪葉片溫度分布(σ=0.027)Fig.7 Temperature distribution on blade with two kinds of fluids(σ=0.027)

液氧的熱力學(xué)效應(yīng)是由流體氣液密度比，熱傳導(dǎo)系數(shù)等物質(zhì)屬性造成的. 由流體物質(zhì)屬性而導(dǎo)出的熱力學(xué)參數(shù)Σ可以很好地表示介質(zhì)的熱力學(xué)效應(yīng)強(qiáng)弱[8-10].Σ的表達(dá)式為

(11)

式中:ρl為液相密度；ρv為氣相密度；Lev為汽化潛熱;α為熱擴(kuò)散率;cp為比熱;T∞為介質(zhì)溫度. 表2給出了液氧和水在不同溫度條件下的熱力學(xué)參數(shù)值，可以看出液氧的Σ值相比水體高出幾個(gè)數(shù)量級(jí)，有顯著的熱力學(xué)效應(yīng).

表2 不同溫度液氧和水Σ值

Tab.2 Value ofΣof liquid oxygen and water at different temperature

液氧水T/KΣ/(m·s-3/2)T/KΣ/(m·s-3/2)801 9442831.9409015 7272932.444100.380 3782987.345110299 57930311.135

2.2 誘導(dǎo)輪空化特性理論預(yù)測(cè)與分析

對(duì)于同一低溫流體介質(zhì)的誘導(dǎo)輪，介質(zhì)溫度不同，所體現(xiàn)出的空化熱力學(xué)效應(yīng)也不同，因此對(duì)于葉輪機(jī)械不同工況下空化特性的理論預(yù)測(cè)就有著重要的工程意義. 如圖8所示，本節(jié)擬基于已知溫度T1、T2條件下的誘導(dǎo)輪空化特性曲線，建立不同溫度條件下誘導(dǎo)輪空化特性曲線的理論預(yù)測(cè)方法，進(jìn)一步分析熱力學(xué)效應(yīng)對(duì)誘導(dǎo)輪空化特性的影響規(guī)律.

傳統(tǒng)空化數(shù)的表達(dá)式為

(12)

圖8 不同溫度誘導(dǎo)輪空化預(yù)測(cè)方法示意圖Fig.8 Schematic diagram of prediction method of inducercavitation at different temperature

由于熱力學(xué)效應(yīng)的影響，空化區(qū)域流體介質(zhì)的實(shí)際飽和蒸氣壓低于給定溫度下的流體飽和蒸氣壓. 因此考慮空化熱力學(xué)效應(yīng)對(duì)傳統(tǒng)空化數(shù)進(jìn)行修正，得到空化相似參數(shù)σc,min為

(13)

式中Pmin為空化區(qū)域最小壓力.

不同工況滿足空化的相似是應(yīng)用式(13)的條件. 對(duì)于誘導(dǎo)輪而言，量綱一的揚(yáng)程系數(shù)ψ/ψ0表征著流道空化程度，可以認(rèn)為在流量系數(shù)φ一定的條件下，對(duì)于以兩種不同溫度流體作為介質(zhì)工作的誘導(dǎo)輪，若他們的量綱的揚(yáng)程系數(shù)ψ/ψ0相同，則此時(shí)兩種工況所對(duì)應(yīng)的空化相似參數(shù)σc,min相同，因此有

σ1+ΔσV1=σ2+ΔσV2

(14)

由克勞修斯-克拉伯龍方程可得：

(15)

式中：J為熱工當(dāng)量；gc為量綱一的常數(shù)；Lev為汽化潛熱；ΔT為飽和蒸氣壓下降ΔP所對(duì)應(yīng)的溫度的變化；Vv和Vl分別為汽相和液相的比體積. 由于Vv?Vl，Vl通常可忽略.

同時(shí)，根據(jù)空化區(qū)域的熱平衡方程

ρvvvLev=ρlvlClΔT

(16)

式中：vv和vl為參與換熱的汽相和液相的體積變化率；Cl為液相比熱. 聯(lián)立式(15)(16)可得：

(17)

式中：Σ為上文提到的熱力學(xué)參數(shù)；α為流體介質(zhì)的熱擴(kuò)散率. 當(dāng)流量系數(shù)φ和量綱一的揚(yáng)程系數(shù)ψ/ψ0相同時(shí)，有[12]

(18)

其中N1，N2為兩種工況的誘導(dǎo)輪轉(zhuǎn)速. 結(jié)合式(13)(17)(18)可得

(19)

聯(lián)立式(14)(19)可解出ΔσV1，ΔσV2的值. 同理對(duì)于待預(yù)測(cè)工況3，將σ3和ΔσV3式(14)(19)可解出二者的值，達(dá)到預(yù)測(cè)的目的.

應(yīng)用上述方法，圖9基于80K和100.3K液氧介質(zhì)誘導(dǎo)輪空化特性曲線理論預(yù)測(cè)了不同溫度液氧介質(zhì)誘導(dǎo)輪空化特性曲線. 對(duì)比108K液氧介質(zhì)誘導(dǎo)輪空化特性預(yù)測(cè)結(jié)果與數(shù)值計(jì)算結(jié)果發(fā)現(xiàn)，理論預(yù)測(cè)結(jié)果與數(shù)值計(jì)算結(jié)果吻合較好. 對(duì)比不同溫度工況下誘導(dǎo)輪空化特性曲線可以看出，隨著溫度的升高，相同揚(yáng)程系數(shù)對(duì)應(yīng)的空化數(shù)逐漸減小，誘導(dǎo)輪空化逐漸受到抑制，熱力學(xué)效應(yīng)越來(lái)越明顯.

圖9 不同溫度液氧誘導(dǎo)輪空化數(shù)值及預(yù)測(cè)結(jié)果Fig.9 Numerical and prediction results of inducercavitation at different temperature

3 結(jié) 論

采用FBM湍流模型和基于熱力學(xué)修正的Kubota空化模型針對(duì)常溫水和液氧兩種介質(zhì)的誘導(dǎo)輪空化性能進(jìn)行了數(shù)值計(jì)算研究. 常溫水的揚(yáng)程性能與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比，結(jié)果表明兩種介質(zhì)誘導(dǎo)輪揚(yáng)程隨空化數(shù)變化趨勢(shì)一致. 相對(duì)于常溫水，以液氧作為介質(zhì)時(shí)液氧誘導(dǎo)輪的臨界空化數(shù)均有所下降. 液氧是一種強(qiáng)熱力學(xué)效應(yīng)介質(zhì)，在液氧空化區(qū)域液相和氣相會(huì)強(qiáng)烈的熱量交換，使得局部溫度下降. 產(chǎn)生這種現(xiàn)象的原因在于液氧的熱力學(xué)效應(yīng). 對(duì)比液氧誘導(dǎo)輪葉片空化和溫度分布可知空化造成了局部的溫降，進(jìn)而導(dǎo)致當(dāng)?shù)仫柡驼魵鈮合陆?，因此抑制了空化的發(fā)展. 通過(guò)引入空化相似參數(shù)σc,min，建立了誘導(dǎo)輪空化特性理論預(yù)測(cè)方法. 該方法可以通過(guò)兩種已知溫度工況下誘導(dǎo)輪揚(yáng)程特性曲線對(duì)其他溫度工況的誘導(dǎo)輪揚(yáng)程特性進(jìn)行預(yù)測(cè). 預(yù)測(cè)結(jié)果與數(shù)值計(jì)算結(jié)果吻合較好，驗(yàn)證了該方法的可行性.