王嘯遠，朱躍釗，陳海軍，王銀峰，范紅途

?

熱虹吸管相變傳熱行為CFD模擬

王嘯遠，朱躍釗，陳海軍，王銀峰，范紅途

(南京工業(yè)大學機械與動力工程學院，江蘇南京，211816)

基于VOF模型和改進的相變模型對熱虹吸管穩(wěn)態(tài)相變傳熱過程進行CFD模擬，根據(jù)溫度分布實驗數(shù)據(jù)和可視化實驗進行驗證，并分析充液率對熱管蒸發(fā)?冷凝傳熱行為的影響。研究結(jié)果表明：CFD模型能夠較準確地描述熱管工質(zhì)的蒸發(fā)和冷凝過程，模擬得到熱管各段均溫的最大偏差為10.2 K，相對誤差2.9%，采用改進的相變模型，溫度分布模擬值與實驗值一致性更好；模擬得到蒸發(fā)段工質(zhì)流型轉(zhuǎn)變過程與可視化實驗定性相符；熱管冷凝段以膜狀冷凝傳熱為主，充液率的改變對其影響不大，充液率主要影響蒸發(fā)段的傳熱行為，進而影響熱管的傳熱性能。該工作有助于對熱管運行機理的理解，并為其CFD模擬提供指導(dǎo)。

熱虹吸管；相變傳熱；CFD；可視化

兩相閉式熱虹吸管(簡稱熱虹吸管)，是依靠重力回流冷凝液的高效傳熱元件，具有超高的導(dǎo)熱性以及優(yōu)良的等溫性，更兼結(jié)構(gòu)簡單和工作可靠的優(yōu)點[1]，在太陽能中高溫熱利用[2]、聚光光伏發(fā)電[3]、余熱回收[4?5]、以及電子元器件[6]的高效散熱等領(lǐng)域具有潛在的價值。熱管通過內(nèi)部工質(zhì)的相變循環(huán)過程進行傳熱，傳熱機理較為復(fù)雜[1, 7?8]。對于熱虹吸管的傳熱行為研究，傳統(tǒng)的數(shù)值模擬方法往往忽略了工質(zhì)的沸騰和流動，因此，與實際過程偏差較大。特別是對于壓力較高的工質(zhì)(如氨)或者高溫熱管工質(zhì)(如鈉)，這類熱管的啟動和相變傳熱行為研究缺乏有效的可視化實驗研究手段，更難直觀分析相變傳熱機理。CFD模擬技術(shù)可直觀地展現(xiàn)工質(zhì)的相變和流動狀態(tài)，有助于熱管運行機理和傳熱行為的深入理解。SCHEPPER等[9]提出了穩(wěn)態(tài)沸騰過程的CFD相變模型，隨后，學者們分別基于該模型利用Fluent軟件對以水為工質(zhì)熱管的傳熱性能[10]、管壁溫度分布[11]以及納米流體的強化傳熱作用[12]進行了模擬研究，但總的來說，熱管相變傳熱過程的CFD模擬研究不多。此外，現(xiàn)有熱管CFD模擬研究[10?12]采用的相變模型是針對穩(wěn)態(tài)沸騰過程提出的，較難適應(yīng)熱管內(nèi)蒸發(fā)?冷凝共存的相變傳熱特征，模擬偏差較大，同時，模擬性能缺乏實測可視化現(xiàn)象的驗證。本文作者擬以ANSYS Fluent 14.5軟件為平臺，對 Schepper相變模型進行改進，建立熱管CFD模型，模擬熱管的相變傳熱行為，結(jié)合溫度分布實驗數(shù)據(jù)和可視化實驗現(xiàn)象對模擬結(jié)果進行較全面的驗證，并根據(jù)模擬結(jié)果深入分析不同充液率下熱管的相變傳熱行為。

1 模型建立及求解

1.1 幾何模型

熱管工質(zhì)以軸向和徑向流動為主，而不是旋轉(zhuǎn)流動，故計算域可簡化為二維幾何模型。采用Gambit 2.4建立耦合固體管壁與汽液流場的銅?水熱管幾何模型，尺寸參照文獻[13]，總長980 mm，內(nèi)徑為25 mm，壁厚3 mm，蒸發(fā)段、絕熱段和冷凝段長度分別為295，305和380 mm。

