王 夕，王 玨，容 易，黃 輝

?

微重力下低溫貯箱內(nèi)推進劑相變仿真模型研究

王 夕1，王 玨2，容 易1，黃 輝1

（1. 北京宇航系統(tǒng)工程研究所，北京，100076；2.中國運載火箭技術(shù)研究院，北京，100076）

低溫推進劑具有沸點低、易汽化的特點，相變是低溫推進劑長時間在軌蒸發(fā)量控制問題中需要考慮的首要影響因素。相變模型對低溫推進劑蒸發(fā)仿真起到重要作用，構(gòu)建合理的相變模型成為低溫流體蒸發(fā)量仿真重要的研究方向?；?種相變理論，采用FLUENT軟件二次開發(fā)的方法，建立基于相平衡和非平衡理論的4種相變模型，開展微重力下液氫推進劑蒸發(fā)的數(shù)值模擬，并與國外探空火箭試驗進行比較和驗證。研究結(jié)果表明：比較4種相變模型對貯箱內(nèi)壓力升高速率預(yù)示的準確性，得出了適用于微重力下低溫推進劑仿真的相變模型。

低溫推進劑；相變；微重力

0 引 言

在液體運載火箭的推進劑中，低溫推進劑具有比沖高、沸點低、難于貯存的特點，因此限制了其長時間的在軌使用。對低溫推進劑長時間在軌蒸發(fā)量控制問題的研究，可以采用試驗和仿真的方法。鑒于低溫和微重力試驗的難度，基于已有商業(yè)軟件的數(shù)值仿真是對貯箱內(nèi)低溫推進劑開展研究的有效途徑。

目前，低溫推進劑蒸發(fā)的數(shù)值研究受到廣泛的關(guān)注。在低溫推進劑蒸發(fā)量的所有影響因素中，相變是重要的考慮因素，相變過程的建模是蒸發(fā)量控制研究的焦點問題，而發(fā)展相變數(shù)值模型的基礎(chǔ)是描述相變物理過程的理論模型。相變的理論模型包括相平衡模型和非平衡模型。相平衡模型中，假設(shè)液相和汽相之間達到平衡時，界面各相之間的溫度相等，處于平衡狀態(tài)[1]；而在非平衡模型中，Schrage[2]使用基于Maxwell速度分布來分析液/汽界面的傳質(zhì)過程，將界面的相變考慮為兩相界面對液/汽分子的捕獲和逃逸的非平衡過程；李震東[3]、Tanasawa[4]等總結(jié)了發(fā)生在液/汽相界面處的相變傳質(zhì)理論，歸納為相平衡模型和非平衡模型；劉秋生[5]等使用簡化的Hertz-Knudsen方程對熱毛細作用主導(dǎo)的流動進行了研究，對微重力下的相變問題具備參考價值。

基于以上兩類相變理論，建立了相應(yīng)的仿真模型，并用于低溫貯箱推進劑蒸發(fā)的仿真中。其中，相平衡模型在低溫貯箱CFD仿真中應(yīng)用廣泛。程向華[6]等根據(jù)相平衡建立了液氧貯箱的二維模型，分析了液氧熱分層的形成過程及原因；Zilliac[7]等對貯箱進行建模，根據(jù)熱力學(xué)平衡原理建立相變模型。相比于相平衡模型在低溫推進劑相變仿真的廣泛使用，非平衡模型的使用很少。Olga[8]等應(yīng)用非平衡模型開展了液氫貯箱內(nèi)徑向噴霧的仿真，但并沒有將非平衡模型與相平衡模型進行對比研究。

實際的相變物理過程在界面處存在溫度跳變和非平衡的瞬態(tài)作用，與相平衡模型存在一定的差異，這些差異對低溫推進劑蒸發(fā)數(shù)值模擬仿真準確性的影響仍未被討論，比較不同相變仿真模型的準確性和適用性對低溫流體仿真具有重要作用。

本文基于已有的相變理論展開了CFD建模工作，通過FLUENT用戶自定義函數(shù)，建立基于不同相變理論的蒸發(fā)模型，并通過適當(dāng)?shù)募僭O(shè)修正蒸發(fā)模型，開展了貯箱內(nèi)液氫蒸發(fā)仿真，與國外Areobee探空火箭試驗[9]進行了對比，并分析了不同的相變模型在低溫流體仿真中的應(yīng)用。

