水升華器升華模式理論分析與數值仿真

廖俊元楊春信楊涵

(北京航空航天大學航空科學與工程學院，北京 100191)

1 引言

隨著中國航天技術的發展及載人登月工程的逐步推進，航天員的艙外活動需求逐漸提高。艙外航天服的熱控系統需要保證航天員在出艙活動時處于安全、舒適的熱環境之中。水升華器是一種利用水的升華潛熱進行散熱的消耗型相變散熱裝置，相比于輻射器等熱控器件具有體積小、質量輕、效率高、微重力下工作可靠等優點，是一種性能優良的空間熱沉，十分適用于應對小散熱面、短時大功耗任務的艙外航天服的熱控系統。

國際對水升華器的研究開展較早，從上世紀60年代初已有關于水升華器作為熱沉應用于航天器熱控系統的相關研究，并制作了如LM-107、LM-209等一系列原理樣機進行了初步的實驗探究，如圖1所示。在上世紀70年代初至上世紀末的期間，水升華器處于大量工程應用的階段，在首次成功應用于美國的阿波羅探月工程中后，還在蘇聯(俄羅斯)的Orlan航天服等航天工程中得到應用。進入新世紀后，為了滿足更高的航天工程需求，水升華器出現了追求防污染、良好動態性能的發展趨勢，對此NASA的研究人員提出了X-38、CIS、SDC、ISDC 4種水升華器的新構型，如圖2所示。

圖1 早期水升華器Fig.1 Early application of sublimator

圖2 新型水升華器Fig.2 New type of sublimator

中國對水升華器的相關研究起步較晚。吳志強等最早對水升華器進行了散熱性能的理論分析及一系列實驗研究；李森等建立了多孔板單孔內水相變過程的物理數學模型，并進行了周期模式的數值模擬；王玉瑩等采用一維熱阻網格模型，通過仿真得到了周期模式下界面位置、水升華器溫度等參數的變化，并對不同加熱邊界的水升華器工作性能進行了實驗研究。航天工程方面，中國關于水升華器的首次成功應用為2008年神舟七號飛船出艙活動的“飛天”艙外航天服。

由于水的蒸發、升華等相變過程發生在水升華器內部的給水腔、多孔板微孔等尺寸微小的區域，實驗僅能測得外部溫度等宏觀數據，因此對于水升華器工作原理的研究主要通過理論分析、數值仿真來進行。但是，目前國內對水升華器的仿真分析較少且主要為一維數值仿真，沒有考慮水升華器不同工作模式的轉換。為了探究水升華器工作模式轉換的規律，分析不同設計參數對水升華器工作特性的影響，本文在總結前人理論分析的基礎上，對水升華器工作模式轉換特性進行量綱分析，通過Fluent數值仿真的方法對水升華器在低熱載荷下的穩態工作模式(升華模式)開展研究，分析不同多孔板結構參數、熱載荷大小下水升華器工作特性的變化規律。

2 水升華器工作過程分析

2.1 水升華器工作原理

本文研究對象為平板型水升華器，其基本結構由供水入口、加熱板、給水腔、多孔板組成，如圖3所示。

圖3 升華模式下的水升華器Fig.3 Water sublimator working in sublimation mode

系統工作過程為:工質水在一定的供水壓力下被送往高真空度的給水腔，水的壓力驟降至其三相點以下而迅速蒸發。蒸發潛熱帶走大量熱量，若此時水溫降至三相點溫度以下，則會凝固形成固態冰。外表面的冰層暴露于高真空環境下，熱載荷將通過升華潛熱排向外太空。

Hamilton Standard公司根據工質水的不同消耗形式，將水升華器工作模式劃分為4種:升華模式、蒸發模式、混合模式、周期模式。當熱載荷較小時，水升華器工作在升華模式。由于水蒸氣流量較小，升華界面對應的飽和蒸氣壓低于水的三相點壓力。根據圖4，此時水不會進入多孔板的毛細孔，而是首先在多孔板與給水腔的界面上形成一定厚度的冰層(圖3)。此時，供水流量等于升華流量，升華帶走的熱量等于底面加熱及供水熱量之和，水升華器內部達到整體的熱質動態平衡。升華模式是水升華器最期望的工作模式[12]。

圖4 水的三相圖Fig.4 Three-phase diagram of water

