常 建，蔡杰進，譚 冰

(華南理工大學 電力學院，廣東 廣州 510640)

空間熱管堆采用熱管作為堆芯向外界熱傳導的媒介，由于缺少重力的作用，毛細芯對冷卻劑的泵送就更重要[1-2]。為更有效地進行傳熱，對毛細芯的相變研究越來越有必要。由于毛細壓力的限制，在毛細芯中所能提供的最大供液流量限制了熱管所能承受的最大熱輸入，過高則會導致毛細芯的燒干。為改善熱管性能，毛細結構需優(yōu)化其泵送能力和滲透率[3-4]。在大多數(shù)情況下，毛細芯結構被認為是一個連續(xù)的多孔介質(zhì)，但越來越多的研究發(fā)現(xiàn)薄膜蒸發(fā)和局部孔隙的馬拉格尼對流等微觀尺度現(xiàn)象的影響在相變過程中起重要作用。一般來說，薄膜蒸發(fā)發(fā)生在幾個微米長的區(qū)域附近的固-液-汽交界處，長期以來一直被認為是在這樣的系統(tǒng)中傳熱的主要模式。薄膜傳熱的高效率是由大的分離壓力梯度以及極低的熱阻共同影響的[5-6]。三相接觸點附近的強烈蒸發(fā)導致了沿半月板界面方向的溫度梯度，進而導致表面張力梯度，從而產(chǎn)生熱毛細對流[7-8]。

本文通過對單個毛細芯孔隙內(nèi)的界面蒸發(fā)過程進行模擬，以對整個毛細芯的相變提供更準確的模型。首先在相場法的基礎上，建立一個從毛細彎液面蒸發(fā)到空氣中的三維模型，考慮界面處的馬拉格尼效應、熱浮力效應、界面處的蒸發(fā)相變過程以及蒸汽擴散產(chǎn)生的反作用力。然后將模擬結果與實驗進行對比驗證。最后對界面處的蒸發(fā)進行敏感性分析，討論壁面過熱度、接觸角和管徑對相變過程的影響。

1 薄膜蒸發(fā)傳熱傳質(zhì)模型

擴散界面法不僅是一個簡單的數(shù)學計算，而且包含著一定的物理意義。它是通過相場變量φ來得到界面層的信息，而不是去直接追蹤兩種流體界面的變化。它把表面張力等效為場變量的梯度與化學勢的乘積，并將它作為一個體積力加入到Navier-Stokes方程中。相場法不但可計算流體界面的對流，而且還保證系統(tǒng)總能量合理減少[9]。

擴散界面法認為界面是一個薄的、物理性質(zhì)連續(xù)變化的區(qū)域，且兩相的物性參數(shù)也在界面連續(xù)變化。本文用相場變量φ區(qū)分兩相，類似于水平集方法。計算時，第1相的值取-1，第2相的值取1，界面處的值從-1到1變化[10]。計算區(qū)域內(nèi)每相所占的份額V可寫成相場變量的函數(shù)：

(1)

式中，下標v、l分別表示氣體和液體。

在模型計算中涉及的物性參數(shù)，包括導熱系數(shù)、黏度、密度以及定壓熱容等，可定義為：

B=BlVf,l+BvVf,v

(2)

式中，B為相界面對應的物性參數(shù)。

λ為混合能量密度，其與表面張力和擴散界面厚度關系如下所示：

(3)

式中，εpf為界面厚度參數(shù)，取決于界面處網(wǎng)格的大小。為提高界面變化的準確度，εpf應盡可能小，以保證相場變量φ及兩相物理性質(zhì)在相界面的均勻變化，但同時為保證數(shù)值計算的收斂性，相界面的厚度也要與在相界面處的網(wǎng)格密度相匹配。通常取界面處最大網(wǎng)格的1/2作為εpf的值，具有較好的模擬結果。

(4)

式中：ρ為密度；n為界面法向量。對n向量求偏導，并整理可得：

(5)

(6)

式中：u為速度場；ξ為相界面平滑函數(shù)：

(7)

在核態(tài)沸騰的過程中，通過相界面所發(fā)生的傳熱傳質(zhì)過程是最為重要和特殊的過程，在相變過程中，熱流密度和汽化潛熱對相變過程的影響十分重要。相變過程的質(zhì)量轉換率可寫為：

