基于有限體積法堿金屬高溫熱管冷態啟動流動換熱數值研究

趙昊城,張澤秦,王成龍,秋穗正,田文喜,蘇光輝

(西安交通大學 核科學與技術學院,陜西 西安 710049)

熱管作為一種非能動熱傳輸裝置,主要通過工質相變傳熱,可以自發運行,輸送工質不需要額外的能量輸入,是一項廣泛應用的創新型輸熱元件[1],熱管具有極高的熱導率及良好的等溫性,能夠在傳熱的同時保持加熱段和冷凝段溫度幾乎均勻[2]。為在高溫環境中長期穩定運行,高溫熱管常選用堿金屬作為傳熱工質[3]。由于高溫熱管工作溫度較高,且工質多為化學性質活潑的堿金屬,實驗研究工作難度較大[4],為此國內外學者提出不同的啟動模型對堿金屬熱管冷態啟動進行模擬。Cao等[5]提出采用自擴散模型描述蒸氣區自由分子流態流動特性;Tournier等[6]開發了適用于堿金屬熱管的冷態啟動分析模型HPTAM,是目前認可度最高的啟動模型,HPTAM完整考慮了吸液芯結構帶來的多孔介質阻力項和慣性阻力項,建立了吸液芯的完整動量方程,對于蒸氣區采用塵氣模型來模擬自由分子流、過渡流及連續流動[7]。Hansel等[8]基于MOOSE有限元框架,采用一維建模兩相可壓縮流動模擬,基于兩流體七方程推導熱管控制方程,對吸液芯部分液相流動采用達西方程描述。國內方面,Zhang等[9-10]開發了基于有限元方法的熱管冷態啟動三階段模型。目前,國內的熱管數值研究大多忽略了吸液芯部分流動,國際熱管數值模擬中吸液芯部分流動相關研究資料較少并缺乏驗證。

本工作針對高溫堿金屬熱管啟動階段工質傳熱和流動進行研究,使用有限體積法(FVM)建立熱管瞬態分析模型,開發熱管瞬態模擬程序,對單根水平鈉熱管冷態啟動過程進行數值模擬,分析熱管溫度和壓力分布以及蒸氣各項參數變化,研究環境溫度和絕熱段長度對熱管啟動特性的影響。

1 數學物理模型

考慮熱管為軸對稱的圓柱結構,選取合適的特征長度,進行軸向和徑向的二維板型計算域簡化,控制區域如圖1所示。圖1中:Qin為蒸發段外界輸入功率;Qout為冷凝段與外界換熱功率;l為熱管整體長度;le、la和lc分別為蒸發段、絕熱段以及冷凝段長度。對該模型主要做出以下假設:1) 忽略熱管中工質受到的體積力;2) 不考慮蒸氣區中不凝氣體的影響;3) 熱物性僅為溫度的函數;4) 蒸氣為一維層流流動。

圖1 熱管模型示意圖Fig.1 Schematic diagram of heat pipe model

1.1 管壁傳熱模型

外壁區傳熱過程使用笛卡爾坐標系中的熱傳導方程[11]描述,控制方程為:

(1)

式中:t為時間;x、y分別為軸向和徑向方向;(ρwacp)wa為管壁體積熱容;Twa為管壁溫度;kwa為管壁導熱系數。

1.2 吸液芯多孔介質流動換熱模型

吸液芯內工質流動的連續方程為:

(2)

式中:ε為吸液芯多孔介質孔隙率;ρwf為液相工質密度;u和v分別為軸向和徑向溫度。

采用達西方程[12]模擬吸液芯內部工質流動:

(3)

式中:μv為液態工質動力黏度系數;K為多孔介質滲透率;p為吸液芯內相對壓力分布;u為液態工質流速。

吸液芯內能量方程為:

(4)

式中:Twi為吸液芯內熱平衡溫度;(ρcp)eff為吸液芯等效體積熱容;keff為等效導熱系數。

等效物性參數由Chi模型[13]計算,具體計算方法為:

(ρcp)eff=ε(ρcp)wf+(1-ε)(ρcp)wi

(5)

(6)

式中:(ρcp)wf為工質體積熱容;(ρcp)wi為吸液芯材料體積熱容;kwf為工質導熱系數;kwi為吸液芯材料導熱系數。

吸液芯內部方程中速度均為流體通過整個多孔介質區域的速度,即表觀速度[14]。

1.3 蒸氣區域模型

根據傳質學研究[15],克努森數Kn常用于表征蒸氣自由分子態與連續態的轉變閾值。熱管處于冷態時,蒸氣區壓力很小,根據Kn定義,此時Kn應是一個很大的值,表明氣體分子間相互作用影響小于氣體分子與孔隙壁面碰撞產生的影響,隨著熱管操作溫度的上升,Kn逐漸減小。本文定義Kn=0.01為轉變閾值,即當Kn<0.01時認為蒸氣區進入連續流動狀態。