當熱載荷繼續增大，水將進入多孔板形成一定厚度的冰層。冰層在升華、熔化的共同作用下逐漸耗盡后，水繼續流入并凝固成冰層，形成了周期性的循環工作過程，即周期模式。若多孔板溫度較高，則冰層不會出現，工質將以液態水蒸發的形式進行散熱，即蒸發模式。水升華器實際工作中可能會出現以上3種模式共存的情況，這種模式被稱為混合模式。

2.2 升華模式分析

吳志強等提出升華模式下水升華器整體處于穩態工作的狀態，可以通過理論分析計算升華熱流、加熱面溫度等參數。假設如下:

1)多孔板、加熱板及工質水的物性為常數；

2)不考慮結構漏熱、輻射漏熱等熱損失；

3)多孔板毛細孔簡化為眾多同徑的豎圓孔；

4)忽略給水腔內部水對流造成的換熱；

5)加熱板、給水腔、多孔板具有相同的面積。

升華模式下水升華器內部的熱質平衡如圖3中展示，圖中紅色虛線箭頭代表熱量的流動，藍色實線箭頭代表質量的流動。定義加熱板熱流密度如式(1):

式中，

q

為底面加熱熱流，W/m；

Q

為水升華器熱載荷，W；

A

為加熱板面積，m。進一步對液態水、冰層分別進行熱平衡分析，可得式(2)、(3):

式中，

q

、

q

分別為水層、冰層的導熱熱流，W/m；

m

為供水流量，kg/s，在升華模式下等于總體升華流量；

t

、

t

分別為供水溫度、水-冰相變溫度，K；

c

為水的比熱，J/(kg·K)；Δ

h

、Δ

h

分別為水-冰相變潛熱及冰-水蒸氣升華潛熱，J/kg。同時，還可以分別寫出水層與冰層的導熱關系如式(4)、(5):

式中，

λ

、

λ

分別為冰和水的導熱系數，W/(m·K)；

t

、

t

為升華溫度、底面溫度，K。需要說明的是，式(4)僅考慮了水層的穩態導熱，而忽略了入口較高溫度供水的影響。王玉瑩在相關仿真中做了類似的處理，其仿真與實驗結果較為一致，可以認為該處理方式的誤差在可接受范圍內。

此外根據式(6)所示，升華模式下水層、冰層的厚度之和與給水腔厚度相等:

式中，

d

、

d

、

d

分別為冰層厚度、水層厚度、給水腔厚度，m。根據式(3)，可推導出水升華器消耗工質水的質量流量如式(7):

式中，

ρ

為水的密度，kg/m；，

S

、

u

分別為供水入口的面積與供水水流速度。由于升華質量流量很小，水蒸氣分子平均自由程遠大于流動特征長度，具有較大的克努森數

Kn

，可以認為處在自由分子流的流動狀態。自由分子流狀態下的質量流量為式(8):

式中，

n

為多孔板單孔個數，

?

為多孔板孔徑，m；

d

為水蒸氣流動長度，m；

m

為單個水分子質量，kg；

k

為玻爾茲曼常數，

p

為真空環境壓力，Pa；飽和蒸氣壓

p

可根據式(9)給出:

此外，將式(8)帶入式(3)中，可得水升華器的升華熱流如式(10):

由于水不進入多孔板，此時水蒸氣的流動長度

d

即多孔板厚度

d

。聯立式(7)～(9)可以數值求解升華溫度

t

，并根據式(4)～(6)計算冰層厚度

d

、水層厚度

d

及底面溫度

t

。若熱載荷大于一定值，使升華表面蒸氣壓

p

大于水的三相點壓力

p

、升華溫度

t

超過三相點溫度

t

，水會進入多孔板轉換為周期模式。存在一個水升華器以升華模式工作的臨界熱載荷，標志升華模式向周期模式的轉換，基于式(7)～(8)可得式(11):

當水升華器的熱載荷小于式(11)時，認為處于升華模式；而當熱載荷大于該值時，則認為水升華器處于周期模式下。此外，升華模式下熱載荷還需要滿足式(12)的限制，避免給水腔被冰層完全填充可能造成的損壞。這個最低熱載荷的數值同樣可以通過數值方法計算得出。

式中，

d

｜代表熱載荷

Q