(8)

式中：q為熱流密度；ΔHvl為汽化潛熱；Ml為液體摩爾質(zhì)量。在相界面處，考慮熱流僅通過導熱方式傳導。將熱流密度通過傅里葉導熱定律展開，式(8)變?yōu)椋?/p>

(9)

式中：T為溫度；k為導熱系數(shù)。上述推導中，主要考慮了流體蒸發(fā)造成的界面移動，由于有出口存在，蒸汽的產(chǎn)生未對界面的移動產(chǎn)生直接影響，但蒸發(fā)時相變帶來的密度變化大，在微小的空間內(nèi)，汽體從出口擴散時對界面的反作用力不容忽視。本文考慮由于相變導致的氣液移動速度，汽體對界面的反作用力fvl可表示為：

(10)

在確定了相變質(zhì)量通量后，可對其他控制方程同樣通過添加源項和界面平滑參數(shù)的方式，修正在相界面處的表達式。

相場控制方程(C-H方程)：

(11)

(12)

式中：t為時間；γ為遷移率。由于各方程所引入的參數(shù)較多，有必要對各參數(shù)進行更詳細的解釋，以幫助理解方程意義。在對氣泡運動的研究中，取常溫下水的表面張力σ為0.076 N/m。

(13)

在微尺度傳熱下，這些溫差變化會大到足以對流場造成實質(zhì)性的影響。此外，由于流體會傳遞熱量，因此溫度場也會受到流場變化的影響。

動量方程(N-S方程)：

(14)

其中：

(15)

式中：Fg為重力項；I為應力張量。

能量方程：

(16)

式中，cp為比定壓熱容。

2 計算模型與網(wǎng)格驗證

通道內(nèi)最小的網(wǎng)格單元邊長為10 μm，其長度遠大于滑移長度，計算中視為不可滑移邊界條件。首先，壁面處沒有熱通量，流體在毛細力的作用下上升，彎液面在重力的作用下重構，并達到穩(wěn)態(tài)，然后添加壁面初始溫度，液體在彎液面通過蒸發(fā)被轉化為蒸汽并通過管中擴散，同時下方的入口處由于毛細作用對流體進行補充，彎液面隨著蒸發(fā)產(chǎn)生波動。表1列出了計算中所需的物性參數(shù)。

表1 數(shù)值模擬計算采用物性參數(shù)Table 1 Physical property parameter used in numerical simulation calculation

圖1為微通道彎液面蒸發(fā)計算模型及自適應網(wǎng)格劃分。數(shù)值模型為采用軸對稱的簡化的二維通道(圖1a)，藍色部分為液體，灰色部分為氣體，初始彎液面構造為圓弧，在相界面初始化過程中，會根據(jù)重力和接觸角對界面進行重構，入口為壓力入口，其值為液面產(chǎn)生的壓力，出口壓力為0。本研究的重點是彎液面處的蒸發(fā)速率，由于彎液面會隨著蒸發(fā)不斷下移，故本文在模擬該模型時，力圖在界面梯度較大的區(qū)域中創(chuàng)建密集網(wǎng)格。但有時很難預測梯度驟變出現(xiàn)的位置。還有一些瞬態(tài)情況，驟變的梯度會移動。這可通過在出現(xiàn)梯度驟變的所有區(qū)域創(chuàng)建細化網(wǎng)格來解決，但通常計算成本很高，較好的解決方法是采用自適應網(wǎng)格劃分。基礎粗化網(wǎng)格用于在一定時間間隔內(nèi)推進求解，隨后使用求解結果基于某個指示函數(shù)來細化網(wǎng)格。最后，使用自適應網(wǎng)格再次模擬時間間隔。

本文采用了三角形網(wǎng)格進行區(qū)域劃分。由于這是一個兩相模型，且主要考慮界面處的蒸發(fā)，因此自適應基于進行相場變量的梯度。圖1d中黑色線為初始界面，紅色線為移動后的界面，顯示了自適應時間間隔結束時的解。由圖1可看出，隨著界面的移動，對界面區(qū)域，即相場變量的梯度變化大的范圍進行了網(wǎng)格細化，達到了所需的計算精度。

