基于熱力學效應修正的誘導輪空化模型研究

李龍賢，丁振曉，吳玉珍

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

0 引言

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

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

1 計算模型數值分析

1.1 研究對象

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

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

圖2 氧泵誘導輪Fig.2 Oxygen pump inducer

泵的設計參數：流量124 L/s，轉速為18 000 rpm，揚程1 135 m，工作介質為低溫工質液氧。

1.2 數值計算條件

1.2.1 網格劃分

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

圖3 計算域網格Fig.3 Calculation field mesh

1.2.2 邊界條件

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

1.2.3 計算域

整個計算域由進口管延伸段、誘導輪、導流支座、離心輪、蝸殼和出口管延伸段6個子域構成，如圖3所示，其中誘導輪和離心輪為轉子域，其他為靜子域。

數值計算利用CFX計算軟件，采用有限體積法對計算區域進行離散，湍流模型選用k-ε雙方程模型，所有控制方程計算采用基于SIMPLE的標準壓力修正算法。

1.2.4 控制方程

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

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

1.2.4.1 連續性方程

(1)

其中

ρ=α1ρ1+αvρv

式中：t為時間；ρ為汽液混合流體的密度；α為體積分數；下標l為液相；下標v為汽相；ρv為真實氣體密度；V為速度矢量；·(ρV)為對速度的散度。

1.2.4.2 動量方程

(2)

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

(3)

1.2.4.3 能量方程

(4)

式中：e為內能；T為溫度；q為與外界的熱交換率；SΦ為耗散函數，表示流場中粘性切應力的所有作用，表示由于變形對流體質點做功得到的能量源項，這些功由于機械作用產生，使流體運動并轉變成內能或熱；SM為由體積力引起的動量源項；V·SM為體積力做功。

1.2.4.4 空泡動力學方程

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

(5)

式中Re，Rc分別為汽泡產生和潰滅的質量輸運源項。

空化動力學方程采用Rayleigh-Plesset方程

pv(T)-p∞=

(6)

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

1.2.4.5 Singhal空化模型

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

(7)

得空泡體積變化率

(8)

空泡質量變化率

(9)

假設每單位體積有NB個空泡，則

(10)

式中αv為汽相的體積分數。

結合式(9)和式(10)得單位質量的工質相間質量輸運率

(11)

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

當p

(12)

當p>pv時，汽泡凝結為液體

(13)

式中：αnuc為空化核的體積分數；Cvap，Ccon分別為汽化和凝結源項的修正系數，通常汽化過程比凝結過程快得多[5]。經驗系數分別取為αnuc=5×104，Cvap=50，Ccon=0.01,RB=1.0×10-6m。

1.2.4.6 基于熱力學效應對空化模型修正

用于描述空化發生可能性的無量綱參數空化數的定義式

(14)

在實際流體中，空化發生時的臨界壓力是跟當地溫度相關的，對于當地溫度流體的空化數

(15)

式中Tc為空泡內的當地溫度，由式(14)和式(15)可得

(16)

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

(17)

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

ρvvvL=ρ1v1Cp1ΔT

(18)

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

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

(19)

關于如何求解B因子，很多學者進行了相關的研究[6-8]，本文采用Franc提出的方法

vv≈(1-α1)Vδc

(20)

v1≈α1Vδc

(21)

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

(22)

式(20)與式(17)聯立可得

(23)

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

(24)

即

(25)

令

(26)

由此可見熱力學效應對空化的影響主要體現在參數∑，該參數作為一個判據，以此來判斷空化過程受熱力學影響的程度。

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

(p

(27)

(p>pv，液化過程)

(28)

1.2.5 工質飽和壓力與飽和溫度關系的擬合

式(25)和式(26)中臨界壓力pv是溫度T的函數。當物質達到飽和狀態時，飽和壓力pv和溫度T之間存在單值函數關系

pv=f(tc)

(29)

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

(30)

式中R為氣體常數。

本文的研究對象氧泵誘導輪采用的真實工質為液態氧，將不同溫度下液氧的臨界壓力數值進行函數擬合，得到飽和壓力pv與溫度T的關系[10]。

設隨機向量x對隨機向量y的多元線性回歸模型為[11]

y=Bx+ζ

(31)

式中：B為未知參數矩陣；ζ為零均值正態分布，x與y采用列向量形式表達。

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

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

(32)

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

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所示。其關系式如下

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

(33)

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

2 計算結果分析

2.1 Singhal空化模型應用于常溫水的數值計算

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

圖5 揚程-NPSH曲線Fig.5 Hydraulic head-NPSH curve

圖6所示為采用Singhal空化模型以水為工質計算得到的不同空化數下的誘導輪葉片空化區分布與可視化試驗結果的對比。空化試驗為本研究團隊在實驗室親自進行并采集到的結果。可以看到采用Singhal空化模型的空化區大小、形狀以及在葉片上的分布與同等工況下試驗測得結果非常相似。計算和試驗表明空化區的空化初生位置位于誘導輪葉片修緣段末端，隨空化數的減小，空化區逐漸向葉片根部和后緣處蔓延。采用Singhal空化模型對泵內空化流場的預測比較準確。

圖6 以水為工質數值計算結果與試驗結果的對比Fig.6 Comparison between calculation and test result with water as working medium

2.2 采用熱力學效應修正后的模型計算結果

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

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

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

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

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

2.3 數值計算結果與試車數據的對比

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

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

圖11 數值計算結果與試車數據的比較Fig.11 Comparison between numerical calculation and trial run result

圖12所示為采用相似計算得到的實際工質泵與水試泵汽蝕性能的比較，相似計算是采用相似定律將試車數據和計算數據折算到與水試相同的工況，并與水試泵汽蝕性能曲線進行對比。可以看到采用液氧為工質的試車數據較水試數據汽蝕性能稍有改善，進口所需的NPSHr(必需汽蝕余量)從26.05 m降到約22 m。這是因為低溫工質相對于常溫水對空化有一定的抑制作用，同一臺泵驅動低溫工質表現出的抗空化性能要優于驅動常溫水的性能。修正前的模型計算所得數據(近似計算處理以后)與試車數據(近似計算處理以后)偏差較大，與水試數據比較接近；修正后的模型與試車數據吻合性較好。

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

3 結論

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

1)當工質的熱力學效應可以忽略不計時(如常溫水)，應用Singhal空化模型計算誘導輪內的空化區分布與試驗觀測到的空化區分布非常接近，計算結果比較準確。

2)直接應用Singhal空化模型計算低溫工質液氧的空化流場，計算結果與試驗數據偏差較大。

3)采用熱力學效應修正后的空化模型應用于氧誘導輪的數值計算，能反映出熱力學效應對空化的抑制作用的趨勢以及空化區與液態流體之間的溫度差異。

4)數值計算結果與試車結果對比，應用修正后的空化模型對參數穩定工況下的計算結果與試車數據吻合度較好。相似計算處理以后，發現采用液氧為工質的泵汽蝕性能較水試泵汽蝕性能略有改善。