1.2 CFD模型

基于ANSYS Fluent 14.5建立熱管的CFD模型，假設(shè)條件如下：1)相變過程壓力恒定，飽和溫度sat為定值；2) 汽液相及管壁物性參數(shù)均為常數(shù)；3) 忽略汽液界面剪切力；4) 汽液流動均為層流。

建立的熱管CFD模型采用VOF多相流模型[14]處理工質(zhì)的流動與傳熱過程，通過CSF(continuum surface force)模型[15]引入表面張力的影響，并對Schepper相變模型進行了改進。

1.2.1 VOF模型

常用多相流模擬方法有歐拉?拉格朗日法和歐拉?歐拉法2種[11]。歐拉?拉格朗日法將主相處理為連續(xù)相，其余相作為離散顆粒追蹤計算，僅適用于顆粒流問題的模擬。對于包含大汽泡流和分層流的熱管模擬，應(yīng)采用將各相均視為連續(xù)相的歐拉?歐拉法。VOF模型是歐拉?歐拉法的一種，其特有的相體積率幾何重構(gòu)算法能夠精確追蹤汽液界面的遷移，實現(xiàn)相變過程熱質(zhì)傳遞現(xiàn)象的模擬[10]。

VOF模型主要對連續(xù)性方程、動量方程和能量方程三大控制方程進行求解。由于熱管流場區(qū)域僅有汽液兩相，同一計算單元內(nèi)兩相相體積率之和為1，故連續(xù)性方程僅對液相進行求解，而動量方程與能量方程均針對汽液混合相而言，控制方程及混合相各物性參數(shù)的計算方法見式(1)~(7)。

連續(xù)性方程：

動量方程：

(2)

能量方程：

混合相物性參數(shù)計算方法：

(4)

(6)

(7)

1.2.2 相變模型

熱管的高效傳熱特性源于腔內(nèi)工質(zhì)的相變傳熱循環(huán)，因此，相變的模擬至關(guān)重要。工質(zhì)相變過程可簡化描述為質(zhì)量和能量的轉(zhuǎn)移，求解控制方程時實時修改質(zhì)量和能量方程的源項即可實現(xiàn)其模擬。源項的修改通過UDF(user-defined function)功能實施，即通過C語言將源項表達式編譯導(dǎo)入ANSYS Fluent 14.5求解器中，使其與控制方程相關(guān)聯(lián)。

SCHEPPER等[9]提出的穩(wěn)態(tài)相變模型如下式 所示。

蒸發(fā)過程(≥sat)：

(9)

(10)

冷凝過程(＜Tsat)：

(12)

(13)

式中：e和c分別為蒸發(fā)和冷凝過程的時間松弛因子；sat為飽和溫度；Δ為汽化潛熱。從式(8)~(13)可知：當高于sat時為蒸發(fā)過程，計算單元內(nèi)液相質(zhì)量減少，汽相增加，顯熱相應(yīng)轉(zhuǎn)化為潛熱；當?shù)陀趕at時為冷凝過程，熱質(zhì)傳遞方向恰好與蒸發(fā)過程相反。

e和c分別控制工質(zhì)的蒸發(fā)和冷凝速率，SCHEPPER等[9]在穩(wěn)態(tài)沸騰過程的模擬研究中均取0.1 s?1，顯然，熱管模擬直接套用該模型，穩(wěn)態(tài)階段蒸發(fā)和冷凝過程的熱質(zhì)傳遞量將難以達到平衡。因此，本文利用UDF功能對c進行修正，促使蒸發(fā)和冷凝過程趨于平衡，修正公式見式(14)，模擬中e保持為0.1不變，c初始值取0.1。

式中：下標表示時間步長數(shù)；e和c分別為蒸發(fā)和冷凝過程的質(zhì)量傳遞總量。

熱管穩(wěn)態(tài)運行階段，蒸發(fā)段工質(zhì)處于核態(tài)沸騰狀態(tài)，沸騰過熱度可由下式[16]求得：

考慮到管壁的熱傳導(dǎo)過程，則sat取值為

