微重力下相變儲能單元融化過程數值模擬

阮世庭， 張濟民， 曹建光， 王江， 徐濤

(上海衛星工程研究所, 上海 201109)

航天器在軌運行時，軌道外熱流周期性變化導致某些載荷的熱環境也呈現周期性變化[1]。相變熱控技術特別適用于具有周期性工作的設備和部件[2]。

當相變材料與發熱部件的接觸界面溫度高于相變材料的熔點時，相變材料融化并吸收發熱部件釋放的熱量，使界面溫度保持在相變點附近；當接觸界面溫度低于相變材料的熔點時，相變材料凝固并釋放潛熱，維持界面溫度基本不變[3]。Leimkuehler等[4-5]選擇純水作為相變材料，針對在低月軌道的航天器和月球車所處的工作環境分別設計了相變溫控單元以及相變熱沉。Ye等[6-7]研究了金屬肋片強化傳熱的相變熱沉內部傳熱與流動過程。Ismail等[8]和Abduljalil等[9]研究肋片對環狀相變材料傳熱性能的影響。Assis等[10-11]研究了球狀相變材料融化與凝固過程，通過量綱分析研究液相質量分數和熱流密度隨時間的關系。

但是，大部分研究人員沒有對微重力環境中相變過程進行數值模擬或實驗研究。基于此，本文通過數值模擬方法探究微重力環境下相變過程中的傳熱與流動特性，并對比重力作用下相變過程，為航天器相變熱控技術提供參考依據。

1 物理模型與數值方法

1.1 物理模型

數值模擬的物理模型如圖1所示。相變儲能單元由基座、肋片、相變材料以及具有一定真空度的空氣組成。氣體和相變材料的體積比為1∶4。空氣的壓力為相變材料在此溫度下的飽和蒸汽壓(約為100 Pa)，預留空氣以便觀察相變材料融化后體積膨脹過程。選用鋁作為基座與肋片材料，增加肋片以強化相變材料與加熱面之間的傳熱，相變材料選用十八烷。計算單元左右兩側設定為對稱面，頂部設為絕熱壁面，底部為溫度加熱面。

圖1 計算單元Fig.1 Computational unit

計算單元尺寸如圖1(b)所示。肋片與相變材料的寬度比為1∶1，肋片高度hfin為22 mm，相變材料腔體高度hPCM為20 mm，氣體高度hair為4 mm，氣體寬度wair為3 mm。十八烷沒有嚴格的相變溫度，定義其相變溫度區間為27～29℃。相對于一個固定的相變點，相變溫度區間更加符合真實的物理現象。十八烷的密度隨溫度變化，當溫度T低于27℃時，密度ρ=814 kg/m3，當溫度高于27℃時，密度表達式為

(1)

式中：ρl為溫度處于Tl時相變材料的密度；β為相變材料的熱膨脹系數。

當27℃29℃時，取β=0.001[12]。鋁、空氣以及相變材料的其他物性參數如表1所示。

計算單元初始溫度為25℃，相變材料處于過冷狀態。加熱面處給定的恒定溫度邊界條件T=48℃，邊界溫度與相變材料平均相變溫度溫差ΔT=20℃。

表1 各物質物性參數

1.2 控制方程

采用焓-多孔介質模型求解相變材料的相變過程，多孔介質區域的每個單元內設置相同的流動阻力。對于全凝固區域和全融化區域，多孔性分別為0和1。Bertrand等[13]將多種數值模擬方法運用于相變過程的數值模擬中，相變過程中考慮自然對流對融化的影響，并且覆蓋2個系列的Prandtl準則數，這2個系列的準則數分別對應金屬和有機材料。結果表明，焓-多孔介質模型能夠很好地應用在具有移動界面的固-液相變問題中。

控制方程如下：

連續方程

(2)

式中:αn為計算單元中第n種流體的體積分數;t為時間。

動量方程

(3a)

(3b)

式(3a)為受重力作用時動量方程；式(3b)為不受重力作用時動量方程；v為速度，P為壓力，μ為動力黏度，g為重力加速度，S為動量源項。