(7)

式中:λ為蒸氣分子平均自由程;D為熱管蒸氣區直徑。

蒸氣自由分子態與連續態的轉變閾值可由Kn[16]得到:

(8)

式中:M為相對分子質量;D為蒸氣空間寬度;Ru為理想氣體常數;μ為蒸氣動力黏度系數;ρ為蒸氣密度。

自由分子態時計算忽略蒸氣區和氣液界面,隨溫度升高,蒸氣區中出現連續流動部分,并隨熱管啟動時間增加逐漸擴展,連續流動部分與自由分子部分分隔界面稱為連續鋒面,連續流態下認為蒸氣處于準穩態運行狀態,當熱管蒸氣區到達該狀態時,采用文獻[16]中給出的一維可壓縮穩態層流流動控制方程進行模擬。這是因為在該階段雖然熱管溫度仍在繼續上升,但蒸氣區內蒸氣流動狀態變化十分緩慢,可使用穩態方程進行處理。在計算控制方程中所需的摩擦因子、動量因子和能量因子時,由于熱管蒸氣區的通道寬度大于高度,視為二維模型進行處理,且在一個面上存在與該面法向方向相同的注入流速,與半多孔介質通道特征相同,因此采用文獻[17]中半多孔介質通道相似解進行處理。

1.4 邊界條件

熱管各邊界所使用的邊界條件不同,各部分邊界具體劃分如圖1所示。

1) 熱管沿軸向劃分為蒸發段、絕熱段和冷凝段,這3部分邊界條件分別設定恒定功率輸入邊界、絕熱邊界和對流換熱邊界:

(9)

式中:hc為外壁面與環境等效換熱系數;Tamb為環境溫度;qin為蒸發段外部熱流密度。

2) 管壁與吸液芯之間主要通過導熱的方式進行傳熱,忽略工質流動造成的對流換熱。

3) 熱管軸向兩端均設為絕熱邊界條件,吸液芯部分添加無滑移邊界條件:

v=0x=0和x=l

(10)

4) 根據蒸氣的不同流態,吸液芯與蒸氣區界面采用不同邊界條件。當蒸氣處于自由分子態時,認為氣液界面絕熱;處于過渡流態時,根據氣體動理論[18]得到界面質量流量。忽略蒸氣區能量和質量的累積,并綜合考慮聲速極限的限制[19]。

未達到聲速極限

(11)

(12)

吸液芯區和蒸氣區通過氣液界面能量方程進行耦合計算:

(13)

式中,kwv為吸液芯中氣液界面導熱系數,定義為吸液芯區域在氣液界面上的導熱系數,在計算過程中根據界面溫度和物性關系式進行計算。

氣液界面邊界上流速與界面質量通量滿足:

(14)

2 程序開發及驗證

采用FVM,基于C語言環境開發了高溫熱管啟動瞬態分析程序HPWF。對熱管管壁及吸液芯部分分別進行二維建模,使用矩形網格進行劃分,如圖2所示。

圖2 程序網格劃分Fig.2 Program grid division

2.1 程序數值方法

在劃分的控制容積中心[11]設置節點,為避免求解出“棋盤式”壓力分布,采用交錯網格[11];分別選取合適的松弛因子,使用超松弛迭代法求解管壁和吸液芯部分壓力方程和能量方程;使用牛頓迭代法求解氣液界面能量守恒非線性方程;蒸氣區參數的常微分方程組通過4階龍格-庫塔法進行求解。

2.2 程序驗證

首先進行程序網格無關性驗證。研究對象設置為長1 m的熱管,對于長徑比較大的熱管,軸向網格劃分對計算結果影響較大,因此在管壁徑向上劃分2層網格,吸液芯徑向上劃分4層網格,在軸向上分別劃分25、40、50和67層網格進行計算,選取蒸氣區最高溫度和計算時間作為對比參數。程序網格無關性驗證結果如圖3所示。當網格數量從300增加到402時,蒸氣區最高溫度變化很小,而計算時間增加25.7%,計算成本大大增加,因此選擇軸向劃分50層網格、網格數量300作為計算工況。

a——蒸氣區最高溫度和計算時間分布;b——氣液界面相變質量通量分布;c——熱管外壁面溫度分布圖3 程序網格無關性驗證Fig.3 Program grid independence verification

為驗證HPWF程序對熱管啟動瞬態分析的正確性和準確性,選用Lee等[20]測試的鈉熱管實驗數據進行驗證。實驗采用長為1 m的鈉熱管,使用可變電壓器控制加熱功率,冷凝段外壁使用空氣自然對流冷卻,將6個K型熱電偶裝入小型屏蔽管,插入熱管中心以測量蒸氣溫度。實驗采用熱管的具體參數列于表1。

