液氫和液氮繞水翼空化流動特性分析

孫鐵志，魏英杰，王　聰

(哈爾濱工業大學 航天學院，哈爾濱150001)

?

液氫和液氮繞水翼空化流動特性分析

孫鐵志，魏英杰，王聰

(哈爾濱工業大學 航天學院，哈爾濱150001)

摘要:為分析液氫和液氮兩種低溫流體介質的空化特性，通過對CFX軟件二次開發，將Schnerr-Sauer空化模型和液氮、液氫隨溫度變化的物性參數嵌入到CFX求解代碼中，同時耦合求解考慮汽化潛熱影響的能量方程，從而在考慮熱力學效應條件下，開展了液氫和液氮繞水翼空化流動的三維數值模擬研究，并將計算結果與試驗數據進行對比，驗證了數值方法的有效性.計算結果表明,熱力學效應對液氫空化區域壓力和溫度參數變化影響更顯著，在液氮空化核心區域內液相體積分數比液氫中的更小，在空泡尾部閉合區域從汽相向液相轉化迅速.汽-液兩相間質量傳輸特性可作為評估空化區域內溫度、壓力以及相體積分數分布的有效依據.

關鍵詞:液氫；液氮；空化流動；熱力學效應；數值計算

低溫流體如液氫、液氧等被廣泛地用于液體火箭發動機推進系統中.發動機渦輪泵高速旋轉時會使得葉片周圍流體介質發生壓降，當壓力降低到當地飽和蒸汽壓強以下時會引起空化的發生[1-2].汽化潛熱的存在使得空化相變過程產生汽-液兩相間熱量的傳遞，即空化熱力學效應[3].常溫水空化過程通常被視為等溫過程，而低溫流體介質物理屬性對溫度變化敏感，使得熱力學效應在低溫流體空化過程中表現較為顯著[4-5].空化現象直接影響著火箭發動機推力的大小，同時也是制約發射成功與否的關鍵因素之一，所以掌握低溫流體空化的流動特性及其預示方法顯得尤為重要.

Stahl等[6]最早提出了B因子理論，分析了熱力學效應對泵的揚程的影響；后來Sarosdy等[7]開展了水和氟利昂空化的對比試驗，發現了氟利昂液體空化時空泡呈泡霧狀態特性；20世紀70年代Hord[8]系統地開展液氫和液氮的低溫空化試驗，從而加深對低溫流體空化特性的直觀認識.但是由于低溫流體工作環境的限制，大大增加了試驗流場參數的測試難度.近些年數值計算成為研究低溫流體空化的主要手段,Hosangadi 等[9]采用Merkle空化模型對液氫和液氮流體空化流動進行了計算，但是在空化閉合區域計算結果與試驗差值較大；Tseng[10]、馬相孚[11]、張小斌[12]和黃彪[13]等對入口黏度、質量傳輸模型參數敏感性等因素對空化流場結果影響進行了分析，但都是二維數值模擬計算.開展三維數值計算可有效獲取空化流場參數，加深對低溫流體空化特性的全面認識.

本文通過CEL (CFX expression language) 語言將Schnerr-Sauer空化模型嵌入到CFX軟件中，并將液氮、液氫的密度、比熱容、熱傳導系數以及飽和蒸汽壓強等隨溫度變化的物性參數引入到求解代碼中；同時考慮空化過程汽化潛熱的影響，并將其以源相的形式添加到能量方程中，從而形成了一套計算低溫流體空化的三維數值方法，為后續開展全模型誘導輪空化提供技術支撐.同時獲得的液氫和液氮空化流場參數為今后液體火箭發動機推進系統的設計提供參考.

1數值計算方法

1.1基本控制方程和湍流閉合方法

對于考慮熱力學效應的低溫空化流動問題計算的控制方程，除了連續性方程和動量方程外，還包括包含能量源項的能量方程，依次如下：

式中：ρm=αlρl+αvρv為混合相的密度；下標v和l分別代表汽相和液相；下標i和j分別代表坐標方向；u為速度；μ為流體的動力黏度；μt是湍流黏度；keff為熱傳導系數；L為汽化潛熱；cp為比定壓熱容.

選取k-ε兩方程模型作為湍流閉合方法，通過求解湍動能和湍流耗散值可有效預測空化流場湍流特性.