下根據導熱關系式及式(6)計算得出的冰層厚度。

綜上，對于尺寸確定的水升華器，其升華模式穩定工作的最高與最低熱載荷均是一個常數，可以通過式(11)～(12)確定。取升華面積0.04 m、多孔板厚度1 mm、給水腔厚度3 mm，可得圖5所示的結果，多孔板孔度、孔徑的增大會使熱載荷范圍變大、總體數值升高。

圖5 升華模式工作范圍Fig.5 Working scope of sublimation mode

2.3 熱載荷公式量綱分析

將式(11)中的小孔個數

n

用孔度

ε

展開并移項，兩邊再同除質量流量

m

，可得式(13):

將上式進行量綱分析可以得到式(14):

式中，將式(13)的左邊將臨界熱載荷

Q

替換為實際熱載荷

Q

，可定義斯坦頓數

St

、雅各布數

Ja

倒數差為式(15):

式中，

ρ

為水蒸氣密度。定義為一個新的無量綱數

Su

如式(17)，包含水升華器的孔度、孔徑、面積等所有設計參數:

結合密度、速度、面積比項，定義無量綱質量流量如式(18):

根據2.2節的內容，當熱載荷大于一定值時水升華器將轉入周期模式。由于此前的分析保留了模式轉換對應的三相點參數，因此熱載荷

Q

取臨界熱載荷

Q

時，可將式(11)轉化為如式(14)所示的無量綱形式，反映臨界熱載荷與水升華器設計參數之間的關系。對于孔度、孔徑及熱載荷分別在0.3～0.7、1～9μm、100～300 W間隨機取值，可得到如圖6所示的水升華器工作模式分布及模式轉換的無量綱特性曲線。對于水升華器不同設計參數組合的無量綱數

Su

，若熱載荷對應的

St

取值在圖中的特性曲線以上，可以認為水升華器處于周期模式下工作，反之則認為進入了升華模式。此外，由式(12)決定的升華模式下邊界曲線同樣可以在圖中給出。

圖6 模式轉換無量綱特性曲線Fig.6 Dimensionless curve of mode transition

3 仿真模型

為了進一步探究水升華器升華模式的具體特性，基于2.2～2.3節理論分析，進行水升華器升華模式的數值建模與仿真計算。

3.1 焓-多孔模型

本文采用焓-多孔模型求解固-液相變問題。該模型以溫度和焓作為共同變量，在全計算區域內建立統一形式的能量守恒方程。溫度和焓作為共同因變量的關系可以式(19)用來表示:

引入糊狀區(Mushy zone)的概念，將溫度處于固相線溫度與液相線溫度之間的區域視為糊狀區，該區域內網格的液相分數視為孔度，通過不斷更新所有網格單元內的孔度追蹤固液界面的位置。液體分數的表達為式(20):

對于水-冰相變問題，其固相線與液相線溫度相等，均為273.15 K。糊狀區的動量表達為式(21):

式中，

β

為液體分數，

ε

為防止分母為0的小量，

A

為糊狀區常數，

V

為速度矢量，

V

為固態離開原區域的漂移速度。基于式(20)～(21)，可得相變問題的控制方程組如式(22)～(24)所示。

連續性方程為:

由于水升華器工作在失重的空間環境中，可以忽略重力的影響，得到動量方程:

能量方程為:

式中，總焓

H

為顯焓

h

與潛熱焓Δ

H

之和，關系如式(25):

式中，

h

為參考焓值，J/kg；

T

為273.15 K的參考溫度；

L

為水-冰相變潛熱，J/kg。

3.2 網格劃分及Fluent求解模型

使用商業軟件ANSYS ICEM CFD及ANSYS Fluent 2020 R1版本作為網格劃分工具及求解器。由于供水流量很小，流動模型選擇層流Laminar。能量方程的求解調用Fluent中基于焓-多孔模型的Solidification/Melting模型，選擇基于壓力的隱式穩態單精度求解器進行求解，當整體熱量達到平衡后計算結束。

仿真分析相同外形的水升華器在不同孔度、孔徑時的升華模式特性，水升華器結構參數及相關常數如表1所示，其中3個工況已在圖7中標出。

表1 水升華器參數及相關常數Table 1 Parameters of sublimator and relative constants