式中：e為實測蒸發(fā)段均溫；o和i分別為熱管外徑和內(nèi)徑；為管壁熱導(dǎo)率；e為蒸發(fā)段長度；為輸入功率。

1.2.3 CSF模型

表面張力是液體表面分子引力不均衡而產(chǎn)生的沿表面作用的張力，可使汽泡或液滴的表面積趨于最小，降低表面自由能。因此，表面張力對汽泡運動以及工質(zhì)的潤濕特性具有重要影響。本文通過引入CSF模 型[15]修改動量源項來模擬表面張力作用，該源項如下式所示：

式中：為表面曲率。

1.3 邊界條件及求解策略

根據(jù)文獻[13]的實驗條件，熱管蒸發(fā)段為恒功率電加熱，冷凝段為水冷。故模擬時蒸發(fā)段施加恒熱流密度進行加熱，冷凝段采用水的對流換熱，輸入熱流密度、冷卻水溫度及對流換熱系數(shù)的取值均由實驗條件而定，絕熱段設(shè)為絕熱壁面，即熱流密度為0。

能量方程與動量方程均以二階迎風格式離散，相體積率和壓力的插值分別采用幾何重構(gòu)法和PRESTO法，速度與壓力耦合關(guān)系由SIMPE算法處理。

計算域初始化溫度和壓力分別為298 K和1.72 kPa，熱管充液率R(工質(zhì)體積與蒸發(fā)段容積之比)通過patch功能設(shè)置。動量方程與連續(xù)性方程的殘差收斂標準設(shè)為10?4，能量方程為10?7。計算初期，時間步長取1×10?5s，隨模擬時間增加而逐步提高至5×10?4s。當熱管蒸發(fā)段輸出功率與輸入功率達到平衡時認為熱管達到穩(wěn)態(tài)，計算完成。

1.4 網(wǎng)格獨立性驗證

為得到網(wǎng)格獨立性的解，本文對R為0.6的熱管進行考核，計算網(wǎng)格數(shù)分別取25 670，39 640，58 410和85 420。網(wǎng)格劃分采用四邊形結(jié)構(gòu)化網(wǎng)格，熱管內(nèi)壁邊界層經(jīng)加密處理，如圖1所示。圖2所示為熱管穩(wěn)態(tài)工作溫度w在不同網(wǎng)格數(shù)下的模擬結(jié)果。從圖2可知：w隨網(wǎng)格數(shù)增多趨于定值，當網(wǎng)格數(shù)目高于 58 410時，網(wǎng)格數(shù)對w的影響可忽略，因此，本文計算網(wǎng)格數(shù)取58 410。

圖1 網(wǎng)格劃分示意圖

圖2 網(wǎng)格獨立性考核

2 可視化實驗

建立的熱管可視化實驗平臺如圖3所示。熱管由全玻璃制成，其尺寸與文獻[13]中相同。實驗中熱管蒸發(fā)段以水浴加熱，冷凝段為玻璃夾套水冷，絕熱段和水冷夾套外用保溫棉絕熱。通過加熱溫度和冷卻條件調(diào)整熱管的工作溫度，使其運行狀態(tài)接近文獻[13]所述工況，并利用高速相機(D7000 日本尼康)記錄穩(wěn)態(tài)工況下熱管蒸發(fā)段的沸騰現(xiàn)象，與模擬結(jié)果對比 驗證。

為確保良好的可視化觀測效果，將可視化熱管蒸發(fā)段浸沒在圓柱形玻璃水槽水浴中，但這樣增加了有效壁溫分布測定的難度，因此，本文未對該可視化熱管的溫度分布進行模擬對比。

圖3 可視化熱管及實驗平臺

3 結(jié)果與討論

3.1 模型驗證

3.1.1 溫度分布驗證

表1列出加熱功率670 W條件下不同充液率熱管蒸發(fā)段、絕熱段和冷凝段穩(wěn)態(tài)均溫的實驗值[13]與模擬值。從表1可以看出：冷凝段模擬值與實驗值的偏差最大，最大偏差為10.2 K(相對偏差2.9%)，原因在于實驗所測溫度為冷凝夾套外壁溫度，而模擬值為熱管管壁溫度。但是，蒸發(fā)段和絕熱段的模擬與實驗值的相對偏差均在2.1%以內(nèi)，遠低于FADHL等[11]的模擬偏差，可見相變模型的改進對模擬結(jié)果準確性的改善作用明顯，溫度分布模擬值與實測值的吻合度更高。

