方肋微通道內(nèi)流動(dòng)沸騰的氣泡動(dòng)態(tài)與傳熱特性分析

申宇，潘振海，吳慧英

（上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院，上海200240）

隨著近年來(lái)微電子技術(shù)和MEMS工藝在諸多工業(yè)領(lǐng)域的飛速發(fā)展，電子器件在尺寸縮小的同時(shí)，功耗卻不斷攀升。這導(dǎo)致電子設(shè)備表面的熱通量激增，嚴(yán)重影響系統(tǒng)性能甚至直接燒毀電子器件[1]。帶肋微通道熱沉的流動(dòng)沸騰換熱在高熱通量芯片熱控制方面具有顯著優(yōu)勢(shì)，得到了廣泛關(guān)注[2-10]：與常規(guī)通道不同，帶肋微通道的換熱比表面積大，結(jié)構(gòu)緊湊，便于封裝；流動(dòng)沸騰利用液體相變潛熱，換熱效率高，并能維持較好的壁面均溫性。此外，相較于平滑微通道，其內(nèi)部的微肋結(jié)構(gòu)還具有強(qiáng)化流體擾動(dòng)和增加氣化核心的雙重功效，進(jìn)一步增強(qiáng)對(duì)流換熱性能并降低沸騰起始溫度。

國(guó)內(nèi)外諸多學(xué)者對(duì)帶肋微通道內(nèi)獨(dú)特的氣液兩相流動(dòng)和相變傳熱機(jī)制展開(kāi)了研究。Lie等[6]通過(guò)實(shí)驗(yàn)對(duì)比了制冷劑FC-72在帶肋微通道內(nèi)兩相流和單相流的換熱效果，發(fā)現(xiàn)兩相流動(dòng)的換熱系數(shù)遠(yuǎn)大于單相流。Krishnamurthy和Peles[7]進(jìn)行了并排帶肋微通道內(nèi)的流動(dòng)沸騰實(shí)驗(yàn)，在通道內(nèi)沿流動(dòng)方向分別觀察到泡狀流、泡狀-彈狀混合流、環(huán)狀流等不同流型，并且發(fā)現(xiàn)帶肋微通道的流動(dòng)沸騰換熱效率明顯高于光滑微通道。Guo 等[8]分別研究了光滑表面和柱狀微結(jié)構(gòu)換熱面的流動(dòng)-噴射復(fù)合式強(qiáng)化沸騰換熱。結(jié)果表明得益于柱間的微對(duì)流運(yùn)動(dòng)，相比光滑表面，相同工況下柱狀微結(jié)構(gòu)的臨界熱通量（CHF）更高，沸騰穩(wěn)定性得到提升。Qu 等[10]研究了進(jìn)口流體過(guò)冷度對(duì)微通道內(nèi)流動(dòng)沸騰的影響，發(fā)現(xiàn)當(dāng)過(guò)冷度增大，低含氣率區(qū)域的壁面?zhèn)鳠嵯禂?shù)會(huì)提高，但是在高含氣率區(qū)域傳熱變化不明顯。為了深入理解流動(dòng)沸騰換熱的諸多物理機(jī)制，學(xué)者們[11-19]采用數(shù)值模擬的方法對(duì)該過(guò)程進(jìn)行研究。楊朝初等[11]模擬了微氣泡的生長(zhǎng)過(guò)程，發(fā)現(xiàn)氣泡一旦成核后會(huì)迅速吸熱膨脹，并在氣泡與通道壁面之間形成薄液膜。彈狀流是微通道內(nèi)非常重要的兩相流型，具有良好的傳熱傳質(zhì)特性[12]。Magnini等[13]基于Hertz-Knudsen 理論和高度函數(shù)法構(gòu)建了高精度的氣液相變模型，對(duì)微圓管內(nèi)彈狀流氣泡的流動(dòng)沸騰過(guò)程進(jìn)行模擬。結(jié)果表明微通道內(nèi)兩相流的壁面?zhèn)鳠嵯禂?shù)比相同條件的單相流高出30%，并發(fā)現(xiàn)薄液膜的蒸發(fā)導(dǎo)熱是壁面?zhèn)鳠釓?qiáng)化的主要原因。徐進(jìn)良等[14]模擬了高熱通量下微圓管內(nèi)觸發(fā)種子氣泡的沸騰傳熱，并分析了不同氣泡觸發(fā)頻率對(duì)壁面冷卻的影響。從以上文獻(xiàn)可以看出，帶肋微通道具有良好的沸騰換熱性能，但是對(duì)流動(dòng)沸騰過(guò)程中所涉及的氣泡運(yùn)動(dòng)規(guī)律以及肋表面的換熱特性還需進(jìn)一步研究。