對給水腔進行二維矩形結構網格的劃分，并且選取了數量為7150、8400、10500的網格進行仿真網格無關性驗證計算。3套網格的工況2計算結果如圖7所示，不同網格數量下的計算結果最大差異不超過0.15%。

圖7 網格無關性驗證Fig.7 Grid independence validation

綜合考慮計算的準確性及計算資源的占用，最終選取600×14共8400數量的網格，網格劃分及邊界條件設置情況如圖8所示。其中，底面熱載荷所在的加熱面設置為等熱流密度邊界，為根據圖6確定的熱載荷區間選取的合適的熱載荷。頂面為升華界面，定義為由溫度決定的熱流邊界，與孔度、孔徑及升華溫度等有關，如式(9)～(10)所示。計算域右側加入一個由式(7)確定的等質量流量入口，同時在頂面施加一個質量源項UDF模擬升華帶走的質量流。邊界條件結果如表2所示。

圖8 網格與邊界條件Fig.8 Mesh and boundary condition

表2 不同工況邊界條件取值Table 2 Boundary condition settings

4 結果與分析

由于Fluent中是基于式(26)計算液體分數的，松弛因子

γ

的存在會導致計算域中出現糊狀區。為避免不確定性，采用

t

=273.15 K等溫面表示水-冰界面如式(26):

三工況整體分布較為相似，取工況2的計算結果進行分析。整體及入口附近液體分數云圖如圖9所示，紅色代表液態水，綠色、藍色代表相變形成的糊狀區及固態冰。可見受到溫度較高的供水的影響，入口附近幾乎沒有形成冰層；而向后直至給水腔盡頭，入口效應逐漸消失，形成了厚度較為一致且水-冰界面與上下邊界大致平行的薄冰層。

圖9 液體分數云圖Fig.9 Contour of liquid fraction

工況2的流場矢量圖如圖10所示。受固液界面位置的影響，水流從入口進入后僅在下側的液態區域中流動，上側的固態冰層區域中不存在流動；隨后水流繼續向前流動遠離入口，速度逐漸降低并最終停止流動。此外，還能觀察到在冰面及底部壁面附近的流動邊界層效應。

圖10 速度矢量圖Fig.10 Vector of velocity

工況2的溫度分布云圖如11所示。可以發現，溫度分布與液體分數分布具有很強的關聯性。同樣將入口附近放大，可見入口的高溫水流前進的同時很快地受到了冷卻，等溫線逐漸傾斜；入口效應逐漸消失后，受到上下熱流邊界的影響，等溫線也逐漸平行于上下邊界。此外，還能觀察到右下角區域在底面加熱作用下達到了全域最高溫度，已超過了293.15 K的供水溫度。

圖11 溫度云圖Fig.11 Contour of temperature

統計各工況計算穩定后的底面溫度、升華溫度與其對應的冰層厚度的仿真結果與理論計算值對比，如表3所示。各工況的溫度結果相對誤差都在0.5%以內，底面溫度仿真結果整體比理論值高1 K左右，而升華溫度仿真結果比理論值稍低。造成這種誤差的原因是計算域頂部升華熱流、升華溫度耦合的邊界條件設置，其熱流是由升華界面平均溫度決定的，因此入口高溫水流對升華界面的邊界有著不可忽略的影響。

表3 溫度、界面位置結果對比Table 3 Comparison of temperature and thickness of ice

各工況冰層厚度的計算結果與理論值相對誤差的絕對值均小于1%。同樣受到入口效應的影響，計算域右側沒有形成冰層，同時溫度的誤差也會影響冰層厚度，因此冰層平均厚度整體較理論值稍低，對應的相對誤差偏向于負值。

5 結論

1)分析推導得到了水升華器工作模式轉換的特性曲線，該曲線能給出不同熱載荷及孔度、孔徑等結構參數下水升華器的工作模式；

2)對不同孔度、孔徑共3種工況進行仿真并與理論值比較，各溫度結果的相對誤差均在0.5%以內、冰層厚度的相對誤差絕對值均在1%以內，仿真模型通過了理論驗證；

3)入口附近較高溫度的供水會造成等溫線傾斜、無冰層形成等入口效應及冰層厚度仿真結果偏低等誤差。

本文僅研究了低熱載荷下的升華模式，可在后續研究中開展對熱載荷較大時的周期模式的理論分析及數值仿真工作。