1 數(shù)值模型

本文參照Aerobee探空火箭液氫蒸發(fā)試驗[9]開展了數(shù)值模擬的對比研究，基于相平衡理論和3種非平衡理論，通過界面假設(shè)和蒸發(fā)假設(shè)，建立了仿真軟件FLUENT下的相變模型，模擬試驗貯箱內(nèi)壓力上升，驗證并比較不同模型的結(jié)果，得到適用于微重力下低溫推進劑蒸發(fā)量預(yù)示的數(shù)值仿真模型。

1.1 控制方程

采用二維軸對稱模型，對液氫貯箱推進劑受熱蒸發(fā)的過程進行分析?？刂品匠虨檫B續(xù)方程、N-S方程和能量方程，對于每一相，有：

式中為密度；為速度；為動力粘度；為時間；為壓力；為能量；為系統(tǒng)加速度；為溫度；為體積力；為能量源項。

液相密度采用Boussinesq近似，其余物性作為溫度的函數(shù)分段插值，氣相則采用理想氣體模型。

（4）

考慮到蒸發(fā)過程流速小，因此流動采用層流模型。

微重力下需要考慮表面張力作用，采用連續(xù)表面張力模型如下：

1.2 相變模型

相變模型是本文主要研究對象，基于兩種假設(shè)，開展對4種相變模型的研究。

a）界面假設(shè)：相變僅發(fā)生在相界面處；

b）蒸發(fā)假設(shè)：只有在液相溫度高于飽和溫度時，才發(fā)生相變，不考慮冷凝的發(fā)生。

1.2.1 模型1

考慮使用相平衡模型[6]，液相的超過飽和溫度的能量將全部轉(zhuǎn)化為相變，模型1方程為

1.2.2 模型2

Schrage[2]推薦的非平衡模型，即Hertz-Knudsen方程為

凈質(zhì)量流率為正值，表明氣液界面上有凈的質(zhì)量從氣相變?yōu)橐合?，即相界面上發(fā)生凝結(jié)；反之，如果凈質(zhì)量流率為負值，則相界面上將發(fā)生蒸發(fā)。

1.2.3 模型3

劉秋生[5]等提出采用非平衡熱力學(xué)方法計算表面張力主導(dǎo)流動下的界面蒸發(fā)流量。蒸發(fā)界面假設(shè)為不變形，采用線性近似后的Hertz-Knudsen方程來描述蒸發(fā)界面的蒸發(fā)流量，即模型3：

1.2.4 模型4

采用相變模型[4,10]，即模型4為

根據(jù)Tanasawa[4]的推薦值及實際計算的結(jié)果比較，取系數(shù)為0.2。

1.3 網(wǎng)格無關(guān)性驗證

對采用的二維軸對稱網(wǎng)格（見圖1）開展網(wǎng)格無關(guān)性驗證（見圖2），使用相變模型1，比較3種不同網(wǎng)格數(shù)量，網(wǎng)格數(shù)分別為5000，10 000和20 000。

圖1 采用的網(wǎng)格

圖2 網(wǎng)格無關(guān)性驗證曲線

3種網(wǎng)格壓力上升的偏差較小，其中與較高網(wǎng)格數(shù)量相比，中等網(wǎng)格數(shù)量的最終壓力差值占總壓力升高的4.17%，該偏差可以被忽略，結(jié)果較為一致?？紤]到網(wǎng)格數(shù)的增加對計算效率的影響，采用中等網(wǎng)格數(shù)量開展數(shù)值仿真。

2 結(jié)果與分析

對4種相變模型開展研究，對比不同模型下貯箱內(nèi)壓力上升速率與試驗實測值的偏差，同時對貯箱內(nèi)的流體行為開展研究。仿真條件如表1所示，選取與試驗相同的初始壓力，并將初始溫度設(shè)置為該壓力下飽和溫度。假設(shè)貯箱為均勻受熱，選取試驗中平均熱流密度作為壁面熱流邊界條件。

表1 數(shù)值模擬條件

Tab.1 Conditions of Simulation

性能參數(shù) 球形貯罐直徑/cm22.86 體積充填率34.3% 平均熱流密度/（W·m-2）473.19 重力/（m·s-2）0.01 初始壓力/MPa0.1241 總受熱時間/s237

由于采用了兩種假設(shè)，本文將對界面假設(shè)和蒸發(fā)假設(shè)的合理性展開討論，并基于兩種假設(shè)，對采用的4種相變模型的對比分析。

