液氫空間貯存過程膜態(tài)沸騰數(shù)值模擬

王嬌嬌 厲彥忠 王 磊

（1 航天低溫推進(jìn)劑技術(shù)國家重點(diǎn)實(shí)驗(yàn)室 北京 100028）

（2 西安交通大學(xué)制冷與低溫工程系 西安 710049）

1 引 言

隨著中國航天技術(shù)的快速發(fā)展,液氫因具有高能量密度和無毒無污染的特性而被認(rèn)為是未來深空探測的首選推進(jìn)劑。由于其飽和溫度低,與環(huán)境溫差較大,無論在地面還是空間工況下,都易發(fā)生相變現(xiàn)象。在空間發(fā)動機(jī)二次啟動之前,需將連接管路和局部件冷卻至低于飽和溫度,以減少空間中液氫的浪費(fèi),使用過冷液氫將管路充分浸泡是較優(yōu)的預(yù)冷方案。研究表明,在此過程中,膜態(tài)沸騰是占主導(dǎo)作用的換熱機(jī)理。因此,為了實(shí)現(xiàn)液氫的安全高效應(yīng)用,揭示在不同重力水平下液氫過冷膜態(tài)沸騰的換熱特性具有重要意義。

相關(guān)學(xué)者針對膜態(tài)沸騰開展了廣泛研究,所采用的加熱結(jié)構(gòu)包括球形、加熱絲和平板。相較而言,關(guān)于過冷膜態(tài)沸騰,特別是微重力條件下的低溫流體沸騰的研究較薄弱。Bradfiel 等[1]和Farahat 等[2]分別使用水和鈉作為實(shí)驗(yàn)流體,研究了球型加熱面的過冷膜態(tài)沸騰現(xiàn)象。Fan 等[3]采用具有超疏水表面的球體進(jìn)行淬火實(shí)驗(yàn)研究了液體過冷度的影響。Dhir 和Purohit[4]開展了過冷膜態(tài)沸騰實(shí)驗(yàn),并考慮了液體過冷度,球體初始溫度和球體物性的影響。為了研究相界面波動情況,Vijaykumar 和Dhir[5]對垂直表面的過冷膜態(tài)沸騰進(jìn)行了實(shí)驗(yàn)研究。此外,Ede 和Siviour[6]針對加熱細(xì)絲在過冷水中膜態(tài)沸騰開展實(shí)驗(yàn),指出換熱量隨著過冷度的增加而增加。Sakurai 等[7]采用不同流體、加熱絲直徑和不同壁面過熱度開展了一系列實(shí)驗(yàn)研究。Rousselet 等[8]的實(shí)驗(yàn)結(jié)果表明,液體過冷度的影響作用隨著壁面過熱度的增加而降低。

隨著計(jì)算機(jī)技術(shù)的發(fā)展,數(shù)值模擬已經(jīng)成為研究相變現(xiàn)象的一種重要方法。研究者為追蹤界面變化提出了多種方法,包括拉格朗日移動網(wǎng)格法,Eularian流體體積（VOF）方法,Eularian 水平集（LS）方法,Lagrangian/ Eularian 界面前跟蹤（FT）方法和Lattice-Boltzmann 方法（LBM）。Son 和Dhir[9]首先使用數(shù)值計(jì)算方法模擬池沸騰現(xiàn)象,他們使用拉格朗日運(yùn)動網(wǎng)格法捕獲界面,但未獲得整個(gè)界面運(yùn)動。VOF 方法作為一種常用方法已被許多研究膜態(tài)沸騰的相關(guān)學(xué)者采用。文獻(xiàn)中多采用二維幾何形狀,結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好。此外,Tomar 等[10]和Ningegowda 等[11]使用水平集和流體體積耦合方法（CLVOF）模擬在平板加熱面上膜態(tài)沸騰換熱過程。近年來,Lattice-Boltzmann（LB）方法在多相流體的模擬中引起了極大的關(guān)注。Sadeghi 等[12]建立了三維LBM 模型模擬水平放置平板加熱面膜態(tài)沸騰。Ma 和Cheng[13]比較了從二維和三維LB 模型,發(fā)現(xiàn)對于穩(wěn)定的膜態(tài)沸騰,計(jì)算結(jié)果幾乎相同。