3.1.2 流型定性驗證

在熱管管壁溫度分布定量驗證的基礎(chǔ)上，進一步對相變流動行為作可視化定性驗證。由表1可知：在文獻[13]測試工況下，R為0.3或0.9時熱管工作溫度均大于373 K，受水浴加熱溫度限制，本文僅針對充液率R為0.6的熱管進行了可視化實驗和模擬流型的定性驗證。圖4所示為熱管(R=0.6)蒸發(fā)段穩(wěn)態(tài)流型CFD模擬結(jié)果與可視化現(xiàn)象的對比，模擬時加熱功率為670 W。從圖4可知：可視化實驗直觀地展現(xiàn)了蒸發(fā)段汽泡長大脫離壁面后的合并、失穩(wěn)和破碎現(xiàn)象，即流型從泡狀流、栓塞流向環(huán)狀流的演變過程，熱管CFD模型較準確地預(yù)測了該流型轉(zhuǎn)變，與實驗結(jié)果基本相符。模擬偏差產(chǎn)生的主要原因可能是熱管CFD模擬蒸發(fā)段采用的恒熱流密度邊界條件與可視化實驗中的對流換熱存在差異。

表1 熱管各段均溫實驗值與模擬值

圖4 蒸發(fā)段流型對比(FR=0.6)

3.2 充液率對熱管相變傳熱行為的影響

3.2.1 汽液兩相分布

工作溫度可定性反映熱管的運行狀態(tài)，在相同的傳輸功率下，工作溫度越低，意味著熱管內(nèi)工質(zhì)的傳熱循環(huán)效率越高，熱管整體傳熱性能越好。表1所示熱管實測溫度顯示：充液率對w有顯著影響，R為0.6時w最低，充液率的增大或降低均會引起w的升高，因此，測試條件下熱管最佳充液率為0.6。

熱管內(nèi)工質(zhì)的流動傳熱行為與充液率密切相關(guān)，為進一步揭示充液率對熱管傳熱性能的影響，有必要對管內(nèi)流型進行分析。圖5所示為穩(wěn)態(tài)工況下，不同充液率熱管(R為0.3，0.6和0.9)內(nèi)兩相分布的模擬結(jié)果。從圖5可以看出：在不同充液率下，熱管冷凝段均為膜狀冷凝換熱，充液率的改變主要影響蒸發(fā)段的傳熱行為。對于R為0.3的熱管，液池區(qū)域為核態(tài)沸騰，但液池僅占據(jù)蒸發(fā)段底部小部分區(qū)域，回流液膜在蒸發(fā)段上部出現(xiàn)了間斷。R為0.6條件下，液池略高于蒸發(fā)段出口，蒸發(fā)段內(nèi)壁完全潤濕，有效沸騰換熱面積增大，熱管處于較佳運行狀態(tài)，工作溫度最低。而當R達到0.9時，蒸發(fā)段大量上浮的汽泡將工質(zhì)液面托至絕熱段，同時，蒸發(fā)段還觀察到成片的膜狀汽泡，表明蒸發(fā)段傳熱條件惡化。

圖5 不同充液率熱管兩相分布

3.2.2 溫度和速度分布

通過對熱管溫度分布和工質(zhì)流動的分析可進一步了解充液率對熱管相變傳熱行為的影響。不同充液率下熱管管壁溫度分布以及管內(nèi)軸線溫度和工質(zhì)流速分布如圖6所示。從圖6(a)可以看出：充液率較低時(R為0.3)，熱管蒸發(fā)段上部形成了局部高溫區(qū)，并且，管內(nèi)軸向均溫性較差。由圖5可知：由于充液率較低，熱管蒸發(fā)段上部無液池浸潤，為不連續(xù)的液膜蒸發(fā)傳熱，換熱系數(shù)遠低于液池區(qū)域的核態(tài)沸騰傳熱，致使局部高溫的形成。蒸發(fā)段的局部干涸又導(dǎo)致蒸汽產(chǎn)生一定的過熱度以及蒸發(fā)段出口蒸汽溫度的振蕩，造成管內(nèi)軸向均溫性較差。