因此本文采用數(shù)值模擬的方法，構(gòu)建氣-液流動(dòng)和相變模型，對(duì)微通道內(nèi)氣泡繞流加熱方肋這一基本的流動(dòng)沸騰問(wèn)題進(jìn)行研究。以氣泡行為、流場(chǎng)演變以及方肋各表面的傳熱特性變化作為主要研究對(duì)象，通過(guò)分析初始?xì)馀蒹w積及入口雷諾數(shù)Re 等不同參數(shù)的影響，探究帶肋微通道內(nèi)沸騰氣-液兩相流動(dòng)及相變傳熱傳質(zhì)的規(guī)律，進(jìn)而為實(shí)驗(yàn)研究提供指導(dǎo)。

1 模型建立

1.1 物理模型

本文構(gòu)造一個(gè)二維微通道物理模型，如圖1所示。微通道長(zhǎng)度L=4mm，寬度D=0.2mm，通道壁面絕熱。邊長(zhǎng)B=0.08mm 的方柱設(shè)置在通道正中心，方柱壁面為恒定熱通量q=5×104W/m2加熱。各壁面均為無(wú)滑移速度條件。微通道入口通入飽和R113制冷劑（Tsat=323.15K），其標(biāo)準(zhǔn)大氣壓下各相物性參數(shù)見(jiàn)表1。入口設(shè)置充分發(fā)展速度條件，出口為標(biāo)準(zhǔn)大氣壓出口條件，并在距離入口0.5mm，靠近通道上壁面的位置初始化一個(gè)氣泡。

以D為特征長(zhǎng)度，定義入口雷諾數(shù)Re[式(1)]。

式中，u為入口平均速度；μl為液體動(dòng)力黏度。定義局部傳熱系數(shù)hx[式(2)]。

圖1 幾何模型和邊界條件

表1 一個(gè)標(biāo)準(zhǔn)大氣壓下R113飽和狀態(tài)的物性參數(shù)

式中，Tw,x為方柱表面局部壁溫。

定義局部努賽爾數(shù)Nux[式(3)]。

式中，λl為液相導(dǎo)熱率。

1.2 網(wǎng)格劃分

本文采用ICEM 軟件劃分網(wǎng)格，部分計(jì)算區(qū)域網(wǎng)格如圖2所示。為了保證計(jì)算域內(nèi)氣液兩相流動(dòng)與沸騰傳熱特性的準(zhǔn)確性，對(duì)微通道壁面以及加熱方柱周圍進(jìn)行了網(wǎng)格加密處理，以確保對(duì)薄液膜的刻畫(huà)。為了驗(yàn)證網(wǎng)格的無(wú)關(guān)性，在此基礎(chǔ)上將原有網(wǎng)格數(shù)量提高4倍，方柱換熱系數(shù)結(jié)果與原網(wǎng)格算例相比差異小于1.5%，因此可認(rèn)為結(jié)果是網(wǎng)格無(wú)關(guān)解。

1.3 控制方程和算法設(shè)置

本研究中的數(shù)值模型基于如下假設(shè)：