本文使用水作為工作液進行了網(wǎng)格無關性研究，如表2所列，生成了4組網(wǎng)格，其對計算結果的影響列于表3。經(jīng)過對比發(fā)現(xiàn)，網(wǎng)格2相比網(wǎng)格1網(wǎng)格數(shù)量減少了1/2，計算結果的誤差在1%左右，滿足網(wǎng)格獨立性的驗證，對下文的計算中均采用網(wǎng)格3作為網(wǎng)格劃分的密度。

a——微通道彎液面蒸發(fā)計算模型示意圖；b——微通道兩相及彎液面處的基礎網(wǎng)格； c——根據(jù)相界面的移動和梯度變化后經(jīng)過細化的界面附近網(wǎng)格；d——根據(jù)加密后的網(wǎng)格計算得到的移動過后兩相界面圖1 微通道彎液面蒸發(fā)計算模型及自適應網(wǎng)格劃分Fig.1 Microchannel meniscus evaporation calculation model and adaptive meshing

表2 網(wǎng)格獨立研究中使用的不同網(wǎng)格Table 2 Different grids used in grid independent research

表3 網(wǎng)格獨立驗證結果Table 3 Grid independent verification result

3 計算結果驗證與討論

Buffone等[11]通過粒子成像系統(tǒng)拍攝不同管徑的石英毛細孔隙內(nèi)的流場和溫度場分布，實驗工質(zhì)為乙醇和丙酮，采用ITO鍍膜進行電加熱，得到了在彎液面附近看到環(huán)形渦旋流場。本文采用實驗中的相似的物性和條件，模擬了豎直放置的微通道彎液面蒸發(fā)，結果如圖2所示。模擬所得彎液面附近的流場與Buffone等實驗結果進行了對比，驗證了數(shù)值模型的正確性。在彎液面附近液面中心的流體沿著相界面向三相接觸線附近流動，推動壁面附近的流體回流，最終在入口補償流的作用下，發(fā)生偏轉，補充位于彎液面中心的流體。這種現(xiàn)象稱之為馬拉格尼流動。液-汽界面處的馬拉格尼對流驅動液體從管中心線的較熱區(qū)域流向壁的較冷區(qū)域，加速了熱對流，促進了蒸發(fā)。入口處的液體流入速度約為1 mm/s，渦旋附近的流速遠大于入口的流速，因此界面附近的流動主要是由溫度梯度導致的馬拉格尼流引起的。剛流入時，流體流線平行于管軸，直到接近彎液面在渦旋的影響下發(fā)生偏轉，最后從彎液面的中心流向接觸線附近。從Buffone等實驗結果也可清晰看到兩個渦旋產(chǎn)生在界面附近，同時補給的流體明顯發(fā)生了偏轉，與模擬所得的流場圖幾乎一致。該結果驗證了模型的正確性，表明建立的基于相場法的界面蒸發(fā)模型可用于蒸發(fā)速率的計算。

圖2 豎直放置的微通道中流場圖Fig.2 Flow field diagram in vertically placed microchannel

不同壁面過熱度下彎液面蒸發(fā)時的溫度場如圖3所示?？煽闯龀跏紩r壁面附近溫度分層分布，因為蒸發(fā)的作用，界面處于飽和溫度，相對于過熱的液體，界面溫度最低，且從界面向下溫度逐漸增高。隨著液體在馬拉格尼流渦旋作用下，界面處的高溫流體被帶入管道的中心，最終在彎液面附近產(chǎn)生較大的溫度梯度，且越靠近壁面，溫度梯度越大，說明壁面附近的蒸發(fā)越強烈[12-14]。對比圖2中的流場圖可看出，溫度場的分布和流場的分布具有一致性。不同于毛細孔隙在層流下的壁面分層傳熱，馬拉格尼渦旋的作用下，壁面處的溫度可更快地傳遞到蒸發(fā)界面，增強了換熱。對于不同的過熱度，具有相似的溫度場分布，過熱度的影響將在接下來的結果討論中分析。

圖3 不同壁面過熱度下彎液面蒸發(fā)時的溫度場Fig.3 Temperature field of meniscus evaporation under different wall superheat degrees

