非均勻熱流下細小通道內相變微膠囊懸浮液傳熱特性

李梓龍 王智彬 賈莉斯 陳 穎 莫松平

(廣東工業大學材料與能源學院 廣東省功能軟凝聚態物質重點實驗室 廣州 510006)

電子科技快速發展,芯片高度集成化帶來產品性能提升的同時,器件熱流密度不斷增大[1-2];另一方面由于芯片上電路或組件布置和集成不均,造成芯片上功率分布不均,不同位置的產熱量存在較大差異,引起局部熱流過高的熱點問題[3-4]。L. S. Maganti等[5]發現Intel?CoreTMi7芯片的表面溫差最大可達40 ℃,Hao Xiaohong等[6]指出熱點的熱流可超過整個芯片平均熱流一倍以上,R. Mahajan等[7]發現功率分布不均引起的熱流非均勻分布現象普遍存在。熱量不能及時散除會導致溫度過高,使電子設備可靠性受到嚴重影響。因此,關于電子器件中非均勻熱流的有效冷卻研究至關重要。

為有效解決“熱障”問題,需采取高效的散熱方式,如微/小通道散熱方法[8]。微通道散熱器具有高效緊湊且傳熱性能好的優點,已被廣泛應用于電子電路[9]、空調系統[10]、太陽能裝置[11]等。為進一步提高微/小通道的散熱能力,有研究者從尋找新型換熱工質角度出發,發現將相變微膠囊(microencapsulated phase change material, MEPCM)與傳統單相流體(基液)結合形成的相變微膠囊懸浮液(microencapsulated phase change material slurry, MEPCMS)在傳熱過程中能同時利用相變潛熱和對流傳熱發揮作用,進而提高傳熱效率。結合微/小通道和MEPCMS的特點,將MEPCMS作為微/小通道散熱工質的相關研究可參見文獻[12]。

關于MEPCMS早期的研究主要圍繞材料選擇[13]和制備方法[14]展開,逐漸開始實驗和模擬研究MEPCMS在管道內的流動傳熱特性。實驗方面,劉東等[15]通過實驗研究針肋式微矩形管程中相變微膠囊懸浮液的換熱特性,結果發現,相同雷諾數Re時,相變微膠囊懸浮液的溫度和管程壁面溫度均低于純水,且隨著Re的增大,降低程度增大。Zhang Guanhua等[16]對懸浮液在恒熱流圓管內流動時的傳熱特性進行實驗研究。結果表明,與水相比,相變材料在固態、固液混合態和液態3種狀態下時,懸浮液的努塞爾數Nu具有不同程度的提高。上述實驗研究均表明MEPCMS的傳熱性能優于基液[17]。隨著計算機技術及計算流體力學的發展,仿真模擬技術發揮著重要作用。以數值傳熱學及計算流體力學為基礎,仿真模擬技術對實際的物理現象以及問題能夠進行精確的求解計算。模擬方面,B. R. Far等[18]采用單相模型模擬懸浮液(芯材為正十八烷,基液為水)在帶有斜鰭的微通道內的層流換熱特性,發現散熱器的冷卻性能有所提高,但歐拉數也同時增加。高冬雪等[19]采用等效比熱模型研究恒熱流圓管內MEPCMS層流對流傳熱,結果發現,懸浮液在管內流動時存在相變區,流體溫度和壁面溫度隨Re的增加而降低,平均Nu隨Re的增大而增大。此外,對MEPCMS在微/小通道內的模擬研究熱度日益上升,但大多數研究[20-21]還停留在單相流模型,忽略了顆粒與基液間的相互作用,認為顆粒在基液中均勻分布,將懸浮液等效化為均質流體;計算精度一定程度上受限于唯象的物性模型,無法揭示顆粒對傳熱特性的影響機制。