由圖6(b)和(c)可知：增大充液率有利于改善熱管管壁和管內(nèi)的溫度分布，然而，充液率的增大為汽泡的發(fā)育提供了更多的空間，汽泡脫離液面的速度也急劇增大。相比于R為0.3的熱管，R為0.6和0.9時熱管蒸發(fā)段液面處工質(zhì)流速從1.5 m/s增加至4.5 m/s左右。另一方面，R為0.6，熱管內(nèi)工質(zhì)最大流速點恰好處于蒸發(fā)段出口，而當R增大至0.9時，液池液面則提升至絕熱段，絕熱段流速較大的工質(zhì)團將阻礙冷凝液及時回流至蒸發(fā)段，這也是其工作溫度相對較高的原因之一。

FR：(a) 0.3；(b) 0.6；(c) 0.9

從圖6還可發(fā)現(xiàn)：不同充液率條件下熱管冷凝段上部均出現(xiàn)了溫度回升現(xiàn)象，這與ALIZADEHDAKHEL等[10]在實驗中觀察到的現(xiàn)象相符，其原因在于，冷凝段頂端冷凝液膜最薄，熱阻較低，液膜在下降過程中逐漸增厚，熱阻增大，故冷凝段頂端換熱系數(shù)高于其余區(qū)域，管壁溫度出現(xiàn)回升。

4 結(jié)論

1) 熱管管壁溫度分布模擬值與實驗值具有較好的一致性，各段均溫的最大偏差為10.2 K，相對誤差2.9%，采用改進的相變模型，模擬結(jié)果準確度更高，同時，模擬得到蒸發(fā)段工質(zhì)流型從泡狀流、栓塞流向環(huán)狀流的轉(zhuǎn)變過程，與可視化實驗定性相符，模型能夠較好地預(yù)測熱管內(nèi)工質(zhì)沸騰與冷凝回流的傳熱行為。

2) 熱管冷凝段均為膜狀冷凝，充液率的改變主要影響熱管蒸發(fā)段的相變傳熱行為：充液率過小易導(dǎo)致蒸發(fā)段局部干涸過熱，均溫性惡化，過大則形成膜狀沸騰，抬高工質(zhì)液面，阻礙冷凝液回流。推薦的熱管充液率為0.6，其液面高度與蒸發(fā)段高度持平，傳熱性能優(yōu)良。

[1] 莊駿, 張紅. 熱管技術(shù)及其工程應(yīng)用[M]. 北京: 化工工業(yè)出版社, 2000: 67.ZHUANG Jun, ZHANG Hong. Heat pipe technology and engineering application[M]. Beijing: Chemical Industry Press, 2000: 67.

[2] LIU Zhenhua, HU Renlin, LU Lin, et al. Thermal performance of an open thermosyphon using nanofluid for evacuated tubular high temperature air solar collector[J]. Energy Conversion and Management, 2013, 73(5): 135?143.

[3] 王子龍, 張華, 張海濤. 聚光光伏系統(tǒng)熱管散熱器運行特性的優(yōu)化[J]. 太陽能學報, 2012(6): 986?992. WANG Zilong, ZHANG Hua, ZHANG Haitao. Optimization of operation characteristics for finned heat pipe radiator in concentrating photovoltaic system[J]. Acta Energiae Solaris Sinica, 2012(6): 986?992.

[4] JOUHARA H, MERCHANT H. Experimental investigation of a thermosyphon based heat exchanger used in energy efficient air handling units[J]. Energy, 2012, 39(1): 82?89.

[5] CIE?LI?SKI J T, FIUK A. Heat transfer characteristics of a two-phase thermosyphon heat exchanger[J]. Applied Thermal Engineering, 2013, 51(S1/S2): 112?118.

[6] TSAI T E, WU H H, CHANG C C, et al. Two-phase closed thermosyphon vapor-chamber system for electronic cooling[J]. International Communications in Heat and Mass Transfer, 2010, 37(5): 484?489.