圖4為不同尺寸的微通道水平放置時在重力作用下渦旋發(fā)生偏移的模擬結果與Hemanth等的實驗結果對比。在模擬中，毛細孔隙的軸線是水平的，因此重力效應會影響彎液面附近的流動方式[15]。圖中觀察到的不對稱流動是浮力和馬拉格尼效應的綜合結果。如圖中實驗結果所示，圖中每個部分均疊加有多個圖像，以提高對比度。通過沿垂直中心平面獲得的長曝光顆粒條紋圖像表征了浮力和熱毛細作用對流動的相對影響，渦旋的偏移與模擬結果一致。圖中的箭頭指向彎液面上停滯點的位置，模擬中該位置處的界面由于補償流的作用溫度相對最低，熱毛細作用最小，與兩側所產(chǎn)生的切應力處于平衡狀態(tài)，當在界面處發(fā)生蒸發(fā)時，它使彎液面區(qū)域液體較離開彎液面的回流液體相對溫度低。溫度的這種差異同時會產(chǎn)生浮力效應，其強度隨管徑的增加而增加。對于研究的最小直徑為75 μm的毛細孔隙，浮力的影響可忽略不計，并且觀察到與豎直的毛細孔隙相同的對稱流型(圖4a)。對于200 μm的毛細孔隙，渦旋的不對稱性開始產(chǎn)生，發(fā)現(xiàn)逆時針底部渦旋大于順時針上部渦旋，這表明存在不對稱的環(huán)形渦旋，前面提到的停滯點不再位于液面的中心，開始向上偏移。浮力效應隨著毛細孔隙直徑的增加而增大，如圖4c、d所示，逆時針渦旋的大小隨著毛細孔隙直徑的增加而變大，停滯點沿彎液面向上移動，浮力增加?？煽闯?，該流動被分為兩個區(qū)域：由表面張力引起運動的上部區(qū)域和由密度差驅動的下部區(qū)域。同時在模擬中還可看到，直徑的增大導致重力在界面的影響也越明顯，圖4a、b中，上下兩側的三相接觸點大致在同一垂直平面上，圖4c、d中，則由于重力的作用，下方界面向前伸展，在加熱時，這種不平衡增大了下方的蒸發(fā)界面，加劇了渦旋的偏移。因此，浮力驅動渦旋的尺寸隨毛細直徑的增加而增大。在界面處產(chǎn)生的相對溫差驅動了浮力引起的流動，其大小取決于蒸發(fā)的質(zhì)量通量。

圖5示出了不同直徑毛細孔隙彎液面處的蒸發(fā)速率。x0為管徑，x為界面與對稱軸的軸向距離。從圖中可看出，對稱軸附近的蒸發(fā)速率均較小，且在約60%的長度區(qū)域內(nèi)，蒸發(fā)速率均變化不大，在靠近壁面的區(qū)域內(nèi)，蒸發(fā)速率快速增長，在接近壁面處，不同管徑的界面具有相同的蒸發(fā)速率。在對稱軸附近，直徑較小的彎液面蒸發(fā)速率較高，其值最多僅有壁面附近蒸發(fā)率的1/10，這是由于較小的彎液面從壁面處由渦旋帶來的熱量更多，其溫度梯度更大，因而小管徑的彎液面局部蒸發(fā)速率更大。從圖5還可看出，蒸發(fā)主要集中在壁面附近，壁面附近的蒸發(fā)速率與管徑尺寸無關。

圖5 不同直徑毛細孔隙界面處的蒸發(fā)速率Fig.5 Evaporation rate at interface of different diameter capillarie interstices

圖6示出管徑為0.4 mm時不同過熱度ΔT和接觸角下毛細孔隙界面處的蒸發(fā)速率。從圖6a可看出，蒸發(fā)速率是過熱度的函數(shù)，呈正比關系，過熱度的增加并未改變彎液面的相對蒸發(fā)速率，即過熱度對渦旋的位置未產(chǎn)生影響，壁面處的蒸發(fā)速率遠大于彎液面中心停滯點。圖3中的溫度場也同樣證明了這一結論，壁面處的界面溫度幾乎恒定，且溫度梯度變化大，蒸發(fā)劇烈；中心處由于補償流的作用，溫度梯度小，蒸發(fā)慢。