離散相模型(discrete phase model, DPM)將基液和顆粒視為相互作用的連續相和離散相,是模擬MEPCMS流動傳熱的有效兩相流模型,A. B. S. Alquaity等[22]使用DPM模型模擬相變顆粒在微通道中的流動傳熱,發現相比單相流模型,DPM模型能更好地預測傳熱過程。目前使用DPM模型對MEPCMS流動傳熱的研究較少,而且未充分考慮顆粒與基液間的相互作用。本文采用DPM模型模擬矩形流道內MEPCMS的傳熱特性,基于兩相的方法更真實地揭示其流動傳熱過程,并重點研究了非均勻熱流密度、顆粒質量分數、進口流速對MEPCMS傳熱特性的影響。

1 模型及數值方法

1.1 模型建立及邊界條件

圖1所示為三維細小矩形通道的模型示意圖,矩形通道的長度L為100 mm、高度H為1 mm、寬度W為0.8 mm,經計算水力直徑Dh約為0.89 mm。u為進口流速(m/s),Tin為進口溫度(K)。

圖1 模型示意圖

本文根據矩形流道內MEPCMS的流動傳熱特點及相關參數范圍,進行如下假設:1)基液為水,流動傳熱過程按穩態不可壓縮處理;2)通道尺寸較小且流速較低,Re較小,流動為層流狀態;3)顆粒質量分數較低,不考慮顆粒與顆粒之間的相互作用;4)顆粒尺寸較小,且密度與水差異較小,不考慮重力和浮升力;5)顆粒視為均質球形,在流動傳熱中按質點處理,不考慮顆粒本身存在的溫度梯度;6)顆粒的相變通過比熱容的變化來體現。

邊界條件如圖1所示,左側為速度進口,顆粒與基液(即懸浮液整體)以相同的速度、溫度均勻注入,右側為壓力出口,周圍壁面為絕熱壁面,底部受熱面沿長度x方向分為3段:起始段La為40 mm、中間段Lb為20 mm、末尾段Lc為40 mm,各段所受恒定熱流密度分別為qa、qb、qc。MEPCM在流動過程中吸收底部的熱量發生相變,最后流體從出口流出。

1.2 控制方程

本文基于歐拉-拉格朗日方法的DPM模型進行數值模擬,基液和顆粒分別被視為連續相和離散相,在Euler坐標系下求解連續相控制方程,在Lagrangian坐標系下對顆粒運動微分方程進行積分來追蹤顆粒軌跡,在各自控制方程中添加相互作用源項來實現相間耦合。連續相的控制方程如下:

連續方程:

(1)

動量方程:

(2)

F=mpfD(up-uf)+Fother

(3)

能量方程:

(4)

式中:u、ρ、μ、p、k、T分別為速度(m/s)、密度(kg/m3)、動力粘度(Pa·s)、壓力(Pa)、導熱系數[W/(m·K)]、溫度(K);F為兩相交換動量源項,N;E為單位質量工質具有的總能量,J/kg;S為兩相交換熱量源項,W;m為質量,kg;fD為顆粒弛豫時間的倒數,1/s;下標f為基液,p為顆粒。

顆粒運動方程如下:

(5)

(6)

Rerel=ρfdp|up-uf|/μf

(7)

(8)

式中:d為直徑,m;mpfD(uf-up)為基液施加在顆粒的曳力項;t為顆粒積分時間步長,s;Rerel為相對雷諾數;CD為曳力系數;對于光滑球形顆粒,a1、a2、a3為文獻[23]中提出的常數值;Fother為顆粒受到的附加力項。本文考慮了4種附加作用力:

Fother=FV+FP+FT+FSL

(9)

(10)

(11)

(12)

(13)

式中:FV為虛擬質量力,N,是離散相相對于連續相作加速運動而產生的附加力;FP為壓力梯度力,N,是流體中的壓力梯度作用在顆粒上的力;FT為熱泳力,N,是由于溫度梯度存在導致顆粒受到液相分子熱運動的作用力不平衡所致[24];DT,p為熱泳系數;FSL為薩夫曼升力,N,是由顆粒周圍流體差異形成的剪切梯度升力;K為常系數,取值為2.594;v為運動粘度,m2/s;dij、dlk、dkl均為流體應變張量,1/s。

考慮顆粒與基液之間的對流換熱,顆粒的能量方程如下:

(14)

(15)