（1）氣液相變溫度恒定為一個(gè)標(biāo)準(zhǔn)大氣壓下的飽和溫度；

（2）氣相和液相的熱物性均設(shè)置為該飽和溫度下的常值；

（3）不考慮液體黏性加熱；

（4）忽略重力作用。

以上假設(shè)的根據(jù)是：①微通道內(nèi)壓降小于5kPa，其對(duì)應(yīng)的液體相變溫度變化小于5%[21]；②最大過(guò)熱度低于13K，其對(duì)應(yīng)的氣相和液相的熱物性變化均小于5%[21]；③最大Re 下，黏性耗散對(duì)流體的加熱僅為壁面加熱的0.1%[22]；④通道封閉數(shù)Co大于5[Co=(σ/gΔρD2)1/2]，重力作用被抑制[23]。

本文使用的流體體積函數(shù)（VOF）模型是一種在固定歐拉網(wǎng)格下追蹤氣液交界面的方法。氣相和液相用網(wǎng)格體積分?jǐn)?shù)a 來(lái)表示：a=1 為氣相，a=0為液相。氣液界面處0＜a＜1，物性為各相物性體積分?jǐn)?shù)的加權(quán)平均值。氣相及液相的質(zhì)量方程如式(4)、式(5)所示，氣液兩相共用的動(dòng)量、能量守恒方程如式(6)、式(7)所示。

式(6)中的Fs表示作用于氣液界面的表面張力，通過(guò)連續(xù)界面力（continuum surface force，CSF）模型將表面力轉(zhuǎn)化為體積力[式(8)]。

式(7)中的能量源項(xiàng)Sh基于Pan 等[18]提出的“飽和界面體積”相變模型，該模型中認(rèn)為蒸發(fā)過(guò)程中氣液界面處維持在飽和溫度，因此計(jì)算可以不依賴于經(jīng)驗(yàn)系數(shù)，其表達(dá)式為式(9)。

對(duì)應(yīng)的式(4)、式(5)的質(zhì)量源項(xiàng)表達(dá)式為式(10)。

本文使用基于有限體積法的Ansys Fluent 18.0離散型非穩(wěn)態(tài)求解器，流動(dòng)模型為層流模型。流場(chǎng)采用PISO（pressure implicit splitting of operator）算法求解耦合的壓力-速度方程，對(duì)時(shí)間項(xiàng)采用隱式格式離散。標(biāo)量梯度的求解采用“基于節(jié)點(diǎn)的格林-高斯”法（Green-Gauss node-based method），采用VOF方法中的“結(jié)構(gòu)重構(gòu)”方法（geo-reconstruct）追蹤氣液界面。壓力插值采用PRESTO 算法（pressure staggering option），動(dòng)量和能量方程的對(duì)流項(xiàng)均采用三階迎風(fēng)格式進(jìn)行離散。為了避免數(shù)值振蕩，通過(guò)設(shè)置CFL 條件實(shí)現(xiàn)對(duì)時(shí)間步長(zhǎng)的控制，Courant 數(shù)設(shè)置為0.05，模擬過(guò)程中最大的時(shí)間步長(zhǎng)為10-7s。質(zhì)量、動(dòng)量和能量方程的殘差分別控制在10-6、10-8和10-10以下。

圖2 計(jì)算網(wǎng)格劃分

1.4 模型可靠性驗(yàn)證

為驗(yàn)證計(jì)算模型的精確性和可靠性，本文首先用R113工質(zhì)對(duì)經(jīng)典的一維Stefan相變問(wèn)題[24]進(jìn)行模擬。圖3 給出了氣液界面位置隨時(shí)間的變化曲線，可以看出，本文的模擬結(jié)果與解析解吻合良好。

圖3 相界面位置隨時(shí)間變化