2.1 界面假設(shè)

界面假設(shè)將相變界定在界面處，是否采用界面假設(shè)將對結(jié)果產(chǎn)生影響。采用/不采用界面假設(shè)代表兩種不同的相變原理。

汽化過程可分為界面蒸發(fā)和沸騰蒸發(fā)兩種方式。當(dāng)模型不采用界面假設(shè)，液相溫度只要溫度超過飽和溫度，在任何位置均可發(fā)生汽化，這樣的模型更接近于沸騰。然而，在實際過程中，液體內(nèi)部發(fā)生沸騰并產(chǎn)生氣泡需要達到一定的活化能，并非達到飽和溫度即開始汽化，在較低熱流密度的條件下相變往往以界面蒸發(fā)為主。熱量傳遞到界面需要一定的時間，當(dāng)界面蒸發(fā)主導(dǎo)相變時，不采用界面假設(shè)將高估汽化傳質(zhì)速率，即使與實際的沸騰過程相比，由于未考慮氣泡生成的過熱度，不采用加密假設(shè)的模型同樣會高估汽化傳質(zhì)速率。

采用模型1和模型4來比較有/無界面假設(shè)的貯箱壓力上升結(jié)果如表2所示。在237 s內(nèi)，無界面假設(shè)的模型得到的貯箱壓力和平均壓力升高速率均高于試驗值，而采用界面假設(shè)的模型得到的結(jié)果更接近于試驗值。結(jié)果表明，不采用界面假設(shè)的模型將高估汽化傳質(zhì)速率，使壓力快速升高；采用界面假設(shè)的模型更為合理，進一步支撐了假設(shè)的合理性。對照試驗以界面蒸發(fā)過程為主，與沸騰過程有本質(zhì)的區(qū)別，汽化過程僅發(fā)生在界面處，采用界面假設(shè)的模型也更符合實際。

表2 界面假設(shè)壓力上升結(jié)果對比

Tab.2 Pressure Rise Result of Interface Assumption

模型終壓/MPa壓力升高速率/(Pa·s-1) 試驗0.7302752557.7 模型1（無界面假設(shè)）0.8424493031.0 模型4（無界面假設(shè)）0.8882393224.2 模型1（界面假設(shè)）0.7569942670.4 模型4（界面假設(shè)）0.7682602718.0

在低溫推進劑貯箱的相變仿真中，是否采用界面假設(shè)，取決于流體的形態(tài)。若貯箱受熱較小，主要形態(tài)為界面蒸發(fā)，則使用界面假設(shè)；若貯箱受熱較大，主要形態(tài)為沸騰過程，則應(yīng)當(dāng)考慮不采用界面假設(shè)。

2.2 蒸發(fā)假設(shè)

蒸發(fā)假設(shè)的實質(zhì)是忽略冷凝的作用。在實際過程中，發(fā)生在氣相內(nèi)部（非界面處）的均勻冷凝一般很難發(fā)生，而非均勻冷凝則發(fā)生在有過冷界面存在的情形下[2]。而試驗中不存在大過冷度的過冷界面，若認為只要低于飽和溫度即發(fā)生冷凝，則可能高估冷凝的作用，造成計算偏差。因此，不采用純蒸發(fā)假設(shè)可能高估冷凝的作用。

采用模型1和模型3來比較有/無蒸發(fā)假設(shè)的貯箱壓力上升結(jié)果如圖3所示。

圖3 蒸發(fā)假設(shè)壓力上升結(jié)果對比曲線

圖3中，最終的壓力值和平均壓力上升速率可用于比較分析。由于試驗中加熱器啟動需要一定的時間，溫度由開始加熱到逐漸穩(wěn)定，因此在試驗初始階段，加熱的熱流更小。而仿真中，熱流邊界條件設(shè)置為平均熱流密度，并保持恒定值，因而圖中曲線在100 s處壓力比數(shù)值模擬結(jié)果低是合理的，最終壓力是較為合理的比較參數(shù)。

由圖3可知，不采用蒸發(fā)假設(shè)時，模型3的最終壓力結(jié)果遠低于試驗值，而模型1的最終壓力同樣低于試驗值。不采用純蒸發(fā)假設(shè)的模型，高估了冷凝的作用，造成凈汽化速率的降低，壓力上升速率減緩，導(dǎo)致了與試驗的偏離。