微通道尺寸：a——75 μm；b——200 μm； c——400 μm；d——762 μm圖4 微通道水平放置時在重力作用下 渦旋發(fā)生偏移的結果Fig.4 Vortex shift results of microchannels placed horizontally under action of gravity

壁面附近的蒸發(fā)速度是壁面上液體厚度的函數(shù)，而這個厚度是由界面的接觸角決定的，為接觸角的函數(shù)。從圖6b可看出壁面處液膜厚度與蒸發(fā)速率的關系，其值相對于接觸角的增加，變化不大。在接近壁面處，接觸角減小導致更好的潤濕性，液體更有利于向壁面流動，補償壁面處的蒸發(fā)，且液膜的厚度更小，蒸發(fā)相對更劇烈，考慮到軸向液膜的長度，則接觸角小的液膜蒸發(fā)速率遠大于大接觸角的情況。本文潤濕性在蒸發(fā)過程中更主要的是提供更薄的液膜，在厚度一致的情況下，蒸發(fā)的速率基本相同，接觸角小意味著有更長的薄膜蒸發(fā)區(qū)域，因此在整個蒸發(fā)表面，有更大的蒸發(fā)速率。

圖6 管徑為0.4 mm時不同過熱度和接觸角下毛細孔隙界面處的蒸發(fā)速率Fig.6 Evaporation rates at capillary interstice interface under different superheat degrees and contact angles with pipe diameter of 0.4 mm

4 半月板上蒸汽擴散通量的不均勻性

通過圖6可知界面處某點的蒸發(fā)速度是液面厚度和溫度梯度的函數(shù)，對于整個界面的蒸發(fā)速率，需對界面作積分以體現(xiàn)在軸向上的液面分布。毛細孔隙界面處蒸發(fā)模型示意圖如圖7所示，d為半徑，r為界面上點到對稱軸距離，r0為彎液面的等效半徑。這里假設重力對液面的影響簡化為圓弧，設彎液面最低點的高度為0，則彎液面的某點處的位置表示為：

圖7 毛細孔隙界面處蒸發(fā)模型示意圖Fig.7 Diagram of evaporation model at interface of capillary interstice

r0=d/cosθ

(17)

(18)

實際中，熱管毛細芯的流體均是選擇潤濕性較好的流體，以達到流體的泵送目的，其接觸角大多在60°以內(nèi)，對圖6b中60°接觸角進行膜厚與蒸發(fā)速率的曲線擬合得到：

m″=22.5+(2.1er/0.11+1.3×10-6er/0.023)/cosθ

(19)

式中，m″為蒸發(fā)速率，g/(m2·s)。

5 結論

本文通過對單個毛細芯微孔內(nèi)的毛細蒸發(fā)過程進行建模，考慮了馬拉格尼效應引起的界面剪切作用和液體區(qū)域的熱浮力效應，研究了微管內(nèi)部的液體在彎液面的傳熱和傳質(zhì)過程，具體結論如下。

1) 在彎液面和管壁的三相接觸處較在彎液面中心處的蒸發(fā)更強。從彎液面中心到壁的溫度梯度導致從管的中心朝向壁的馬拉格尼流。

2) 重力的作用導致液體中的渦旋不對稱，這種不對稱對小管徑影響不大，在大于0.4 mm的情況下更顯著。

3) 液面中心處的蒸發(fā)速率均較小，在靠近壁面的區(qū)域內(nèi)，蒸發(fā)速率快速增長，在接近壁面處，不同管徑的界面具有相同的蒸發(fā)速率。

4) 壁面附近的蒸發(fā)速度是壁面上液體厚度和過熱度的函數(shù)，而該厚度是由界面的接觸角決定的，為接觸角的函數(shù)。

5) 隨著孔隙在毛細芯中的液體有更多的比表面積，文中可看出接觸線附近蒸發(fā)速率最大，因此孔隙減小，蒸發(fā)速率增大。但較小的孔隙會導致較差的滲透率，發(fā)生燒干，需綜合考慮孔隙減小對蒸發(fā)的影響。