[7] 辛公明, 張魯生, 王鑫煜, 等. 自濕潤流體內(nèi)螺紋重力熱管傳熱特性[J]. 工程熱物理學報, 2013, 34(12): 2363?2365. XING Gongming, ZHANG Lusheng, WANG Xinyu, et al. Heat transfer characteristics of internal helical micro fin gravity heat pipe with self-rewetting fluid[J]. Journal of Engineering Thermophysics, 2013, 34(12): 2363?2365.

[8] 曲偉, 王煥光. 高溫及超高溫熱管的相容性和傳熱性能[J]. 化工學報, 2011, 62(S1): 77?81. QU Wei, WANG Huangguang. Compatibility and heat transfer of high and super high temperature heat pipes[J]. Journal of Chemical Industry and Engineering (China), 2011, 62(S1): 77?81.

[9] SCHEPPER S C K D, HEYNDERICKX G J, MARIN G B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker[J]. Computers and Chemical Engineering, 2009, 33(1): 122?132.

[10] ALIZADEHDAKHEL A, RAHIMI M, ALSAIRAFI A A. CFD modeling of flow and heat transfer in a thermosyphon[J]. International Communications in Heat and Mass Transfer, 2010, 37(3): 312?318.

[11] FADHL B, WROBEL L C, JOUHARA H. Numerical modelling of the temperature distribution in a two-phase closed thermosyphon[J]. Applied Thermal Engineering, 2013, 60(1/2): 122?131.

[12] ASMAIE L, HAGHSHENASFARD M, MEHRABANI- ZEINABAD A, et al. Thermal performance analysis of nanofluids in a thermosyphon heat pipe using CFD modeling[J]. Heat Mass Transfer, 2013, 49(5): 667?678.

[13] NOIE S H. Heat transfer characteristics of a two-phase closed thermosyphon[J]. Applied Thermal Engineering, 2005, 25(4): 495?506.

[14] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201?225.

[15] BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J]. Journal of Computational Physics, 1992, 100(2): 335?354

[16] REAY D A, KEW P A. Heat pipes[M]. 15th ed. Oxford: Butterworth-Heinemann, 2006: 65.

(編輯 陳愛華)

CFD modeling of phase change heat transfer behaviors in thermosyphons

WANG Xiaoyuan, ZHU Yuezhao, CHEN Haijun, WANG Yinfeng, FAN Hongtu

(School of Mechanical and Power Engineering, Nanjing Tech University, Nanjing 211816, China)

A computational fluid dynamic (CFD) model based on the volume of fluid (VOF) multiphase model and an improved phase change model were developed to simulate the phase change heat transfer process in thermosyphons at steady state, and the modeling results were validated by measuring wall temperature distribution and visualization experiments. Besides, the influence of filling ratio (FRs) on the evaporation-condensation behavior in thermosyphons was analyzed. The results show that this CFD model can reproduce both the vaporization and condensation processes well. The predicted temperature distribution show better consistency with experimental results using the improved phase change model, the maximum deviation of average temperature at different sections being 10.2 K with a relative deviation of 2.9%. Furthermore, the predicted flow regimes are generally in agreement with the observed phenomenon in visualization testing. Filmwise condensation is the main heat transfer mode at condenser area in thermosyphons which is not affected by FRs. However, the variation of FRs considerably affects the heat transfer behavior at evaporator and thus changes the heat transfer performance of thermosyphons. This study is conducive for understanding the operating mechanism of thermosyphons and meanwhile provides some guidance for its CFD modeling.

thermosyphon; phase change heat transfer; computational fluid dynamic (CFD); visualization

10.11817/j.issn.1672-7207.2017.05.035

TK172.4

A

1672?7207(2017)05?1391?07

2016?07?21；

2016?09?09

國家自然科學基金資助項目(51276086)；國家科技支撐計劃項目(2014BAJ01B06) (Project(51276086) supported by the National Natural Science Foundation of China; Project(2014BAJ01B06) supported by the National Science and Technology Support Program of China)

朱躍釗，博士，教授，博士生導(dǎo)師，從事熱科學與工程等研究；E-mail: zyz@njtech.edu.cn