表1 鈉熱管結構參數Table 1 Structure parameter of sodium heat pipe

熱管啟動過程中,蒸氣溫度變化與穩態蒸氣溫度分布的模擬結果與實驗數據對比如圖4所示。由圖4可見,模擬結果與實驗數據符合較好,最大相對誤差約為9.8%。

a——啟動瞬態;b——穩態運行圖4 模擬結果與實驗數據的對比Fig.4 Comparison of simulation result and experimental data

3 計算結果

對單根水平鈉熱管冷態啟動瞬態過程進行計算分析,熱管蒸發段長為0.4 m,絕熱段長為0.2 m,冷凝段長為0.4 m,熱管外徑為30 mm,管壁厚度為4 mm,吸液芯厚度為2 mm,材料為316L不銹鋼,蒸發段外壁面邊界熱流為定常輸入,輸入功率設定恒定的1 000 W,換算熱流密度為53 051.65 W/m,冷凝段外邊界換熱系數設為100 W/(m2·K)。

3.1 啟動特性分析

圖5示出熱管冷態啟動壁面溫度軸向分布。啟動總用時3 000 s,壁面軸向溫差最終穩定在22 K。由于各段給定邊界條件均勻,啟動階段任一時刻蒸發段,以及初始階段和準穩態穩態運行階段絕熱段、冷凝段溫度分布均勻。

圖5 壁面溫度軸向分布Fig.5 Axial distribution of wall temperature

啟動開始時蒸發段溫度快速升高。t=200 s時,連續鋒面進入絕熱段,軸向溫度分布曲線上出現了很大的溫度梯度;t=300 s時,絕熱段整體溫度到達776 K溫度平臺,連續鋒面進入冷凝段,熱管整體的溫度曲線呈三階梯式分布,蒸發段與冷凝段連續流區域的溫差仍較大,約為60 K;t=600 s后,溫度曲線分布形狀不再發生改變;t=800 s時,軸向最大溫差維持在穩定值。連續流動鋒面在冷凝段移動過程中,已到達連續流態的冷凝段部分平臺溫度略有下降,下降幅度約為15 K,這是由于冷凝段升溫過程中,軸向傳熱量被限制在較低的聲速極限并保持不變,而程序假設換熱系數為較大的定值,當冷凝段升溫區域增加時,外界環境換熱量增大導致冷凝段溫度降低。

圖6示出熱管啟動過程中蒸氣區理想熱流和聲速極限隨時間的變化。在t=150～700 s階段,實際傳熱受到聲速極限限制。t=570 s之前,熱管的理想熱流量上升速度很快;當t=570 s時,熱管的理想軸向熱流量達到最大值14 090.3 W,而此時聲速極限將傳熱量限制在600.1 W;t=570～720 s時,理想熱流量迅速下降,此階段蒸氣區整體達到連續流態,冷凝段氣液界面溫度快速升高,直到熱管整體溫差保持恒定。

圖6 理想熱流與聲速極限的變化Fig.6 Variation of ideal heat flux and limit of sound velocity

圖7示出熔化前沿與連續鋒面相對位置隨時間的變化,熔化前沿是指吸液芯中固液交界面所在的位置。t=50 s時,熔化前沿出現在蒸發段與絕熱段交界面上;t=150 s時蒸氣區中出現連續流動區域,連續鋒面開始沿軸向推進。熱管軸向熱流受到聲速極限限制,熔化前沿和連續鋒面沿軸向以恒定速度推進,第2階段末期(t=400～650 s),推進速度略有降低,最終在t=600 s和t=650 s時熔化前沿和連續鋒面先后到達冷凝段末端,熱管進入準穩態運行階段。

圖7 熔化前沿和連續鋒面位置的變化Fig.7 Position variation of melting front and continuous front

設置蒸發段端部壓力為0作為參考壓力,圖8示出吸液芯內部壓力相對值分布。由圖8可見:熱管吸液芯徑向壓力變化不到1 Pa,軸向壓力梯度較大,穩態壓降可達到47 Pa;當連續鋒面到達冷凝段末段后,吸液芯整體壓降達到最大值,啟動繼續進行,整體壓降約7.73 Pa,這是因為工質的動力黏度系數隨溫度升高而降低,相同速度條件下壓力梯度略有降低。

圖8 吸液芯壓力的變化Fig.8 Variation of pressure of wick