其次，Han和Shikazono[25]實(shí)驗(yàn)研究了微通道內(nèi)彈狀流的薄液膜厚度變化，本文針對(duì)液膜厚度與文獻(xiàn)結(jié)果進(jìn)行對(duì)比。在Ca=μu/σ=0.0029 時(shí)，測(cè)量值為10.2μm，模擬值為10.7μm；Ca=0.044 時(shí)，測(cè)量值為63.2μm，模擬值為66μm，模擬值與實(shí)驗(yàn)值比較接近。

最后，本文模擬了圓管微通道內(nèi)彈狀氣泡的流動(dòng)沸騰，并與文獻(xiàn)[13]中基于Hertz-Knudsen 理論的蒸發(fā)模型進(jìn)行對(duì)比，結(jié)果如圖4、圖5 所示。可以看出本文的氣泡位置和局部傳熱系數(shù)均與文獻(xiàn)結(jié)果吻合較好，說(shuō)明本文采用的模型和方法能夠準(zhǔn)確模擬微通道內(nèi)的相變流動(dòng)問(wèn)題。

圖4 氣泡前端位置隨時(shí)間變化

圖5 壁面?zhèn)鳠嵯禂?shù)沿軸向分布

2 計(jì)算結(jié)果及分析

2.1 氣泡動(dòng)態(tài)及流場(chǎng)變化

基于上述模型，本文通過(guò)改變?nèi)肟诶字Z數(shù)（Re=60～360）和氣泡初始大小（d=0.04mm、0.08mm、0.12mm）進(jìn)行了一系列模擬，以研究流體流速和氣泡體積對(duì)方肋微通道內(nèi)兩相流動(dòng)和強(qiáng)化傳熱的影響。

圖6 顯示的是初始直徑d=0.08mm 的氣泡在不同Re 下，流過(guò)加熱方柱時(shí)微通道內(nèi)速度場(chǎng)、溫度場(chǎng)以及氣泡形狀的演變過(guò)程。其中，氣泡的輪廓在速度云圖中用黑線標(biāo)識(shí)，在溫度云圖中用白線標(biāo)識(shí)，帶箭頭的線為流線。

從圖6左側(cè)圖片可以看到隨著Re增大，受流體慣性力作用，氣泡變形程度加劇。同時(shí)分析速度場(chǎng)和流線圖，Re 為60 和120 時(shí)方柱后表面的尾渦區(qū)內(nèi)有兩個(gè)對(duì)稱的渦旋；而當(dāng)Re提高到240后，尾渦區(qū)內(nèi)的渦旋開(kāi)始周期性地脫落，增強(qiáng)了下游流體的擾動(dòng)。當(dāng)氣泡從方柱與壁面的狹縫穿過(guò)時(shí)，氣液界面掃掠方柱上表面的過(guò)熱液體，形成一層薄液膜，擠壓熱邊界層。對(duì)比氣泡進(jìn)出方柱前后，尾渦區(qū)的流場(chǎng)和溫度場(chǎng)發(fā)生了明顯改變。

圖7 分別是四種Re 下，氣泡的當(dāng)量直徑deq隨氣泡質(zhì)心軸向位置xc變化的曲線。其中，氣泡的當(dāng)量直徑deq定義為同等體積下球形氣泡的直徑；氣泡質(zhì)心的軸向位置xc可用式(11)計(jì)算。

式中，N為計(jì)算域的網(wǎng)格數(shù)，xi是第i個(gè)網(wǎng)格的軸向坐標(biāo)，ai和Vi分別為該網(wǎng)格的氣相體積分?jǐn)?shù)與網(wǎng)格體積。

圖6 氣泡在4種不同Re下，通過(guò)加熱方柱過(guò)程中三個(gè)位置的絕對(duì)速度及溫度云圖

圖7 不同Re下氣泡當(dāng)量直徑deq隨質(zhì)心xc位置的變化