本研究采用數(shù)值計(jì)算方法研究不同重力水平下液氫膜態(tài)沸騰特性與規(guī)律,以揭示微重力條件下氫的氣泡運(yùn)動行為和沸騰傳熱強(qiáng)度。研究工作不僅為氫在空間中的安全高效應(yīng)用提供理論指導(dǎo),而且對兩相沸騰換熱有關(guān)的科學(xué)技術(shù)也具有意義。

2 數(shù)值模型

2.1 研究對象

針對液氫膜態(tài)沸騰,本研究建立了如圖1 所示的計(jì)算模型,細(xì)加熱絲浸泡于過冷液氫池中,液池上方為氫氣枕,幾何參數(shù)及初始條件如表1 所示。首先對ΔTsub=2 K,ΔTsup=30 K 在常重力工況進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證,結(jié)果展示于表2 中,最終選取120 萬網(wǎng)格的模型進(jìn)行計(jì)算。

表1 幾何參數(shù)及初始條件Table 1 Geometric parameters and initial conditions

表2 常重力工況下網(wǎng)格無關(guān)性驗(yàn)證Table 2 Mesh independence analysis

圖1 計(jì)算模型示意圖Fig.1 Schematic diagram of calculation model

2.2 計(jì)算設(shè)置

采用CFD 數(shù)值計(jì)算軟件Fluent 19.2 進(jìn)行計(jì)算,基于壓力的瞬態(tài)求解器,使用VOF 多相模型捕獲氣液相界面運(yùn)動情況,并選擇SIMPLE 算法求解壓力和速度方程。密度、動量和能量方程通過二階迎風(fēng)格式離散化,體積分?jǐn)?shù)通過Geo-Reconstruct 格式離散化。梯度和壓力分別通過基于Green-Gauss Cell 和PRESTO 格式離散化。氫的物性參考流體特性數(shù)據(jù)庫（REFPROP）,并將氫氣的密度視為理想氣體。時(shí)間步長從10-6s 變化到10-4s 以滿足收斂要求,并且設(shè)置收斂標(biāo)準(zhǔn)將能量方程的殘差控制在10-6以下,其他項(xiàng)控制方程殘差在10-4以下。其中連續(xù)方程、動量方程和能量方程如下:

式中:α為體積分?jǐn)?shù),ρ為密度,S為體積源項(xiàng),v為速度項(xiàng),E為內(nèi)能項(xiàng),P為壓力項(xiàng),Q為能量轉(zhuǎn)移速率。

每一相的體積分?jǐn)?shù)α取值在0 到1 之間,且加和為1。

在模型中采用連續(xù)表面力（CSF）模型計(jì)算表面張力對流體界面運(yùn)動的作用,該模型將表面張力轉(zhuǎn)化為體積力源項(xiàng)Fvol加載至動量方程進(jìn)行求解計(jì)算。

其中,界面曲率κ可由以下公式計(jì)算。

式中:κl和κg為表面曲率,由界面上表面法線的局部梯度計(jì)算而來。

2.3 相變模型

選取Lee 模型計(jì)算氣液相間的質(zhì)量、能量轉(zhuǎn)移作用,并通過UDF 的形式,向連續(xù)方程與能量方程加載對應(yīng)的源項(xiàng)。數(shù)值計(jì)算中,相變源項(xiàng)將作用于所有流體網(wǎng)格,通過比較流體網(wǎng)格與飽和溫度的相對大小確定相變量,其判斷依據(jù)如下:

當(dāng)T≥Tsat時(shí)（液相蒸發(fā)）

當(dāng)T＜Tsat時(shí)（氣相冷凝）

伴隨著質(zhì)量的傳遞,能量轉(zhuǎn)移速率可按下式確定:

式中:Ci表示相變系數(shù),1/s;用來調(diào)節(jié)每次迭代的相變量大小。Ci的取值應(yīng)找到一個(gè)合適的范圍,偏小會出現(xiàn)計(jì)算誤差,而偏大則會導(dǎo)致計(jì)算發(fā)散。其取值與幾何尺寸、網(wǎng)格大小及時(shí)間步長等很多因素有關(guān)。圖2 展示了壁面面積熱流量隨相變系數(shù)的變化情況,由圖可知,計(jì)算中Ci合適取值范圍在103—104之間。