采用蒸發(fā)假設(shè)的模型1和模型3的最終壓力值與試驗值更接近，結(jié)果更準確。通過有/無蒸發(fā)假設(shè)的對比分析，在沒有過冷壁面存在的條件下，采用蒸發(fā)假設(shè)的相變傳質(zhì)模型更合理。

2.3 4種相變模型的結(jié)果分析

在界面假設(shè)和蒸發(fā)假設(shè)下，對4種蒸發(fā)模型開展了對比研究，以分析模型的準確性以及相平衡模型和非平衡模型的差別。

4種模型結(jié)果對比如表3所示，4種模型的壓力結(jié)果如圖4所示。

表3 4種模型結(jié)果對比

Tab.3 Results of Four Models

模型最終壓力/MPa平均壓力升高速率/（Pa·s-1）偏差 試驗0.730 2752557.7— 模型10.756 9942670.44.4% 模型20.657 5342250.8-12% 模型30.734 5462575.70.7% 模型40.768 2602718.06.3%

圖4 4種模型壓力結(jié)果

由表3可知，模型2的平均壓力升高速率偏差最大，為-12.0%。由于在相變模型中考慮了壓力，而壓力是較為敏感的參數(shù)，同時考慮壓力和溫度也增加了模型的復(fù)雜性，其計算結(jié)果可能出現(xiàn)較大的偏差。

模型4的壓力升高速率結(jié)果偏差為6.3%，比模型2的偏差更小，比其他兩種模型略為偏大，略高于試驗貯箱壓力結(jié)果。模型1是相平衡模型，而其他3種模型均為非平衡模型。模型1的偏差為4.4%，結(jié)果優(yōu)于模型2、模型4，略高于試驗貯箱壓力結(jié)果。

模型3的壓力升高速率與試驗的偏差最小，偏差為0.7%，最適合用于微重力下液氫貯箱蒸發(fā)量的預(yù)示。

通過4種模型的對比和分析發(fā)現(xiàn)，模型1、模型3及模型4的偏差均小于7%，對微重力下低溫推進劑貯箱蒸發(fā)量預(yù)示的偏差均可被接受。其中，模型3考慮了在熱毛細力主導(dǎo)下的界面蒸發(fā)，對微重力下的情形更為適用，其結(jié)果相比于其他3種模型具備明顯的優(yōu)勢，偏差最小，適用于微重力下低溫推進劑貯箱蒸發(fā)量預(yù)示和對貯箱壓力升高的數(shù)值模擬，具備良好的準確性。

在國內(nèi)外對低溫推進劑蒸發(fā)的研究中，使用非平衡相變模型的較少，一般使用相平衡模型對蒸發(fā)過程進行數(shù)值模擬。根據(jù)以上的分析，相平衡模型具備在趨勢上預(yù)測貯箱壓力的能力，其仿真結(jié)果偏差在可接受的范圍內(nèi)，這些研究的結(jié)果并不會在趨勢上受蒸發(fā)模型的影響。

在微重力這樣的特殊環(huán)境中，模型3相比于相平衡模型和其他非平衡模型具備準確性優(yōu)勢。因此，在特定環(huán)境下，采用與環(huán)境相適應(yīng)的非平衡模型將進一步提高仿真的準確性。

3 結(jié) 論

本文利用國外Aerobee探空火箭試驗數(shù)據(jù)對微重力下液氫貯箱內(nèi)相變仿真開展CFD建模工作，根據(jù)不同的相變理論，通過兩種假設(shè)，對不同相變模型開展比較分析，得到如下結(jié)論：

a）在較低熱流、貯箱內(nèi)以界面蒸發(fā)為主導(dǎo)的情況下，宜采用界面假設(shè)；在沒有過冷壁面存在的條件下，對低溫推進劑受熱蒸發(fā)相變的數(shù)值模擬宜采用蒸發(fā)假設(shè)。

b）本文比較了4種相變模型對封閉液氫貯箱內(nèi)壓力升高的預(yù)測，其中包括1種相平衡模型和3種非平衡模型，結(jié)果表明，模型1、3、4均能在趨勢上對微重力下液氫封閉貯箱壓力升高趨勢進行預(yù)測，可有效模擬貯箱中的流體形態(tài)和溫度分布。相平衡模型1具備對貯箱內(nèi)壓力升高趨勢的預(yù)測能力。