圖9示出吸液芯穩態流速分布。由圖9可見,徑向流速很小,最大流速僅為0.01 mm/s,冷凝段末端徑向流速很大,達到6 mm/s。由達西方程可知,在吸液芯滲透率很小的情況下,一個很小的壓力波動會引起流速巨大變化,而吸液芯內部流速變化對溫度分布的影響很小,因此實際上壓力和溫度分布都無明顯變化。絕熱段處流體的軸向流速絕對值最大,速度剖面呈梯形,而最大軸向速度也僅約3 mm/s。圖9中負流速表示液相流動方向與軸向方向相反,軸向速度梯度為0.007 5 s-1。

圖10示出蒸氣區完全進入連續流態后不同時刻熱管蒸氣各項參數分布。在相變質量的累積引起的蒸氣流加速壓降和摩擦壓降共同作用下,蒸發段中蒸氣壓力急劇下降,溫度下降幅度也達到最大,約為3.8 K;絕熱段中由于摩擦壓降作用,壓力近似線性下降,且溫度繼續降低,由于鈉蒸氣膨脹引起密度減小,流速和馬赫數略有增加;冷凝段中蒸氣相變引起的質量損失導致速度降低,而摩擦壓降的影響始終存在,因此壓力無法完全恢復。

a——啟動1 112.40 s;b——啟動1 785.60 s;c——啟動2 418.20 s;d——啟動4 993.20 s圖10 熱管蒸氣參數的變化Fig.10 Variation of steam parameter of heat pipe

3.2 環境溫度敏感性分析

表2列出環境溫度對熱管啟動運行特性的影響。由表2可見,環境溫度降低10 K,穩態工作溫度也下降約10 K,完全啟動時間縮短約50 s。

表2 環境溫度對運行特性的影響Table 2 Influence of ambient temperature on operating characteristics

不同初始溫度穩態蒸氣壓力和流速的分布如圖11所示。由圖11可見,鈉蒸氣壓力隨環境溫度增加而增大,蒸氣流速隨環境溫度增加而降低,這是因為鈉蒸氣密度隨溫度升高而迅速升高,界面相變量保持恒定,導致蒸發段流速增大幅度下降。

圖11 穩態蒸氣參數分布Fig.11 Distribution of steady state steam parameter

3.3 絕熱段長度敏感性分析

表3列出不同絕熱段長度下熱管啟動運行特性。隨絕熱段長度增加,熱管啟動時間延長,氣相壓降和液相壓降增加,工作平均溫度相等。

表3 絕熱段長度對運行特性的影響Table 3 Influence of insulation section length on operating characteristics

不同絕熱段長度下穩態液相壓力和軸向流速分布如圖12所示,圖12中負流速表示液相流動方向與軸向方向相反。由圖12a可見,蒸發段和冷凝段流速變化相同,而絕熱段液相流速不變。這是因為液相不可壓縮,密度保持不變,因此界面質量通量相同的情況下,流速不發生變化,因此絕熱段長度并不影響液相最大流速以及蒸發段和冷凝段的速度變化趨勢。由圖12b可見,3種熱管液相壓力分布形狀一致,且蒸發段和絕熱段壓力變化形狀相同,而絕熱段壓降隨長度增加對應線性增加。

圖12 吸液芯穩態參數分布Fig.12 Parameter distribution of steady state of wick

4 結論

本文提出的模型采用Kn劃分熱管啟動不同階段,并采用達西方程描述吸液芯中液態工質的流動,對于蒸氣區采用準穩態一維可壓縮層流假設建立控制方程。本文建立了熱管啟動模型并開發了熱管冷態啟動瞬態分析的HPWF程序,使用公開文獻實驗數據驗證程序,最大相對誤差約為9.8%。本文模擬研究了高溫鈉熱管啟動過程中吸液芯壓力變化,并分析了啟動階段蒸氣各項參數的分布及變化關系,得到的主要結論如下。

1) 熱管冷態啟動過程中,溫度分布整體呈現為三階梯式,外壁面軸向溫差先增大后減小,最大溫差達到了487 K,到達穩態后工作溫度約為841.777 K,壁面溫差約為22.51 K。

2) 吸液芯內部壓力相對值隨工質熔化部分的擴展而增大,最大軸向壓降達到50 Pa,當工質完全熔化后,軸向壓降隨溫度升高略有降低,最終穩態工況下軸向壓降穩定在47 Pa。

3) 在蒸氣區部分處于連續流態時,熱管受到聲速極限限制,軸向傳熱量很小,當整體到達連續流態時,蒸氣溫度快速升高,氣液界面溫差迅速減小,最終突破聲速極限。

4) 環境溫度升高使得熱管穩態工作溫度升高,啟動時間延長,蒸氣區整體絕對壓力升高,最大流速增加,但對軸向溫差無明顯影響;絕熱段長度增加延長了熱管啟動時間,使得吸液芯和蒸氣區內壓降增大,而對穩態工作溫度和外壁面軸向溫差影響較小。