1.2空化模型

式中,ρv隨溫度變化的屬性通過自定義函數引入到CFX中，計算過程中ρv對應于當地溫度下的密度值.

Schnerr-Sauer空化模型[14]主要是基于Rayleigh-Plesset汽(氣)泡動力學方程(1)的推導，得到描述汽(氣)泡生長和潰滅的基本方程：

(1)

式中：pv為泡內壓強；p為氣泡周圍流場壓強；R為汽泡半徑；ρl為液體密度；ζ為液體表面張力.忽略方程(1)右側的二次項和表面張力項可得表征空化過程氣泡半徑RB的變化方程：

基于上述分析得到Schnerr-Sauer空化模型中的蒸發源相和凝結源相表達式為：

(2)

(3)

式中,pv(T)為遠處流場當地溫度下的飽和蒸汽壓強.

1.3計算模型

液氫和液氮繞水翼空化計算模型與Hord試驗一致，試驗結構及尺寸如圖1所示.根據試驗條件，計算域入口采用速度入口，出口為壓力開口，流域邊界及水翼設置為絕熱不可滑移壁面，且各個邊界位置初始液相體積分數設置為1，汽相體積分數為0.同時并在翼型表面分別設置5個溫度監測點和5個壓力監測點，以便將數值計算結果與試驗數據進行對比.

圖1　計算域及邊界條件

圖2給出了計算模型網格劃分，計算域內采用H-型和C-型網格，以提高計算效率；同時在水翼附近網格進行加密處理，以有效捕捉翼型周圍流場參數，網格總數量為104萬，且在計算過程中對網格無關性進行了驗證.

圖2　計算網格劃分

2結果與分析

根據Hord試驗數據隨機選取4種不同工況，其入口流速Uint、遠場溫度T∞和入口空化數σint如表1所示.計算過程中保持計算網格、湍流模型、邊界條件以及試驗參數等條件一致.通過對比兩種介質下空化區域兩相分布特性、翼型周圍壓力和溫度數據等綜合參數，分析液氫和液氮空化特性的差異.

表1　計算工況

為更好地對比分析熱力學效應下液氫和液氮空化特性，圖3給出了兩種介質中汽化壓強和液-汽密度比隨溫度變化曲線.針對本文計算工況，當液氫溫度在20.5 K時，液-汽密度比為50，當溫度降低1 K時，汽化壓強減小28 kPa；當液氮溫度在83 K時，液-汽密度比為95，當溫度降低1 K時，汽化壓強減小18 kPa[15].

表2對比了兩種溫度下液氮和液氫的比定壓熱容、密度比、熱傳導系數以及汽化潛熱等物性參數.

表2液氫和液氮部分物性參數比較

Tab.2Comparison of the physical properties of liquid hydrogen, liquid nitrogen and water

介質比定壓熱容/(J·(kg·K)-1)液/汽密度比熱導率/(W·(m·K)-1)汽化潛熱/(kJ·kg-1)液氫(20K)981657100446液氮(83K)207595134190

圖3　液氫和液氮隨溫度變化的物理屬性

2.1壓力及溫度分布特性

圖4給出了翼型表面壓力和溫度分布的數值計算結果與試驗數據對比，其中定義壓降Δp=p-pv(T∞)，溫降ΔT=T-T∞.考慮到Hord試驗測量過程中溫度不確定誤差為0.2 K，液氫和液氮的壓力測量誤差分別為6.90 kPa和10.34 kPa，可知數值計算得到的溫降和壓力分布與試驗結果吻合較好，從而驗證了數值計算方法的有效性.當水翼周圍發生空化時，由于汽化潛熱的影響使得空化過程不斷從周圍流體吸收熱量，從而導致空化區域的溫度降低，空化強度的差異使得空化區域溫度分布不同，即體現為液氫和液氮發生空化時水翼周圍表現為Δp<0和ΔT<0，而在常溫水空化時溫降和壓降可忽略不計.在液氫254C和260D工況中，254C工況在空化閉合區域壓力恢復到遠場壓力梯度較小，低壓區域較大；液氮290C和296B工況中，雖然296B工況入口空化數較小，但其空化區域較小，最大壓降低于290C工況，可見當入口空化數較小時熱力學效應抑制作用更明顯，所以在液氮空化中入口空化數不能充分體現空化強度的大小.