圖2 相變系數(shù)驗(yàn)證Fig.2 Verification of phase transition coefficient

3 模型驗(yàn)證

利用文獻(xiàn)[13]中的實(shí)驗(yàn)數(shù)據(jù)對模型進(jìn)行驗(yàn)證,選取1 g 和0.235 g 兩個(gè)重力工況,實(shí)驗(yàn)數(shù)據(jù)和模擬結(jié)果比較如圖3 所示,兩者吻合較好。因此,本研究建立模型可用于不同重力下氫膜態(tài)沸騰換熱預(yù)測研究。

圖3 數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)比較Fig.3 Comparison between numerical simulation results and experimental data

4 結(jié)果分析

基于上述模型,首先對常重力工況下加熱細(xì)絲浸沒在過冷液氫池膜態(tài)沸騰過程中氣泡運(yùn)動行為及換熱特性進(jìn)行揭示。在此基礎(chǔ)上通過分析氣泡脫離直徑、氣泡生長周期及壁面熱流,進(jìn)一步研究了重力水平和液體過冷度對膜態(tài)沸騰換熱特性的影響。

4.1 氣泡行為與換熱特性

圖4 展示了在ΔTsub=2 K,ΔTsup=30 K 和1 g 條件下加熱細(xì)絲浸沒在過冷液氫池中發(fā)生膜態(tài)沸騰過程氣泡的形成、生長及脫落過程。首先加熱絲表面形成一層氣膜,在浮力作用下氣體在細(xì)絲上方積聚,隨著氣體量的增多,在細(xì)絲頂部可依次觀察到氣泡的形成及脫落過程。氣泡的運(yùn)動行為進(jìn)一步影響了壁面換熱特性,圖5 為穩(wěn)定膜態(tài)沸騰過程中壁面熱流隨時(shí)間的變化情況。由圖可知,在氣泡連續(xù)形成與脫落過程中,壁面熱流隨之經(jīng)歷了周期性變化。隨著氣膜的增厚,由于氣體熱阻增大壁面熱流隨之降低,當(dāng)氣體積聚在細(xì)絲上方開始形成氣泡時(shí),加熱絲下方氣膜厚度減小,熱流隨之增加,直至氣泡脫離。

圖4 膜態(tài)沸騰過程氣泡運(yùn)動行為FIg.4 Bubble movement behavior in film boiling process

圖5 膜態(tài)沸騰過程熱流隨時(shí)間變化Fig.5 Variation of heat flow in film boiling process with time

4.2 重力水平影響分析

圖6 展示了在ΔTsub=2 K,ΔTsup=30 K 和0.5 g條件下,含氣率、熱流和氣泡行為隨時(shí)間的變化。如圖所示,在此重力水平下,氣泡的行為與常重力下的相同,并且含氣率和熱流隨時(shí)間規(guī)律的周期波動。在相界面處,蒸發(fā)和冷凝的兩種作用同時(shí)發(fā)生,如果蒸發(fā)作用大于冷凝作用,則含氣率增加。相反,如果冷凝效果更強(qiáng),則含氣率降低。蒸發(fā)作用的能量來自于氣膜的熱傳導(dǎo),冷凝作用是由過冷液體和飽和溫度之間的溫度差驅(qū)動。由于液池較大,該溫度差幾乎保持恒定,因此,含氣率主要取決于加熱絲傳導(dǎo)至氣液界面的熱量。當(dāng)氣膜增厚時(shí),換熱量大大降低,進(jìn)而導(dǎo)致液池中含氣率降低。隨著氣體在浮力作用下積聚在上部形成氣泡并脫落時(shí),加熱絲下方的氣膜變薄,換熱量與含氣率開始增加,直到氣泡從界面釋放。與正常重力相比,在0.5 g 的條件下,氣泡的生長時(shí)間變長,平均換熱量變小。

圖6 0.5 g 工況下膜態(tài)沸騰過程熱流和含氣率隨時(shí)間變化Fig.6 Variation of heat flow and gas holdup in film boiling process with time at 0.5 g