式中:hp為顆粒表面的流體-顆粒傳熱系數,W/(m2·K),根據文獻[25]中顆粒的努塞爾數Nup表達式獲得;cp為比定壓熱容,J/(kg·K);A為表面積,m2;Prf為基液的普朗特數。

能量方程源項即兩相交換熱量源項,通過確定顆粒經過每一個控制網格單元時熱能的變化進行計算:

(16)

1.3 相變微膠囊懸浮液的物性

本文基液為水,顆粒采用F. Dammel等[26]制備的微膠囊顆粒,芯材為正二十烷(n-Eicosane),殼材為聚甲基丙烯酸甲酯(polymethylmetharcylat, PMMA),給出的顆粒及懸浮液的物性參數如表1所示。本文中水的物性參數隨溫度變化而變化[27],除比熱容外顆粒的其它物性參數均不隨溫度變化。

表1 材料物性參數

單個顆粒直徑為10 μm,微膠囊顆粒中芯材的質量占比為85%,顆粒的物性公式如下[28-30]:

y=mc/ms

(17)

(18)

(19)

(20)

式中:y為核殼質量比;下標c、s分別表示芯材、殼材,芯材的密度為固態和液態時的平均密度,顆粒的導熱系數根據式(20)計算;使用等效比熱容法來處理顆粒的相變,即將顆粒的比定壓熱容視為溫度的函數[28-30]:

cp,p=

(21)

式中:hsf為顆粒的相變焓值,kJ/kg,取值為247.3;T1、T2、Tmr分別為相變起始溫度、相變終止溫度、融化溫度,K,取值分別為 309.55、313.05、3.5。

1.4 求解策略

采用基于有限體積法的Fluent 2020R2軟件,對控制方程進行離散化并迭代求解,分別使用SIMPLE算法計算壓力-速度耦合,Second Order Upwind離散動量方程和能量方程,設置控制方程(連續性方程、動量方程、能量方程)的殘差收斂標準分別為10-5、10-7、10-9。

1.5 網格無關性驗證及模型驗證

為確保數值模擬結果可靠,進行網格無關性檢驗。將通道底部加熱面的3段熱流密度均設置為5 W/cm2,進口流速為0.4 m/s,進口溫度為298.15 K,劃分4種網格如表2所示。由表2可知,640 000的總網格數約為360 000的1.8倍,但兩者流體進出口溫差相差僅為2.4%;綜合考慮計算成本和精確度,選用總網格數為360 000的網格進行模擬計算。

表2 網格無關性檢驗

由于MECPMS的強化對流換熱機理相似,為保證數值方法的合理性,本文利用不同形狀通道內MEPCMS的對流換熱實驗對DPM模型進行驗證。文獻[31-32]分別研究了MEPCMS在圓管和矩形通道中的對流傳熱特性,其中Zeng Ruolang等[31]的通道尺寸和工況條件為:圓管長度為1.46 m、內徑為4 mm,顆粒粒徑為8.2 μm,Re為1 283(層流),熱流密度為2.495 8 W/cm2,進口溫度為283.75 K。Rao Yu等[32]的通道尺寸和工況條件為:矩形通道長為150 mm、寬為2 mm、高為4.2 mm,顆粒粒徑為4.97 μm,質量流量為0.05 kg/min(層流),熱流密度為1.923 W/cm2,進口溫度為297.15 K。采用本文的數值模型和求解方法模擬文獻[31-32]的研究,結果如圖2所示。由圖2(a)可知,DPM模型對MEPCMS局部努塞爾數Nux的預測與實驗數據[31]吻合較好,平均相對誤差小于13%;且DPM模型對兩種質量分數cm的壁面溫度分布趨勢預測也與實驗數據[32]吻合較好,最大相對誤差為2.9%(如圖2(b))。因此,本文采用DPM方法得到的模擬結果具有可信度。

圖2 模型驗證

2 結果與分析

為更好地反映對流傳熱特性,引入局部努塞爾數Nux、平均壁溫Tw,av(K)、平均傳熱系數hav[W/(m2·K)]和進出口壓降Δp(Pa),計算如下:

(22)

(23)

(24)

(25)

Δp=pin-pout

(26)