對比圖4中液氫和液氮壓力和溫度分布可知，在液氮空化低壓區壓力恢復到遠場壓強的壓力變化梯度明顯較大，且在空化閉合區域出現壓力峰值，這可能是空泡閉合位置汽液之間質量傳輸特性的差異引起的；對比溫降數據，可以看出液氮在空化區域溫降較大.在空化熱力學敏感介質中，最大壓降和溫降百分比是評價熱力學影響的重要參數之一，根據圖4的計算結果，表3總結了4種工況下最大壓降(Δp|=(pmin-pv(T∞))/pv(T∞))和溫降(|ΔT|=(Tmin-T∞)/T∞)百分比數據.4種工況下雖然空化數接近，但是液氫最大壓降百分比超過40%，最大溫降百分比大于6%，二者數值均約為液氮的兩倍，可見熱力學效應明顯改變空化區域流場參數的分布，且熱力學效應對液氫空化特性的影響更加顯著.

圖4　數值計算結果與試驗數據對比

%

2.2兩相分布特性及質量傳輸對比

為進一步對比分析兩種介質空化特性的差異，圖5給出了254C和290C工況穩態計算空化流場特性對比.圖5(a)兩種工況下三維空泡形態沿水翼跨度方向基本呈均勻對稱分布，均形成了穩定的空泡，從圖中也可初步看出二者空泡形態稍有差異，在254C液氫工況計算的空泡較厚.同時圖5(a)也給出了水翼表面壓力場分布，可見壓力沿水翼跨度方向(圖1中z軸)基本呈對稱均勻分布，壓力在水翼前緣最大，在水翼靠近前緣位置形成低壓區，并沿流動方向壓力逐漸恢復.從圖中可以看出290C工況低壓區域分布較大，在空化閉合區域壓力梯度變化較大，在空泡內部水翼表面壓力等值線呈U型分布，可見在三維計算中可獲得壁面效應對空泡形態和壓力分布的影響.圖5(a)也給出了流場內部速度矢量分布圖，可以看出在靠近水翼前緣位置速度變化梯度較大，在空泡內部及閉合位置沒看到明顯的回射流現象.

為更有效對比空泡特征，圖5(b)給出了254C和290C工況在水翼中截面壓力分布(下半區域)和液相體積分數分布(上半區域).從壓力云圖分布可以進一步看出254C工況低壓區域覆蓋范圍較小，壓力沿在水翼周圍變化梯度不明顯.在低溫狀態下空化區域壓力分布實際反映著空泡內兩相分布特性，空泡內在液相向汽相轉化過程中，由于發生溫降，使得空泡壓強與當地溫度下飽和蒸汽壓強直接相關；同時液氮高密度比(見圖3)和空化時溫降值較大(見圖4)，從而導致液氮空化流場壓力梯度變化明顯.

從圖5(b)可知兩種工況下空泡長度基本一致，但是空泡形態輪廓和液相體積分數分布有所差異.兩種工況下水翼表面液相最小體積分數都大于零，這主要是由于低溫液氫和液氮屬性對溫度變化敏感，同時液-汽密度比值較小，使得在質量傳輸過程中液相向汽相轉換不充分，這與參考文獻[9,13,16]中結果一致.在液氫254C工況中空泡形態近似呈橢圓形，液相最小體積分數為0.45，在空化閉合區域汽相向液相轉換的凝結過程緩慢，表現為液相體積分數變化梯度較小.而在液氮290C工況中空泡最大厚度發生在空化下游區域，空泡內液相最小體積分數比液氫中小，最小值為0.2，可見在水翼周圍空化強度較大，同時在空泡尾部閉合位置汽相向液相轉化的體積分數變化梯度較大，這也進一步解釋了液氫空化閉合區域壓力變化緩慢和液氮空化空泡尾部出現壓力峰值的原因.同時根據Hord實驗觀測，液氮和液氫空泡呈不透明的泡霧狀，這與水空化空泡汽液界面清晰有明顯不同，計算得到圖5(b)中體積分數分布特性有效地解釋了液氮/液氫空泡呈泡霧狀現象問題.