隨著重力水平的進(jìn)一步降低,浮力作用更加減弱,在0.1 g 重力條件下熱流、含氣率和相界面隨時(shí)間變化情況如圖7 所示。由圖可知,在此重力水平下,膜態(tài)沸騰過程中沒有氣泡的形成和脫離,熱流和含氣率隨著氣膜在加熱絲表面的震蕩而波動。當(dāng)加熱絲被厚的氣膜包裹時(shí),換熱量降低到最小值,之后在冷凝作用下含氣率降低,一旦含氣率降至一定程度時(shí),換熱量立即增加,導(dǎo)致加熱絲再次被厚氣膜包裹。因此,在0.1 g 的重力水平和2 K 的過冷度下,沒有氣泡的生成與脫離現(xiàn)象。

圖7 0.1 g 工況下膜態(tài)沸騰過程熱流和含氣率隨時(shí)間變化Fig.7 Variation of heat flow and gas holdup in film boiling process with time at 0.1 g

為揭示加熱絲表面的過冷膜態(tài)沸騰換熱特性,圖8 展示了氣泡脫離直徑、氣泡生長時(shí)間和換熱量隨重力水平（從0.01 g 到1 g）的變化情況。隨著重力水平的降低,氣泡脫離直徑和氣泡生長時(shí)間明顯增加,流體與壁面間的換熱量相應(yīng)降低。一旦重力水平降至一定程度,氣膜在加熱絲表面不停晃動,不會形成氣泡。這是由于液體存在過冷度,氣膜不會無限增厚,而在重力水平較低的時(shí)候,不會產(chǎn)生足夠的氣體形成預(yù)期的大氣泡。對于ΔTsub=2 K,ΔTsup=30 K的計(jì)算工況,臨界重力在0.1 g 至0.15 g 之間。

圖8 不同重力水平下熱流與氣泡脫離直徑和生長時(shí)間變化情況Fig.8 Change of heat flow and bubble separation diameter and growth time under different gravity levels

4.3 液體過冷度影響分析

為了分析液體過冷度的影響作用,圖9 展示了在0.1 g 重力水平氣泡脫離直徑、氣泡生長時(shí)間和換熱量隨液體過冷度的變化情況。隨著液體過冷度的增加,氣泡的生長時(shí)間明顯增加,氣泡脫離直徑幾乎保持相同。根據(jù)圖6 和圖7 可知,在穩(wěn)定膜態(tài)沸騰過程中,氣體含氣率在一定范圍周期性波動。圖10 展示出了在不同液體過冷溫度下液池內(nèi)含氣率的時(shí)間平均值變化情況。含氣率隨過冷度的增加而減小,當(dāng)氣體減少到一定程度時(shí),就不會形成氣泡。這是因?yàn)橐后w過冷度決定了冷凝效果,當(dāng)冷凝作用足夠強(qiáng)時(shí),不能產(chǎn)生足夠的氣體在0.1 g 重力水平下形成氣泡。結(jié)合重力水平的影響,可知在低重力水平且高液體過冷度工況中,更容易發(fā)生氣膜附著在加熱絲上晃動,而沒有氣泡產(chǎn)生脫落的現(xiàn)象。

圖9 0.1 g 工況不同液體過冷度下熱流與氣泡脫離直徑和生長時(shí)間變化情況Fig.9 Change of heat flow and bubble separation diameter and growth time under different liquid undercooling at 0.1 g

圖10 0.1 g 工況不同液體過冷度下含氣率變化情況Fig.10 Variation of gas holdup under different liquid undercooling at 0.1 g

5 結(jié) 論

為揭示微重力工況下液氫膜態(tài)沸騰機(jī)理,本研究建立了加熱細(xì)絲浸沒在過冷液氫池中的數(shù)值計(jì)算模型,從氣泡運(yùn)動行為和換熱特性兩個(gè)方面開展研究,主要結(jié)論如下:

（1）液體過冷度和重力水平是影響液氫過冷膜態(tài)沸騰傳熱特性的兩個(gè)重要因素,其中液體過冷度決定了液池中的氣體量,而重力水平則影響著氣泡的大小以及氣泡脫離氣液界面的能力。

（2）針對高重力水平、低液體過冷度的工況,氣泡在加熱絲上方不斷產(chǎn)生與脫落。隨著重力水平降低,氣泡脫離直徑和生長時(shí)間增加,流體與壁面間的換熱量隨之降低。

（3）針對低重力水平、高液體過冷度的工況,氣膜附著在加熱絲表面不斷晃動,而不會發(fā)生氣泡形成與脫落現(xiàn)象。