式中:q、Dh、Tw、Tb分別為熱流密度(W/cm2)、通道水力直徑(mm)、壁面溫度(K)、流體質量加權平均溫度(K)。

2.1 相變微膠囊顆粒軌跡分析

為便于分析,從受熱底部(-y方向)視圖給出u=0.40 m/s、cm=10%時,不同熱流分布下顆粒的比熱容和軌跡沿流動方向(+x方向)的變化,如圖3所示。

圖3 懸浮液中顆粒的比熱容和軌跡 (-y視圖)

由圖3可知,由于壁面周圍流速低主流區域流速高,顆粒起始相變區域形狀(x-z平面)呈拋物線型。隨著熱點區域向前移動和局部熱流密度增大,底部顆粒可更快達到相變起始溫度,且完成相變所需距離變短。

為分析熱泳對顆粒行為及換熱的影響,圖4從+z視圖給出u=0.40 m/s、cm=10%、qa-qb-qc=9-5-5 W/cm2時顆粒比熱容和軌跡變化。由圖4可知,考慮熱泳效應后底部的顆粒有向上偏移現象,這是由于顆粒受到與溫度梯度方向相反的熱泳力,使其從高溫區域向低溫區域運動,致使顆粒相變起點推遲、長度延長。表明顆粒與水的相互作用一定程度上會對懸浮液的換熱過程有影響。此外,大部分顆粒由于流速過快導致還未開始相變便流出通道。

圖4 懸浮液中顆粒的比熱容和軌跡(+z視圖)

2.2 非均勻熱流密度的影響

在u=0.40 m/s、Tin=298.15 K、cm=10%時,采用非均勻熱流進行加熱,記熱流密度分布為qa-qb-qc,本文研究中背景熱流為5 W/cm2,局部熱流為7 W/cm2和9 W/cm2,局部熱流所處位置可位于起始段(La)、中間段(Lb)、末尾段(Lc)中的一段。

圖5、圖6分別為u=0.40 m/s、cm=10%時,通道內不同區域的平均壁溫Tw,av和平均傳熱系數hav隨熱流分布的變化。由圖5可知,局部熱流所處位置會影響該區域及其后面區域的Tw,av,而不影響其前面區域的Tw,av,這是由于流體流過局部熱點區域后吸收熱量溫度上升,對沿程壁面的冷卻效果有所影響;且隨著局部熱流密度增大,Tw,av幾乎呈線性增加。由圖6可知,隨著熱點區域熱流密度的增加,hav在熱點區域的上游區域保持不變,在熱點區域隨之增大,而在熱點區域的下游區域逐漸減小。這是由于加大熱流密度會加速懸浮液中顆粒相變,換熱效果得到提高,但Tw比Tb增加明顯(如圖7(a)),一定程度上削弱了換熱效果。與均勻熱流條件(qa=qb=qc=5 W/cm2)相比,非均勻熱流會影響溫升情況,改變區域間的平均壁面溫差,使高溫區域范圍發生變化;同時會改變區域平均傳熱系數的相對大小。當qa-qb-qc=9-5-5 W/cm2時,相比均勻熱流條件,La、Lb、Lc各區域的Tw,av分別提高2.31、1.77、1.27 K,hav分別提高1 091 W/(m2·K)、降低2 339 W/(m2·K)、降低1 300 W/(m2·K)。

圖5 非均勻熱流條件下Tw,av曲線

圖6 非均勻熱流條件下hav曲線

圖7 非均勻熱流條件下Tw、Tb曲線和Nux曲線

2.3 顆粒質量分數的影響

采用底部熱流密度分布為qa-qb-qc=9-5-5 W/cm2,u=0.40 m/s,選取cm為5%、8%、10%進行模擬,模擬結果如圖8所示。

由圖8(a)可知,相同x處MEPCMS的溫度比水低,且cm越大溫度越低,這是由于懸浮液中相變微膠囊顆粒受熱發生相變,吸收熱量減緩溫度上升,cm越大,表明顆粒數目越多,相變時吸收熱量越多,溫升進一步得到縮小。由圖還可知,在起始段水和MEPCMS溫度差異較小,這是由于此時相變微膠囊還未開始發生相變,隨著顆粒吸收了足夠多的熱量,達到相變溫度,顆粒開始發生相變,小部分熱量以顯熱被吸收,大部分熱量被顆粒以相變潛熱吸收,因此溫度上升緩慢。與純基液相比,在cm=10%時,懸浮液最高能夠減小壁面溫升8.79%和流體溫升15.14%。