從圖7可以看出，當(dāng)氣泡到達(dá)方柱中心時(shí)（xc=2mm），氣泡體積開(kāi)始陡增，說(shuō)明此時(shí)相變蒸發(fā)已經(jīng)開(kāi)始。從曲線斜率的變化可以看出，受過(guò)熱液體分布的影響，氣泡在方柱和尾渦區(qū)域（xc=2～2.5mm） 體積迅速增長(zhǎng)，而到通道下游（xc＞2.5mm）體積增長(zhǎng)則不明顯。此外，氣泡在低Re條件下的增長(zhǎng)速度要遠(yuǎn)遠(yuǎn)大于高Re 的情況，尤其是在尾渦區(qū)內(nèi)。出現(xiàn)該現(xiàn)象的原因主要有兩點(diǎn)：一是高Re 下，氣泡流速快，在過(guò)熱區(qū)停留時(shí)間短，蒸發(fā)吸熱過(guò)程不充分；二是在高Re 下，流場(chǎng)以渦旋脫落的方式將尾渦區(qū)過(guò)熱液體帶向通道下游，對(duì)冷熱液體進(jìn)行混合，降低了尾渦區(qū)液體的過(guò)熱度，從而抑制了該區(qū)域的相變蒸發(fā)。

2.2 換熱效果分析

2.2.1 入口雷諾數(shù)Re對(duì)方柱傳熱的影響

圖8對(duì)比了4種Re下，方柱四個(gè)加熱面平均傳熱強(qiáng)化率隨時(shí)間的變化。其中，平均傳熱強(qiáng)化率定義為(htp-hsp)/hsp，反映了兩相流下表面平均傳熱系數(shù)htp與相同條件下單相流hsp的相對(duì)大小。定義量綱為1 時(shí)間τ=t/(B/u)，當(dāng)氣泡前端靠近方柱時(shí)τ 設(shè)定為0 時(shí)刻。從圖中可以看出，當(dāng)氣泡穿過(guò)方柱時(shí)，方柱4個(gè)表面的傳熱均受到不同程度影響。其中上表面的傳熱強(qiáng)化最為顯著：在Re=60下，最大的換熱系數(shù)超過(guò)單相流的4.5 倍。這是因?yàn)楫?dāng)氣泡穿過(guò)方柱時(shí)，氣泡在表面張力的作用下擠壓周圍液體，在方柱的表面和氣泡界面之間形成一層薄液膜（圖4）。此時(shí)薄液膜內(nèi)部流動(dòng)緩慢，薄液膜的蒸發(fā)導(dǎo)熱成為了主要換熱方式，其熱阻較小，強(qiáng)化傳熱顯著[15]。而隨著Re 增大，流體慣性力開(kāi)始占據(jù)主導(dǎo)，液膜厚度變厚，導(dǎo)熱熱阻增大，導(dǎo)致強(qiáng)化傳熱效果迅速下降。方柱前表面和后表面也形成部分液膜，因此具有類似變化規(guī)律。而對(duì)于方柱下表面，雖然未與氣泡直接形成薄液膜，但由于氣泡堵塞上通道而使下通道的液體流動(dòng)加快，增強(qiáng)了下表面的對(duì)流傳熱強(qiáng)度。

圖8 方柱各表面兩相流傳熱強(qiáng)化率隨量綱為1時(shí)間τ的變化

進(jìn)一步對(duì)方柱4個(gè)表面的傳熱強(qiáng)化率做積分平均，圖9(a)比較了4種Re下方柱整體傳熱強(qiáng)化率隨時(shí)間的變化；對(duì)氣泡前端靠近方柱直到氣泡尾端離開(kāi)方柱這段時(shí)間進(jìn)行積分平均，圖9(b)統(tǒng)計(jì)了方柱的時(shí)均-整體傳熱強(qiáng)化率與Re的關(guān)系。從這兩個(gè)圖都可以看出，在氣泡流經(jīng)方柱的過(guò)程中，Re越小，氣泡流動(dòng)沸騰對(duì)方柱整體傳熱強(qiáng)化效果越明顯。

圖9 氣泡流經(jīng)方柱過(guò)程方柱整體換熱強(qiáng)化率的變化