圖5　液氫和液氮空化流場特性對比

注：彩圖見電子版(http://hit.alljournals.cn)(2016年第8期)

空化過程汽液兩相間質量傳輸率對流場分布特性具有重要影響，為深入分析上述計算結果，圖6給出了254C和290C工況水翼表面蒸發源相和凝結源相的分布曲線.可以看出兩種介質中蒸發源相m-明顯大于凝結源相m+，從而證明了低溫流體中汽化過程速度要比凝結過程快得多，這與常溫水中的結論一致.

圖6　液氫和液氮空化質量傳輸率對比

在圖6(a)中最大蒸發率發生在x=0.003 5位置(圖中①和②)，結合圖4壓力和溫度變化曲線，可知最大溫降和壓降發生在空化區域最大蒸發率位置.圖6(b)中290C工況凝結源相峰值(圖中③位置)明顯較大，且曲線變化陡峭，即體現為液氮空化時在空泡尾部閉合區域汽-液兩相間傳輸劇烈，從而引起該區域壓力出現峰值(見圖4(a))、液相體積分布變化迅速(見圖5(b)).不同工況的計算結果均表明液氮的m-和m+大于液氫，圖6(a)液氮在①位置值約為液氫在②位置的4倍，而根據圖3中兩種介質屬性，290C工況T=83 K時，液-汽密度比為95，液氫254C工況T=20.5 K時，液-汽密度比為50，兩種工況操作溫度下液氮液汽密度比約為液氫的2倍，可見液汽密度比不能完全反映出空化過程兩相間質量傳輸，結合式(2)～(3)可知空化區域兩相密度、體積分數以及壓力分布等因素綜合影響著液汽兩相間質量轉化的大小.

3結論

1)建立了液氫和液氮繞三維水翼空化的數值計算方法.對比空化區域壓力和溫度分布的數值結果和試驗數據，在測量誤差范圍內兩者體現了較好的一致性，從而驗證了計算方法的有效性.

2)熱力學效應對兩種介質空化區域溫度和壓力影響程度不同.熱力學效應引起液氮和液氫空化區域發生壓降和溫降，但是在空化數相近時，液氫的最大溫降和壓降百分比為液氮的2倍.

3)流體介質物理屬性影響著空化區域汽液兩相分布特性.液氫在空化區域液相體積分數較大，且在空泡閉合區域液汽兩相間轉化緩慢；液氮在空泡尾部閉合區域液相體積分數變化迅速.

4)汽液兩相間質量傳輸特性體現著流場參數的分布.最大壓降和溫降發生在蒸發率最大位置；液氮中空泡閉合位置的凝結源相的突變反映此空化區域壓力在恢復到遠場壓力過程中變化梯度較大，同時引起該區域液相體積分布變化迅速.

參考文獻

[1] BATCHELOR G K. An introduction to fluid dynamics[M]. UK: Cambridge University Press, 2000.

[2] JOSEPH D D. Cavitation in a flowing liquid[J]. Physical Review E, 1995, 51(3): R1649-R1650. DOI：10.1103/PhysRevE.51.R1649.

[3] WATANABE S, HIDAKA T, HORIGUCHI H, et al. Analysis of thermodynamic effects on cavitation instabilities[J]. Journal of Fluids Engineering, 2007, 129(9): 1123-1130. DOI:10.1115/1.2754326.

[4] TSENG C C, WEI Yingjie, WANG Guoyu, et al. Modeling of turbulent, isothermal and cryogenic cavitation under attached conditions[J]. Acta Mechanica Sinica, 2010, 26(3): 325-353. DOI：10.1007/s10409-010-0342-7．

[5] 馬相孚. 低溫流體空化特性數值研究[D].哈爾濱：哈爾濱工業大學，2013.

MA Xiangfu. Numerical study on cavitating flow of cryogenic fluids[D]. Harbin: Harbin Institute of Technology, 2013.

[6] STAHL H A, STEPANOFF A J. Thermodynamic aspects of cavitation in centrifugal pumps[J]. ASME Journal of Basic Engineering, 1956, 78(8): 1691-1693.

[7] SAROSDY L R, ACOSTA A J. Note on observations of cavitation in different fluids[J]. Transactions of the ASME, 1961, 83(3):399-400.DOI:10.1115/1.3658979．