c）非平衡模型3仿真結(jié)果更接近試驗，對微重力下液氫蒸發(fā)的數(shù)值模擬具備良好的準確性，在微重力環(huán)境下的仿真中具備優(yōu)勢。在特定環(huán)境下，采用與環(huán)境相適應(yīng)的非平衡模型將進一步提高仿真的準確性。

[1] 吳雙應(yīng), 曾丹苓, 黃雙, 等. 微層蒸發(fā)模型的非平衡熱力學(xué)分析[J]. 熱科學(xué)與技術(shù), 2005, 4(4): 283-288.

Wu Shuangying, Zeng Danling, Huang Shuang, et al. Non-equilibrium thermodynamic analysis of micro-layer evaporation model[J]. Journal of Thermal Science and Technology, 2005, 4(4): 283-288.

[2] Schrage R W. A theoretical study of interphase mass transfer[M]. NewYork: Columbia University Press, 1953.

[3] 李震東, 趙建福, 魯仰輝, 等. 池沸騰現(xiàn)象中熱毛細對流的成因[J]. 空間科學(xué)學(xué)報, 2008, 28(1): 38-43.

LI Zhendong, Zhao Jianfu, Lu Yanghui, et al. Origin of thermocapillary convection in pool boiling[J]. Chinese Journal of Space Science, 2008, 28(1):38-43.

[4] Tanasawa I. Advances in condensation heat transfer[J]. Advances in Heat Transfer, 1991, 21: 55-139.

[5] 劉秋生, 汪洋, 紀巖. 蒸發(fā)相變與界面流動耦合機理研究[J]. 工程熱物理學(xué)報, 2010(10): 1751-1754.

Liu Qiusheng, Wang Yang, Ji Yan. Coupling mechanism of evaporation phase-change and interfacial flow[J]. Journal of Engineering Thermophysics, 2010(10): 1751-1754.

[6] 程向華, 厲彥忠, 陳二鋒, 等. 新型運載火箭射前預(yù)冷液氧貯箱熱分層的數(shù)值研究[J]. 西安交通大學(xué)學(xué)報, 2008, 42(9): 1132-1136.

Cheng Xianghua, Li Yanzhong, Chen Erfeng, et al. Numerical investigation of thermal stratification in liquid Oxygen tank for new-style launch vehicle during ground precooling[J]. Journal of Xi'an Jiaotong University, 2008, 42(9): 1132-1136.

[7] Zilliac G, Karabeyoglu M A. Modeling of propellant tank pressurization[R]. AIAA, 2005-3549, 2005.

[8] Kartuzova O, Kassemi M, Agui J H, et al. A CFD model for the Multipurpose Hydrogen Test Bed (MHTB) ground-based self-pressurization and pressure control experiments[C]. Atlanta: 11th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, 2014.

[9] Knoll R H, Smolak G R. Weightlessness experiments with liquid hydrogen in aerobee sounding rockets, uniform radiant heat addition-flight 1[R]. NASA TM X-484, 1962.

[10] ANSYS FLUENT 13.0 Theory Guide[M]. Canonsburg, PA, ANSYS. Inc., 2011.

Computational Research on Phase Change Model forCryogenic Propellant in Microgravity

Wang Xi1, Wang Jue2, Rong Yi1, Huang Hui1

(1. Beijing Institute of Astronautical Systems Engineering, Beijing, 100076; 2. China Academy of Launch Vehicle Technology, Beijing, 100076)

Boil off is a characteristic of cryogenic propellant, and the phenomenen of phase change is the major factor of cryogenic propellant storage on orbit. As the phase change model plays a key role in cryogenic fluid Computational Fluid Dynamics(CFD) simulation, the usage of reasonable phase change model will be important. Based on four different phase change theory, the commercial software FLUENT are used and four different phase change models are presented by user defined secondary developing code. A CFD research on liquid hydrogen evaporation in microgravity is present, comparing with NASA sounding rocket experiment. Finally, a reasonable phase change model is commended by comparing the pressurizing rate of the four models, which can be used in cryogenic fluid management simulation.

Cryogenic propellant; Phase change; Microgravity

1004-7182(2018)01-0036-05

10.7654/j.issn.1004-7182.20180107

V511

A

2016-12-14；

2017-02-20

王 夕（1989-），男，博士，工程師，主要研究方向為運載火箭總體設(shè)計