能量方程

(4)

(5)

動量方程中的源項S=-A(γ)v，其中A(γ)為Brent等[14]定義的“多孔函數”。定義A(γ)使得動量方程能夠模擬多孔介質中流動的Carman-Kozeny(卡爾曼-科澤尼)方程。定義如下：

(6)

式中：ε為恒等于0.001的計算小量，用來消除分母為0時產生的振蕩；C為反映融化前沿形態的模糊區常數，一般取104～107，本文取C=105[15]。

采用Fluent14.5進行求解，使用雙精度求解器和SIMPLE算法進行壓力-速度的耦合。選取3種不同網格數量進行網格獨立性驗證，分別是2 700、4 500和6 200。圖2(a)顯示3種不同的網格數量下，相變材料液相質量分數隨時間變化情況，從圖2(a)可以看出，不同網格數對相變過程的影響差別很小，本文選取網格數為2 700進行數值模擬，具體網格劃分見圖2(b)。非穩態時間步長等于比特征長度除以比特征時間，本文取時間步長Δt=0.005 s。每個時間步長內，確保連續性方程和動量方程殘差小于10-6，能量方程殘差小于10-10。

圖2 網格無關性驗證和網格劃分Fig.2 Grid independence verification and grid partition

1.3 實驗驗證

因為實驗條件限制，無法搭建微重力條件下的實驗臺。搭建地面實驗臺，觀察受重力影響時相變儲能單元相變過程。實驗裝置如圖3所示。實驗開始前，先對實驗件進行抽真空處理，抽完真空后，往實驗件中填充體積分數為80%的相變材料。待實驗件充裝完畢后，將其放置在水箱中。接著，打開泵和閥門，從恒溫水槽內通入48℃熱水，從而達到恒定溫度邊界條件。通過相機拍攝實驗結果。

恒溫水槽采用DC-1030低溫恒溫槽，溫度范圍為-10～90℃，相機為佳能EOS70D單反相機。

選取相變過程中8個時間節點作對比，分別是0、Δtmelt/7、2Δtmelt/7、3Δtmelt/7、4Δtmelt/7、5Δtmelt/7、6Δtmelt/7、Δtmelt，Δtmelt為相變材料完全融化的時間。實驗與數值模擬如圖4所示，其中白色區域是固態相變材料，黑色部分是液態相變材料，由于部分固態相變材料存在于液態相變材料后面，影響觀察結果，導致部分固-液相變界面比較模糊，不過仍可以從圖4看出，數值模擬結果與實驗結果相吻合，數值模擬結果具有參考價值。

圖3 實驗裝置Fig.3 Experimental setup

圖4 相變材料融化過程數值模擬與實驗結果對比Fig.4 Comparison of numerical simulation and experimental results for melting process of phase change material

2 結果分析與討論

為了更加深入了解微重力環境下，相變材料融化過程中各時刻狀態變化，將同一計算單元進行受重力影響和受微重力影響2種情況的數值模擬，對比分析2種情況下相變材料融化過程的異同。給定恒定的溫度邊界條件T=48℃，初始溫度T0=25℃。

2.1 融化速率與熱流密度

圖5和圖6是2種情況下，相變材料液相質量分數和邊界熱流密度對比情況，邊界熱流密度為加熱壁面傳遞給相變儲能單元的熱流密度。

從圖5可以看出，當相變儲能單元受重力作用時，其完全融化耗時70 s；當相變儲能單元受微重力作用時，其完全融化耗時90 s。相對于受重力作用，當相變儲能單元受微重力作用時，相變材料融化速率明顯下降。融化速率受邊界熱流密度影響，當相變儲能單元受微重力作用時，其邊界熱流密度明顯小于相變儲能單元受重力作用時的熱流密度。

當時間t<15 s時，2種情況下液相質量分數和邊界熱流密度沒有明顯的區別；當時間t>15 s且t<20 s時，2種情況下液相質量分數沒有明顯的差異，但是受重力作用的相變儲能單元邊界熱流密度明顯大于受微重力作用的相變儲能單元邊界熱流密度，說明此階段中，受重力影響的相變儲能單元內相變材料吸收更多的熱量，并將這部分熱量轉化為相變材料的顯熱。