2.2.2 氣泡初始直徑對(duì)方柱傳熱的影響

圖10 比較了Re=120 下3 種量綱為1 初始直徑的氣泡（= deq/D），其質(zhì)心位置xc恰好位于方柱中心位置（即xc=4mm）時(shí)氣泡的形狀、方柱壁面過(guò)熱度以及對(duì)應(yīng)的局部努賽爾數(shù)Nux的分布情況。

圖10 不同初始體積的氣泡穿過(guò)方柱的對(duì)比

3 結(jié)論

本文基于VOF 方法，采用“飽和界面”相變模型，以R113 為工質(zhì)對(duì)帶肋微通道內(nèi)氣泡流動(dòng)沸騰進(jìn)行了二維數(shù)值模擬，分析了氣泡體積和入口雷諾數(shù)Re 對(duì)氣泡增長(zhǎng)速率、液膜厚度和方柱表面?zhèn)鳠岬戎匾獋鳠醾髻|(zhì)參數(shù)的影響規(guī)律。得到如下結(jié)論。

（1）氣泡穿過(guò)方柱時(shí)，在氣泡界面和方柱表面之間形成一層薄液膜，擠壓溫度邊界層；氣泡對(duì)尾渦區(qū)的流動(dòng)結(jié)構(gòu)存在擾動(dòng)作用。

（2）氣泡的蒸發(fā)相變主要發(fā)生在方柱和尾渦區(qū)內(nèi)，且Re越小，氣泡的體積增長(zhǎng)越快。

（3）薄液膜的蒸發(fā)導(dǎo)熱是傳熱強(qiáng)化的主導(dǎo)機(jī)制，氣泡的流動(dòng)沸騰對(duì)方柱表面換熱有明顯的強(qiáng)化作用；對(duì)于相同體積的氣泡，隨著Re 增加，液膜厚度逐漸變厚，強(qiáng)化傳熱作用迅速減弱。

（4）對(duì)于相同Re，氣泡體積越大，所形成的薄液膜區(qū)覆蓋面積就越大，液膜蒸發(fā)帶來(lái)的強(qiáng)化換熱也越顯著；而小氣泡的兩相流換熱效果同單相流相比并無(wú)明顯差異。

符號(hào)說(shuō)明

B—— 方柱邊長(zhǎng)，mm

Co—— 通道封閉數(shù)，(σ/gΔρD2)1/2

Ca—— 毛細(xì)數(shù)，μu/σ

cp—— 比定壓熱容，J/(kg·K)

D—— 通道寬度，mm

d—— 氣泡直徑，mm

Fs—— 體積表面張力，N/m3

g—— 重力加速度，m/s2

hfg—— 氣化潛熱，kJ/kg

hx—— 對(duì)流換熱系數(shù)，W/(m2·K)

L—— 通道長(zhǎng)度，mm

Nux—— 局部努賽爾數(shù)

n—— 單位法向量

p—— 壓強(qiáng)，Pa

q—— 熱通量，W/m2

Re—— 入口雷諾數(shù)，Re=ρluD/μl

Sm—— 質(zhì)量源項(xiàng)，kg/(m3·s)

Sh—— 能量源項(xiàng)，W/m3

t—— 時(shí)間，s

U—— 速度矢量，m/s

u—— 軸向速度，m/s

a—— 體積分?jǐn)?shù)

λ—— 熱導(dǎo)率，W/(m·K)

μ—— 動(dòng)力黏度，Pa·s

ρ—— 密度，kg/m3

σ—— 表面張力系數(shù)，N/m

τ—— 量綱為1時(shí)間，tu/B

x—— 軸向坐標(biāo)，mm

下角標(biāo)

c—— 質(zhì)心

eq—— 當(dāng)量

l—— 液相

sat—— 飽和狀態(tài)

sp—— 單相流

tp—— 兩相流

w—— 壁面

v—— 氣相