由圖8(b)可知,水和MEPCMS在入口段的Nux較大,然后沿流動方向下降,這是由于入口段邊界層逐步變厚所致;到中間段時,Nux會出現大幅下降,再略有上升后保持基本穩定,這是由于熱流密度變化所致。距離入口一小段區域后,懸浮液的Nux均高于水,且隨著顆粒濃度增大,Nux比水高出的幅度越明顯。這是由于微膠囊顆粒相變強化了傳熱過程。與純基液相比,在cm=10%時,懸浮液Nux最大提升幅度為8.26%,Lc段平均提升幅度為10.16%。

2.4 相變微膠囊懸浮液進口流速的影響

采用底部熱流密度分布為qa-qb-qc=9-5-5 W/cm2,cm=10%,選取u為0.40、0.45、0.50、0.55 m/s進行模擬,模擬結果如圖9所示。

圖9 不同進口流速時Tw、Tb曲線和Nux曲線

由圖9可知,隨著流速的增加單位時間通過微通道的流量隨之增加,單位時間內可帶走的熱量更多,所以無論是純基液還是懸浮液,整體溫升均減小,Nux均增大。當流體流經熱流密度較小的中間段時,Tw減小和Tb曲線斜率降低。由于顆粒融化吸熱效應,相同速度的MEPCMS要比水對通道壁面的冷卻效果更好;但可以看出,隨著速度的增加,相同速度時兩種流體的Tw曲線之間的距離趨于縮短,當流速遞增時,對Tw降低的幅度依次為8.79%、8.37%、7.96%、7.21%,即MEPCMS增強換熱的效果有所降低。這是因為懸浮液的流速過快,大部分顆粒還未相變就流出流道,不能充分發揮其相變吸熱能力,更多的熱量被流體工質以顯熱形式吸收。

熱流分布為qa-qb-qc=9-5-5 W/cm2時,u和cm對進出口壓降Δp、hav的影響如圖10所示。由圖10可知,相同顆粒質量分數時,隨著進口流速提高,流體流量增加可帶走更多熱量,hav增大,同時沿程流動阻力增加,Δp增大;另一方面,相同進口流速時,隨著顆粒質量分數增大,更多顆粒發揮相變潛熱優勢吸收更多熱量,hav增大,同時顆粒與基液的相互作用增強,Δp增大。

圖10 不同進口流速和質量分數時Δp和hav曲線

3 結論

本文采用DPM模型對細小三維矩形通道內MEPCMS的傳熱特性進行研究,并重點剖析熱流分布、質量分數和進口流速的影響,得到如下結論:

1) 相變微膠囊顆粒相變行為可強化懸浮液在通道內的對流換熱,在熱流密度分布為9-5-5 W/cm2、進口流速為0.40 m/s、質量分數為10%時,懸浮液最高能夠減小壁面溫升8.79%和流體溫升15.14%;懸浮液Nux最大提升幅度為8.26%,出口段平均提升幅度為10.16%。

2) 在研究范圍內,當局部熱流位置越向前時,能更早發揮顆粒相變吸熱能力,且隨熱流密度增大顆粒越早完成相變過程;同時顆粒的軌跡會受熱泳力影響向上移動偏離高溫底部。

3) 從高熱流區域過渡到低熱流區域,局部努塞爾數先減小再增大,最后基本保持穩定,且相同位置處懸浮液的局部努塞爾數更大。

4) 增加顆粒質量分數和懸浮液的進口流速均能提高平均傳熱系數;但同時加大了沿程流動阻力和兩相相互作用的影響,導致進出口壓降增大。

本文受廣東省自然科學基金(2019A1515012119)資助。(The project was supported by the Natural Science Foundation of Guangdong Province (No. 2019A1515012119).)