圖5 液相質量分數隨時間變化Fig.5 Liquid fraction versus time

圖6 邊界熱流密度隨時間變化Fig.6 Heat flux density versus time

2.2 速度分布

圖7為2種情況下相變儲能單元內相變材料的速度分布，取液相質量分數φ=0.2、0.4、0.6、0.8和1.0等5個時刻表示整個融化過程。

從圖7可以看出，有無重力作用對相變儲能單元內部速度分布起決定性作用。當相變儲能單元受重力作用時，肋片處液相相變材料內部產生自然對流，并且對流換熱是主要的換熱形式[6]。當相變儲能單元受微重力影響時，液態相變材料速度由相變材料膨脹產生，其數量級大概為1.0×10-5m/s，無對流換熱，熱量主要通過熱傳導傳遞。

2.3 溫度分布

圖8為2種情況下相變儲能單元內相變材料的溫度分布圖，同樣取液相質量分數φ=0.2、0.4、0.6、0.8和1.0等5個時刻表示整個融化過程。

相變過程初始階段，相變儲能單元內部溫度分布并沒有明顯的區別。當融化進行到一定階段，二者開始出現差異。當相變儲能單元受重力影響時，內部自然對流使得溫度較低、密度較大的相變材料進入底部，底部出現局部低溫區域。當相變儲能單元受微重力影響時，內部溫度分布無自然對流影響，局部低溫區域出現在相變材料頂端，由于相變儲能單元頂部預留部分氣體，此部分氣體比熱容很小，溫度上升快，所以相變材料低溫區域出現在相變儲能單元中上部區域。

2.4 固-液兩相分布

2種情況下相變儲能單元固-液兩相分布如圖9所示，當相變儲能單元受重力影響時，融化的相變材料膨脹從頂端溢出，受重力作用后，液相相變材料覆蓋于頂部；當相變儲能單元受微重力影響時，融化的相變材料從頂端膨脹溢出，無重力作用時，向空間擴散。受重力作用時，未融化的相變材料下沉；無重力作用時，未融化相變材料懸浮于已融化相變材料中間，不出現下沉。

圖7 相變材料融化過程中速度分布 Fig.7 Velocity distribution for melting process of phase change material

圖8 相變材料融化過程中溫度分布Fig.8 Temperature distribution for melting process of phase change material

圖9 相變材料融化過程中固-液兩相分布Fig.9 Solid-liquid distribution for melting process of phase change material

3 結 論

通過數值模擬方法對相變過程在重力和微重力作用的2種情況進行研究，并對受重力作用的相變過程數值模擬結果進行實驗驗證。數值模擬結果表明：

1) 相對于受重力作用，當相變儲能單元無重力作用時，相變材料融化速率明顯下降。融化速率受邊界熱流密度影響，當相變儲能單元受微重力作用時，其邊界熱流密度明顯小于相變儲能單元受重力作用時的熱流密度。

2) 當相變儲能單元受重力作用時，肋片附近液相相變材料內部產生自然對流，并且對流換熱是主要的換熱形式。當相變儲能單元受微重力影響時，液態相變材料速度是由于相變材料膨脹產生，無對流換熱，熱量主要通過熱傳導傳遞。

3) 當相變儲能單元受重力影響時，內部自然對流使得溫度較低、密度較大的相變材料進入底部，底部出現局部低溫區域；當相變儲能單元受微重力影響時，內部溫度分布無自然對流影響，局部低溫區域在相變儲能單元中上部。

4) 當相變儲能單元受重力影響時，融化的相變材料膨脹從頂端溢出，受重力作用后，液相相變材料覆蓋于頂部；當相變儲能單元受微重力影響時，融化的相變材料從頂端膨脹溢出，微重力作用時，向空間擴散。受重力作用時，未融化的相變材料下沉；無重力作用時，未融化相變材料懸浮于已融化相變材料中間，不下沉。