[8] HORD J. Cavitation in liquid cryogens II-hydrofoil: NASA Contractor Report: CR-2156[R]. Washington, DC: NASA. 1973.

[9] HOSANGADI A, SCIENTIST P, AHUJA V, et al. Numerical study of cavitation in cryogenic fluids[J]. Journal of Fluids Engineering, 2005, 127(2): 267-281. DOI:10.1115/1.1883238．

[10]TSENGC C, SHYY W. Modeling for isothermal and cryogenic cavitation[J]. International Journal of Heat and Mass Transfer, 2010, 53(1): 513-525. DOI:10.1016/j.ijheatmasstransfer.2009.09.005．

[11]MA Xiangfu, WEI Yingjie, WANG Cong, et al. Simulation of liquid nitrogen around hydrofoil cavitation flow[J]. Journal of Ship Mechanics, 2012,16 (12): 1345-1352.

[12]ZHANG Xiaobin, WU Zhao, XIANG Shijun, et al. Modeling cavitation flow of cryogenic fluids with thermodynamic phase-change theory[J]. Chinese Science Bulletin, 2013, 58(4/5): 567-574. DOI: 10.1007/s11434-012-5463-x．

[13]HUANG Biao, WU Qin, WANG Guoyu. Numerical investigation of cavitating flow in liquid hydrogen[J]. International Journal of Hydrogen Energy, 2014, 39(4): 1698-1709. DOI:10.1016/j.ijhydene.2013.11.025．

[14]SCHNERR G H, SAUER J. Physical and numerical modeling of unsteady cavitation dynamics [C]//Proceedings of the 4th International Conference on Multiphase Flow. New Orleans:[s.n.], 2001. [15]LEMMON E W, HUBER M L, MCLINDEN M O. Reference fluid thermodynamic and transport properties [DB/M]. NIST Standard Reference Database 23, Version 7.0, 2002.

[16]SUN Tiezhi, WEI Yingjie, WANG Cong, et al. Three-dimensional numerical simulation of cryogenic cavitating flows of liquid nitrogen around hydrofoil[J]. Journal of Ship Mechanics. 2014, 18(12): 1434-1443.DOI：10.3969/j.issn.1007-7294.2014.12.003．

(編輯張紅)

doi：10.11918/j.issn.0367-6234.2016.08.024

收稿日期:2015-05-16

基金項目:中央高校基本科研業務費專項資金資助(HIT.NSRIF.201159；黑龍江省自然科學基金(A201409)；哈爾濱市科技創新人才研究專項資金資助(2013RFLXJ007)

作者簡介:孫鐵志(1986—)，男，博士研究生； 魏英杰(1975—)，男，教授，博士生導師；

通信作者:王聰，alanwang@hit.edu.cn

中圖分類號:TJ763

文獻標志碼:A

文章編號:0367-6234(2016)08-0141-06

Study on the behavior of cavitating flows around the hydrofoil in liquid hydrogen and nitrogen

SUN Tiezhi, WEI Yingjie, WANG Cong

(School of Astronautics, Harbin Institute of Technology，Harbin 150001, China)

Abstract:The objective of this study is to analyze the cavitation characteristics in liquid hydrogen and nitrogen. The aim was realized by implanting the Schnerr-Sauer cavitation model and the physical properties of liquid hydrogen and liquid nitrogen at different temperatures into the CFX solver code, and coupling the energy equation considering the latent heat. Then the three-dimensional numerical simulation of cavitating flows was conducted around a hydrofoil in liquid hydrogen and nitrogen, and the experimental results of the pressure and temperature were utilized to validate the numerical strategy. The results show that the thermodynamic effects have more pronounced impact on the pressure and temperature in the cavitation region of liquid hydrogen. The liquid phase volume fraction in liquid nitrogen is smaller in the core cavitation region than that of liquid hydrogen，and the rate of phase transition from vapor to liquid is large in the closure region. The mass transfer rate between liquid and vapor can be used to evaluate the temperature, pressure and phase volume fraction inside the cavity effectively.

Keywords:liquid hydrogen; liquid nitrogen; cavitating flows; thermodynamic effects; numerical simulation

王聰(1966—)，男，教授，博士生導師