池沸騰孤立氣泡生長過程中微液層蒸發(fā)影響的實驗和模擬耦合分析

潘豐，王超杰，母立眾，賀纓

（大連理工大學能源與動力學院，遼寧大連116024）

引 言

沸騰傳熱是復雜的相變傳熱過程，幾十年來，針對核態(tài)沸騰中的高效換熱機制，學者們提出了諸多機理模型，包括瞬態(tài)導熱、微對流和液層蒸發(fā)機理等。其中，瞬態(tài)導熱理論認為熱量的傳遞主要發(fā)生在氣泡脫離階段，氣泡的脫離周期性地排開加熱表面上的過熱液體使得冷液體占據(jù)加熱表面并與之發(fā)生熱交換[1]。微對流理論認為沸騰換熱的本質(zhì)為被強化了的強制對流換熱過程，其高效的換熱性能來源于由氣泡的周期性脫離而帶來的強化效應，而并非相變過程本身[2]。而液層蒸發(fā)機理則認為，在氣泡生長初期，不斷生長的氣泡將周圍液體推開，在此過程中少量液體留在加熱壁面上形成很薄的一層微液層[3]。同時，在氣泡外側的熱邊界層內(nèi)則形成一層大液層[4]。氣泡生長所需要的能量主要來源于加熱表面上氣泡底部及附近的過熱液體層的相變潛熱[5-6]。Tange等[7]利用背光法和PIV測速技術測量了沸騰過程中近壁面的流場，據(jù)此對瞬態(tài)導熱、微對流和液層蒸發(fā)三種換熱機理的傳熱系數(shù)進行了評估。他們認為隨著沸騰熱流的增大，液層蒸發(fā)占沸騰換熱的比例持續(xù)增大，而微對流所占的比例逐漸減小。在高熱通量區(qū)域(0.65 MW/m2)，液層蒸發(fā)所占的比例達到90%以上。

圖1 氣泡底部微液層與大液層區(qū)域示意圖Fig.1 Schematic diagramof the microlayer and macrolayer underneath a vapor bubble

根據(jù)液層厚度的不同，氣泡生長過程中底部液體層可分為三個區(qū)域，由內(nèi)向外依次為：吸附層，微液層，大液層。其中吸附層為納米量級，其換熱通常忽略不計。如圖1所示，氣泡中心到氣泡底部氣液界面與加熱表面接觸點的距離定義為基圓半徑。存在于氣泡基圓半徑內(nèi)側的為微液層，尺度較小，為微米量級。而大液層則為基圓半徑以外到熱邊界層δt與氣液界面交界位置Rt的液體層，其厚度為幾百微米到毫米量級[8]。Sakashita等[9]測量了近臨界熱通量(CHF)處加熱表面附近的空泡率，實驗結果與大液層蒸發(fā)理論有較好的一致性：在高熱通量區(qū)域，大液層蒸發(fā)消耗的熱量占沸騰過程總傳熱量的94%。相比于大液層區(qū)域，盡管氣泡底部的微液層區(qū)域較小，但由于該液體層在厚度方向上熱阻較小，平均熱通量高達1 MW/m2，在低熱通量沸騰區(qū)域內(nèi)，其換熱量不容忽視[6]。H?hmann等[10]使用熱色液晶(TLC)對沸騰過程中加熱表面溫度進行測量，結果表明氣泡生長過程中大液層區(qū)域的溫度基本穩(wěn)定，而在寬度約50μm的微液層區(qū)域出現(xiàn)了明顯的溫降。類似地，Stephan等[11]通過紅外相機觀測了FC-84和FC-3284沸騰過程中不銹鋼表面上的瞬態(tài)溫度分布，同樣發(fā)現(xiàn)加熱表面上的低溫區(qū)域出現(xiàn)在三相接觸線附近的微小區(qū)域內(nèi)；進一步的計算表明這一部分的換熱量達到了總換熱量的40%～60%。Nakabeppu等[6]則認為，氣泡生長所需的能量約50%來源于微液層，其余則來自于氣泡周圍的大液層，且這一比例在不同的過熱度下基本保持不變。陳宏霞等[12]對光滑單晶硅表面上不同熱通量下單氣泡沸騰過程中的壁面溫度演變規(guī)律進行了對比分析，結果發(fā)現(xiàn)在高熱流沸騰過程中，氣泡底部形成了較大的干斑區(qū)域，導致汽化核心中心位置溫度持續(xù)較高，而低溫區(qū)域出現(xiàn)在大液層附近，因此此時相變換熱主要發(fā)生在氣泡兩側的大液層區(qū)域；而在低熱通量區(qū)間內(nèi)，氣泡底部存在一層均勻的微液層，加熱表面上對應地呈現(xiàn)出較大的低溫區(qū)域，同時壁面溫度分布表現(xiàn)為“中間低，兩頭高”。同時，他們在低熱通量下的沸騰模擬結果也表明微液層的蒸發(fā)在氣泡生長初期的貢獻量高達90%以上[13]。而隨著氣泡的生長，微液層所占比重逐漸減小，在氣泡脫離前微液層的蒸發(fā)占總換熱量的52%。因此，觀測加熱壁面上氣泡底層的液層動態(tài)分布有助于揭示沸騰傳熱機理，提出沸騰強化換熱新手段[14]。

由于吸附層和微液層尺度較小，常規(guī)的觀測手段只能觀測到大液層區(qū)域，而沸騰過程中的微液層區(qū)域則處于氣泡宏觀基圓半徑內(nèi)側的微觀區(qū)域。近年來，借助先進的實驗觀測手段和逐漸完善的數(shù)值模型，沸騰過程中氣泡底部的液層厚度變化及蒸發(fā)過程可以被較好地捕捉到。Nakabeppu等[6]借助微電子芯片(MEMS)傳感器記錄了氣泡底層的溫度場的動態(tài)變化，由此獲得了氣泡底部微液層區(qū)域的演變以及干斑區(qū)域的擴張和再潤濕過程。他們測得基圓半徑內(nèi)側瞬態(tài)熱通量峰值高達2.8 MW/m2，并由此得到微液層厚度約為3μm。相比之下，基圓半徑外側熱通量則不到0.4 MW/m2。Gong等[15]在ITO透明加熱材料上進行了沸騰實驗，使用非接觸式的共焦光學探針，測量得到微液層的厚度為5～50 μm，大液層的厚度約為50～500μm。Zou等[8]應用飛秒級光源技術，產(chǎn)生了能在過冷沸騰過程中保持若干小時的單個蒸氣泡，并對不同潤濕性加熱表面上的氣泡底層的微液層進行了測量。結果表明，親水表面上的氣泡底層干斑區(qū)域較大且液層較薄，而疏水表面上微液層區(qū)域較大且液層較厚。利用激光消光法[16]和激光干涉法[17-19]對乙醇/去離子水沸騰實驗時氣泡底部液層結構的動態(tài)測量表明，微液層厚度均在1～5μm左右。同時，基圓半徑在氣泡生長初期不斷外擴，氣泡底層的微液層隨之迅速增加。隨著氣泡的生長，微液層由于不斷蒸發(fā)而逐漸變薄，并在中心區(qū)域出現(xiàn)干斑。而后，在氣泡生長末期，干斑區(qū)域隨著氣泡基圓半徑的收縮而重新潤濕，最終氣泡脫離時干斑區(qū)域與微液層區(qū)域均變?yōu)?。這些詳實的觀測結果為深入了解沸騰換熱機理，精確分析沸騰傳熱過程奠定了理論基礎。

目前，代表性的微液層蒸發(fā)模型[20-22]基本基于實驗觀測和氣泡動力學分析。其中應用最為廣泛的是Stephan等[23]基于雷諾潤滑理論提出的靜態(tài)微液層蒸發(fā)模型。該模型通過求解一個四階非線性常微分方程來計算微液層厚度與相變熱流。在此基礎上，Jiang等[24]考慮了微液層的橫向流動并結合接觸線的移動和接觸角的變化提出了動態(tài)微液層蒸發(fā)模型。相比之下，Sato等[25]提出的靜態(tài)液層蒸發(fā)模型因具有參數(shù)少且簡便有效等優(yōu)點而在近年來被經(jīng)常使用[13,26]。其模型基于Utaka等[16]的線性初始液層厚度假設提出：

該模型中唯一需要確定的參數(shù)為液層的初始斜率Cslope，對于去離子水在石英板上的沸騰過程，Utaka等[16]測得的取值為4.46×10-3，陳宏霞等[13]和王燁等[26]在模擬中也采用了這一數(shù)值。Sato等[25]在模擬銅表面上的沸騰過程時通過多次迭代匹配的方法，最終確定的數(shù)值為2.45×10-2。事實上，該參數(shù)在物理本質(zhì)上表征著微液層前端微觀接觸角[18]的大小，因而與工質(zhì)物性和加熱表面的潤濕性相關。對于不同的工質(zhì)和加熱面材料，該斜率的取值需要更多數(shù)據(jù)予以確認。

由于氣泡底層結構復雜，且伴隨著液層的蒸發(fā)和移動，尤其在不銹鋼、銅等傳統(tǒng)非透明金屬表面上，微液層的厚度更加難以測量。而相比之下，氣泡的一些宏觀參數(shù)如基圓半徑和氣泡生長速率等則較容易獲得。因此，本研究希望通過實驗中測得的宏觀氣泡動態(tài)參數(shù)來確定微液層與大液層的蒸發(fā)在沸騰傳熱過程中的影響。在此基礎上，通過建立沸騰界面?zhèn)鳠崮Ｐ蛠眍A測低熱通量條件下銅表面上沸騰過程中氣泡底部的微液層厚度，從而為進一步實現(xiàn)完整的沸騰過程模擬提供理論支持。

1 實驗裝置與數(shù)據(jù)處理

1.1 實驗裝置

圖2 單氣泡池沸騰實驗裝置Fig.2 Schematic diagram of the experimental apparatus for single bubble boiling system

池沸騰實驗裝置示意圖如圖2所示。其中沸騰池為有機玻璃制成的透明容器，其尺寸為160 mm×160 mm×210 mm。容器頂部預留有排汽孔，容器內(nèi)氣壓與環(huán)境保持相同。采用去離子水作為沸騰工質(zhì)，實驗前首先通過真空泵去除溶解在去離子水中的不凝氣。主加熱器置于容器底部，材質(zhì)為銅，其上表面為加熱表面，與去離子水直接接觸，直徑為10 mm。銅芯周圍由聚四氟乙烯(PTFE)包裹以減少熱量損失。銅芯底部固定有電阻加熱棒，通過控制電阻棒兩端電壓可以實現(xiàn)不同的加熱功率。銅芯內(nèi)部不同深度(2.5 mm和7.5 mm)處固定有T形熱電偶，據(jù)此可以對加熱表面溫度和熱通量進行預估[27]。為了降低加熱表面粗糙度以減少沸騰過程中氣泡間的干涉，實驗前首先采用不同型號砂紙(p180～p3000)對加熱表面進行逐級打磨，隨后在光滑加熱表面中心加工單一孔穴以便產(chǎn)生單一汽化核心。當加熱器溫度逐漸升高時，加熱面中心位置處的汽化核心會率先活化產(chǎn)生氣泡。

實驗過程中首先通過輔助加熱器將去離子水加熱至接近飽和狀態(tài)。隨后打開主加熱器，通過電阻棒對銅芯進行加熱。實時監(jiān)測加熱棒內(nèi)溫度分布，直至2 min內(nèi)加熱棒內(nèi)溫度變化不超過0.1℃認為達到穩(wěn)態(tài)。逐級、緩慢調(diào)整加熱功率至加熱表面上汽化核心處有氣泡產(chǎn)生，通過溫度采集儀和高速相機分別記錄加熱表面溫度變化和氣泡生長圖像。繼續(xù)調(diào)整主加熱器兩端電壓即可得到不同熱通量下沸騰過程中加熱表面的氣泡生長過程。在整個過程中，通過熱電偶實時監(jiān)測沸騰池內(nèi)主流體的溫度，調(diào)節(jié)輔助加熱器功率，使池內(nèi)去離子水維持在飽和狀態(tài)。

本實驗系統(tǒng)中所用高速相機型號為Photron Fastcam Mini UX50,其在最大分辨率1280×1024下最大拍照速度為2000 fps，降低分辨率后其速度最高可達160000 fps，曝光時間為1.01～3.9μs。所匹配的顯微鏡頭為體視顯微鏡，型號為Nikon SMZ800N，其放大倍數(shù)為5～480倍。T型熱電偶溫度采用Agilent 34970A結合34901A數(shù)據(jù)采集模塊進行讀取，響應速度10-5s,采用單通道測溫時其最大掃描頻率為60次/秒，多通道測溫時掃描頻率會隨之降低。

1.2 氣泡生長圖像序列處理

為獲取加熱表面上氣泡動態(tài)參數(shù)的變化，需首先對由高速相機獲得的圖片序列進行預處理以準確提取氣泡邊緣。如圖3(a)所示，通過去除背景和高斯濾波后得到背景較為光滑連續(xù)的圖像；然后調(diào)整對比度增強背景與氣泡間的差異后進行圖像二值化處理即可得到較為清晰的氣泡形態(tài)；進一步將由光照引起的高亮部分進行填充，再進行邊緣檢測即可得到完整的氣泡輪廓。圖3(b)給出了氣泡輪廓與原圖像的比對結果，可以發(fā)現(xiàn)，本方法較好地還原了孤立氣泡在加熱表面上的生長過程。

1.3 誤差分析

本實驗中氣泡體積計算的主要誤差來源于氣泡形狀提取時的圖像處理。氣泡生長圖像中610個像素點代表了真實加熱表面上10 mm的區(qū)域，因此每個像素點對應的長度應為0.016 mm/pix，在氣泡邊界提取過程中引入的誤差不超過3個像素點[28]即0.048 mm。氣泡在脫離時的直徑約為1.9 mm，高度為2.2 mm。氣泡的體積采用微元法根據(jù)二維氣泡形狀計算得到，如圖4所示，將氣泡分割成N個微圓柱體(N=Hb/Δh)，其中Δh取一個像素的高度，將所有微圓柱體體積Vi疊加得到整個氣泡的體積，即：

圖3 氣泡輪廓提取流程以及所提取到的氣泡輪廓與原始圖像的對比Fig.3 Theextracting procedure of bubble outlines and the comparison of the obtained bubble outlines with the original images

圖4 基于微元法的氣泡體積計算示意圖Fig.4 The measurement of the bubble volume from the bubble growing images

根據(jù)多參數(shù)單次間接測量誤差傳遞公式[29]，間接測量得到的氣泡體積的相對誤差約為：

氣泡基圓半徑由氣泡圖像序列獲得，平均值為0.2 mm，因此，其測量相對誤差為ΔRb/Rb=0.048/0.2=24%。提取基圓半徑時，由于采用同一分割閾值對一個氣泡周期內(nèi)的全部時序列圖像進行處理，這一測量誤差不會影響其在一個周期內(nèi)的位置變化趨勢。

加熱器內(nèi)的溫度測量采用T型熱電偶，測溫不確定度為0.1 K。銅芯內(nèi)的總熱通量可根據(jù)其兩端電壓、電流和直徑得到，其相對誤差為[27]：

在后續(xù)模擬中，固體加熱器底部采用實驗中測得的溫度作為定溫邊界條件，熱通量的相對誤差對沸騰實驗有影響，但不影響后續(xù)的模擬分析。

2 實驗結果與模擬分析

2.1 典型氣泡生長周期內(nèi)的氣泡動態(tài)參數(shù)的變化

實驗中首先對加熱表面上的孤立氣泡生長過程進行了觀測。利用高速相機從沸騰池側面對氣泡的動態(tài)行為記錄，拍照頻率為3200 fps，圖像分辨率為1/62 mm/pix。打孔前光滑加熱表面粗糙度Ra約為0.15μm。采用去離子水作為沸騰工質(zhì)。

圖5所示為一個典型周期內(nèi)氣泡生長過程中的氣液界面位置。加熱表面過熱度約為4.82 K，熱通量5.5 kW/m2。在此工況下，氣泡脫離直徑為1.89 mm，氣泡周期為46 ms，無等待周期。可以看出在整個氣泡生長過程中，基圓半徑的變化非常顯著。由于氣泡底部基圓半徑的變化對加熱表面換熱特性有著重要的影響[6]，研究中根據(jù)氣泡底部基圓半徑的變化，將氣泡的生長過程分為基圓半徑迅速擴張-維持-迅速收縮三個階段。圖6(a)給出了氣泡生長過程中高度和直徑的變化?？梢钥吹剑跉馀萆L前8 ms內(nèi)，由于氣泡底部基圓半徑迅速向外擴張，氣泡在橫向和縱向均迅速增大且尺寸相當，形狀更接近半圓形，這段時間約占氣泡總周期的20%。在基圓半徑相持階段(8～40 ms)，氣泡底部與加熱表面接觸線基本變化不大，氣泡在縱向和橫向的生長速度均有所減緩[圖5(b)]；20 ms以后，氣泡在縱向的尺寸略大于橫向，縱向拉伸逐漸明顯。當氣泡臨近脫離時，基圓半徑在極短的時間內(nèi)(3 ms)迅速縮減為0，氣泡呈現(xiàn)明顯的倒立水滴狀，氣泡半徑基本不變，而縱向尺寸拉伸到最大，直至脫離加熱表面[圖5(c)]。

圖5 氣泡生長過程中不同階段的氣泡形狀及底部基圓半徑的變化Fig.5 The variation of vapor-liquid interface and bubble root during the bubble growing period

在此過程中，隨著氣泡受力情況的不斷變化氣泡的形狀不斷改變，其中氣泡的脫離過程主要由氣泡受到的浮力Fb和由表面張力引起的壁面附著力Fσ決定[24]：

由式(6)可看出，壁面能提供的最大附著力主要由氣泡底部基圓半徑Rb和接觸角θ的變化決定。圖6(b)、(c)分別給出了氣泡體積和氣泡底部與加熱表面接觸角的變化。在氣泡生長初期8 ms內(nèi)接觸角變化較快，由120°左右迅速減小至40°左右，隨后保持在35°左右小范圍波動，直到氣泡臨近脫離時減小到約30°。而氣泡體積在氣泡生長初期增長較快，8 ms時速率達到最大，這意味著氣泡生長初期氣液相變較多，而此時氣液宏觀接觸面積較小[圖5(a)]，因此相變主要發(fā)生在實驗較難觀測到的微液層區(qū)域。隨后氣泡體積增速逐漸減緩直至脫離時氣泡體積達到最大，約3.7 mm3。根據(jù)式(5)和式(6)計算得到的氣泡受到的浮力和壁面能提供的最大附著力如圖6(d)所示。可以看到，在43 ms之前，氣泡受到的浮力始終小于壁面能提供的最大附著力，因此氣泡附著在壁面生長。43 ms時，氣泡受到的浮力與壁面附著力持平。隨后，由于氣泡體積繼續(xù)增大，氣泡受到的浮力隨之增大，因此氣泡被迅速拉高。同時氣泡底部基圓半徑和接觸角也隨之迅速減小，氣泡受到的壁面附著力也隨之減小，因此氣泡在很短的時間內(nèi)脫離加熱表面。

2.2 基于氣泡動態(tài)參數(shù)的近壁面處相變機理分析

池沸騰過程中界面處的傳熱包含多種傳熱機理，包括微液層的蒸發(fā)，大液層以及熱邊界層外氣液界面的蒸發(fā)和自然對流等[30]。自然對流一般發(fā)生在氣泡投影面積外側，而相變過程則主要發(fā)生在氣泡底部的微液層和大液層區(qū)域。氣泡基圓半徑定義為氣泡中心到氣泡底部氣液界面與加熱表面接觸點的距離。所在位置為微液層與大液層的分界點，如圖1所示?；鶊A半徑Rb的變化可以從實驗中直接測得，而熱邊界層厚度δt則可用式(7)計算[31]：

式中，g為當?shù)刂亓铀俣龋蚻、αl、β分別為液體工質(zhì)運動黏度、熱擴散率、體積膨脹系數(shù)，Tw和Tsat分別為加熱表面溫度和工質(zhì)飽和溫度。

圖6 氣泡生長過程中氣泡高度、直徑、體積、接觸角以及浮力和表面張力隨時間的變化(ΔT=4.82 K)Fig.6 The variations of bubble height,diameter,volume,contact angle,buoyancy and surface tension with time during the bubble growth period

圖7 實驗中測得的氣泡生長速率、氣泡底部基圓半徑和大液層區(qū)域隨時間的變化Fig.7 Comparison of the bubble growth rates with the bubble root radii and macrolayer regions measured in the experiment

為進一步探究沸騰過程中尤其是在氣-液-固三相共存的熱邊界層內(nèi)的相變機理，實驗中對三組不同氣泡生長過程中的基圓半徑Rb和大液層區(qū)域(Rt-Rb)的變化進行了分析，如圖7所示?？梢钥吹剑M數(shù)據(jù)中氣泡基圓半徑均在氣泡生長初期迅速增大，隨后呈現(xiàn)一定的波動，而在氣泡脫離前則快速減?。欢笠簩訁^(qū)域主要與熱邊界層內(nèi)的氣泡尺寸相關：隨著氣泡體積的逐漸增大，大液層區(qū)域在氣泡的整個生長過程中均呈現(xiàn)增大的趨勢，在氣泡脫離前達到最大。可以看到，氣泡生長過程中微液層區(qū)域與大液層的變化趨勢截然不同。因此，通過對比分析氣泡生長速率與實驗中氣泡基圓半徑、大液層區(qū)域隨時間的變化可定性地判斷微液層區(qū)域與大液層區(qū)域在氣泡生長中的貢獻，即，倘若大液層的蒸發(fā)是氣泡生長最主要的蒸汽來源，氣泡的生長速率應該與其變化規(guī)律基本一致；若微液層蒸發(fā)占主導，則氣泡的生長速率應該與基圓半徑變化規(guī)律相近。

進一步，通過實驗中氣泡體積變化可以得到氣泡的生長速率V′，如圖7(a)所示。對比分析氣泡生長速率與實驗中氣泡基圓半徑、大液層區(qū)域隨時間的變化，可以看到，氣泡生長速率并不隨著大液層區(qū)域的增大而增大，而與氣泡基圓半徑的變化趨勢基本相同：當氣泡基圓半徑增大時，氣泡生長速率隨之增大，反之亦然。進一步利用相關性公式(8)計算氣泡生長速率與基圓半徑和大液層半徑的相關性，可以得到圖7的三個周期中氣泡生長速率與大液層區(qū)域變化兩組數(shù)據(jù)之間的相關系數(shù)分別為-0.1965、-0.3195和-0.4554，其中負號表示二者呈現(xiàn)負相關；而與氣泡基圓半徑變化之間的相關系數(shù)分別為0.8017、07654和0.7451。

因此，可認為氣泡生長速率的變化與基圓半徑的變化呈現(xiàn)顯著的正相關，這表明孤立氣泡沸騰過程中的相變換熱主要發(fā)生在基圓半徑內(nèi)側即微液層區(qū)域內(nèi)，而大液層區(qū)域的蒸發(fā)在當前實驗工況下則貢獻較少。

2.3 基于孤立氣泡沸騰過程中宏觀參數(shù)的微液層厚度預測

2.3.1 沸騰換熱模型 既然氣泡底部的微液層蒸發(fā)是孤立氣泡沸騰過程中最重要的換熱機制，其厚度動態(tài)變化則成為理解沸騰過程的重要參數(shù)。由于其尺度相對較小(一般小于50μm)，厚度通常難以直接觀測，實驗中通常借助激光干涉等技術對其進行測量[18-19]。本文嘗試結合數(shù)值模型，通過匹配實測的氣泡生長速率，來對氣泡生長過程中底部微液層的厚度進行預測。由于加熱器關于中心軸線嚴格對稱，汽化核心位于其頂面中心處，計算中采用5°的二維軸對稱楔形區(qū)域作為計算區(qū)域，如圖8(a)所示。區(qū)域的半徑設置為5 mm，高7.5 mm。加熱器內(nèi)熱傳導控制方程為：

式中，T為加熱器內(nèi)溫度分布，αw為固體的熱擴散系數(shù)，對于銅材料為1.14×10-4m2/s。計算域內(nèi)側為旋轉軸，模擬中采用對稱邊界，外側為隔熱層，忽略其換熱。底部根據(jù)實驗結果設置為固定溫度。計算域上表面為加熱表面。汽化核心正下方2.5 mm處設置有溫度觀測點，模擬中用此點的計算溫度與加熱器內(nèi)的溫度變化實測值進行匹配。

圖8 計算區(qū)域和沸騰換熱機理示意圖Fig.8 Schematics of the calculated domain and the heat transfer mechanisms on the boiling surface

對于加熱表面上熱邊界層內(nèi)的換熱，根據(jù)實驗中氣泡界面在不同時刻的不同位置，模擬中將其劃分為三個區(qū)域,由氣泡中心向外依次為干斑區(qū)、微液層區(qū)和自然對流區(qū)，如圖8(b)所示。其中，氣泡基圓半徑外側為自然對流區(qū)，基圓半徑的位置則由實驗測量結果確定。而在基圓半徑Rb以內(nèi)，加熱表面由微液層覆蓋，微液層區(qū)域中心最內(nèi)側為干斑區(qū)域。Gao等[18]通過激光干涉法對乙醇氣泡沸騰過程中干斑區(qū)的演化進行了測量，結果表明，在微液層蒸干前，干斑區(qū)域半徑隨著氣泡生長時間的推移呈現(xiàn)對數(shù)變化，其他相關實驗[17,19]與模擬結果[32]也基本支持這一結論。同時在氣泡生長后期，由于氣泡底部基圓半徑發(fā)生收縮，三相區(qū)域內(nèi)會出現(xiàn)較為明顯的橫向回流[6,33]，導致加熱表面被重新潤濕，如圖5(c)所示，干斑區(qū)域隨之減小。根據(jù)Jiang等[24]的模擬結果，液層的回流速度與基圓半徑的收縮速度基本呈正比例變化，因此，模型中液層回退的距離為：

式中，Rb,0和Rd,0分別為上一時刻的基圓半徑和干斑區(qū)域半徑。當氣泡脫離時，干斑區(qū)域半徑與基圓半徑重合，微液層完全蒸干，微液層區(qū)域為0。據(jù)此，可以假設干斑區(qū)域變化為：

式中，τd為氣泡脫離時刻，tshrink為基圓半徑開始迅速收縮的時刻，可根據(jù)氣泡受到的浮力與壁面的附著力來確定[圖6(d)]；Cdry和m為根據(jù)實驗數(shù)據(jù)擬合得到的經(jīng)驗參數(shù)，其取值與工質(zhì)物性、壁面過熱度、壁面潤濕性等因素有關。對于m的取值，Gao等[18]給出的參考值為0.6～1.2之間。微液層存在區(qū)域為氣泡底部三相接觸線到基圓半徑之間，因此微液層區(qū)域半徑為：

模擬中加熱器上表面的熱通量根據(jù)各區(qū)域上的換熱狀況確定。在氣泡中心最內(nèi)側的干斑區(qū)域，加熱表面直接與過熱蒸汽相接觸，換熱量很小，忽略不計。而在微液層蒸發(fā)區(qū)，加熱表面的熱量通過熱傳導進入微液層，液層在氣液交界面位置蒸發(fā)產(chǎn)生蒸汽。根據(jù)熱平衡原理，微液層蒸發(fā)消耗的熱量與熱傳導進入氣泡的熱量相同，可根據(jù)傅里葉導熱定律計算[34]：

式中，Tw為加熱表面溫度，Tsat為飽和液體溫度。而在微液層外側，由于氣泡與加熱表面之間的液體較厚，熱量從加熱表面直接傳遞進入主流體中，其熱通量根據(jù)Holman自然對流關系式評估[35]：

式中，Tw(r,t)為t時刻加熱表面上r處的溫度；Rd為干斑區(qū)域半徑，即三相接觸線位置；Rb為氣泡基圓半徑，即微液層最外側。

由于干斑區(qū)域和自然對流區(qū)域不發(fā)生相變，因此這個過程中產(chǎn)生的蒸汽量(由此可計算氣泡生長率)可根據(jù)微液層的蒸發(fā)量來確定，即[36]：式中，Hfg為汽化潛熱，ρv為蒸汽密度。

2.3.2 微液層厚度的預測 由于微液層的分布在目前實驗中無法直接測定，根據(jù)陳宏霞等[12]的實驗結果，計算中初步假定其厚度均勻。其取值需要在計算中進行多次迭代匹配來確定，計算過程如圖9所示。模擬前首先給定微液層厚度的初始值，范圍為1～5μm。根據(jù)預設的微液層厚度結合實驗中測定的氣泡動態(tài)參數(shù)通過式(15)可計算得到加熱表面的熱通量，以此作為加熱表面邊界條件進行加熱器的導熱計算則可獲得當前時刻的加熱表面溫度分布。迭代多個周期直至計算中加熱表面溫度達到穩(wěn)定時，根據(jù)式(16)可計算得到一個周期內(nèi)的蒸汽產(chǎn)生量和氣泡生長率并與實驗所獲氣泡生長率進行比較，若與實驗不符合，則調(diào)整液層厚度進行下一迭代步的計算。氣泡生長率誤差函數(shù)設置為：

圖9 加熱表面?zhèn)鳠嵊嬎闩c實驗測量參數(shù)間的迭代匹配策略Fig.9 The matching strategy between the simulation of boiling surface heat transfer with the measured bubble dynamics

根據(jù)誤差函數(shù)計算結果，采用二分法調(diào)整預設的微液層厚度重新進行多個周期的迭代計算，直到計算得到的蒸汽產(chǎn)生量與實驗結果間的誤差函數(shù)小于預設的最小值ε。另一方面，也比較模擬中熱電偶測點位置處的溫度與實驗測得的結果，當誤差在0.1 K以內(nèi)時認為與實驗結果相符，迭代完成，此時的微液層厚度即為當前實驗條件下的預測結果。2.3.3 模擬結果 基于OpenFOAM平臺3.0.x版本下的laplacianFoam求解器對加熱器內(nèi)部溫度場[式(9)]進行求解。采用有限體積法對控制方程進行離散，非穩(wěn)態(tài)項采用隱式離散，擴散項采用中心差分進行離散，矩陣求解選用預處理共軛梯度法。計算開始前首先將加熱器上表面邊界條件設置為自然對流，熱通量由式(14)決定，下表面過熱度根據(jù)實驗測量值固定為4.95 K。時間步長設為0.3125 s。

圖10 不同網(wǎng)格數(shù)量下熱電偶處的過熱度變化Fig.10 The local superheat at the position of thermal couple T TC in simulations with different mesh sizes

分別采用50×60、125×150和250×300三種不同數(shù)目的網(wǎng)格對加熱表面上的自然對流換熱過程進行了模擬，其對應最小網(wǎng)格尺寸分別為0.1、0.04和0.02 mm。圖10給出了三種網(wǎng)格下模擬中熱電偶位置處的過熱度TTC的變化?？梢钥吹剑斁W(wǎng)格尺寸為0.1 mm時，熱電偶處的過熱度略小于網(wǎng)格加密后的計算結果，而當網(wǎng)格尺寸為0.04 mm時，計算結果不再隨著網(wǎng)格的加密而發(fā)生變化。因此，模擬中最終將計算區(qū)域劃分為125×150個網(wǎng)格。此外，在1.0 s時過熱度逐漸趨于平緩，從0.9 s到1.0 s過熱度降低小于10-3K，因此認為此時換熱過程已經(jīng)達到穩(wěn)態(tài)。模擬中將此時的加熱器內(nèi)的溫度分布作為初始場進行加熱表面的換熱計算。

圖11 氣泡底層基圓半徑、干區(qū)和微液層區(qū)域半徑變化的模擬結果Fig.11 The simulated variations of bubble root,dry area and microlayer region radiiwith time

根據(jù)干斑區(qū)域計算公式(11)，結合所測一個周期內(nèi)氣泡基圓半徑，可得干斑和微液層區(qū)域的變化如圖11所示。其中右上角插圖給出了氣泡基圓半徑Rb，干斑區(qū)域半徑Rd和微液膜區(qū)域半徑Rmic的位置關系。經(jīng)過多次迭代計算后，模擬中最終得到的微液層厚度為3.43μm，與文獻[8,18-19]的結果基本吻合。此時計算得到的一個周期內(nèi)氣泡生長速率如圖12所示，與實驗結果的相對誤差為22%?？梢钥吹?，與實驗結果相比，模擬得到的氣泡生長速率基本合理。這也進一步表明微液層在孤立氣泡的沸騰過程中起主導作用。

圖12 氣泡生長速率的模擬結果和實驗結果的對比Fig.12 The comparison of the simulated and experimental bubble growth rate

圖13則給出了模擬過程中加熱表面的平均過熱度以及熱電偶位置處的過熱度隨計算時間的變化。在經(jīng)過20個氣泡周期的換熱計算后，加熱表面溫度變化基本達到穩(wěn)定，平均過熱度約為4.82 K。由圖13(b)可以看出，穩(wěn)定后加熱表面在一個周期內(nèi)的溫度變化與微液層區(qū)域的變化密切相關：當微液層區(qū)域增大時，加熱表面平均溫度降低；而在氣泡脫離前，隨著微液層區(qū)域的迅速減小，壁面溫度迅速升高。模擬中得到的汽化核心下方2.5 mm處熱電偶的過熱度TTC基本在4.87～4.88 K范圍內(nèi)波動，與實驗參考值4.88 K基本相符。

圖14給出了氣泡生長不同階段內(nèi)不同時刻加熱表面上的過熱度分布?？梢钥闯?，在氣泡生長初期的8 ms之內(nèi)[圖14(a)]，加熱表面溫度基本呈現(xiàn)下降趨勢，而在氣泡生長后期[圖14(b)]逐漸恢復到初始狀態(tài)。其中，加熱表面中心位置溫度波動最為劇烈：在前8 ms內(nèi)，由4.83 K迅速降低到4.28 K；隨后由于干斑區(qū)域的出現(xiàn)，其溫度逐漸回升。相比之下，距離汽化核心較遠的位置則波動較小。圖中“○”表示該時刻的三相線位置，其內(nèi)側為干斑區(qū)，外側為微液層區(qū)域；而“◆”則表示該時刻的基圓半徑位置，其外側為自然對流區(qū)域?？梢钥吹?，加熱表面上的最低溫度出現(xiàn)在兩標記之間的區(qū)域即微液層覆蓋區(qū)；在更內(nèi)側的干斑區(qū)域，加熱表面與過熱蒸汽直接接觸，因此溫度較高；而在外側的自然對流區(qū)，由于無相變發(fā)生，換熱較差。因此，氣泡底部最終形成以三相接觸線和氣泡基圓半徑為邊界的環(huán)形低溫帶[11,37]。

模擬中通過與實驗得到的氣泡生長速率進行匹配，間接得到的氣泡底層微液層厚度為3.43μm。其他研究者直接測量得到的微液層厚度分別為:1～4μm(Zou等[8])，3μm(Jung等[19])，4μm(Gao等[18])，本研究得到的微液層厚度在合理范圍。本文的預測結果進一步驗證了微液層蒸發(fā)機理。對于銅表面上初始液層厚度斜率Cslope[16]的取值，在過熱度為5 K左右時，本文采用的參考值為0.02～0.03，與Sato等[25]的推薦值0.0245基本相符。

圖13 加熱表面及加熱表面以下2.5 mm處過熱度隨時間變化的模擬結果(a)和一個周期內(nèi)的加熱表面過熱度隨時間的變化(b)Fig.13 The simulated variations of the mean superheat on the heated surface and thelocal superheat at the position of 2.5 mm below the surface(a)and the mean superheat of heated surface within one bubble period(b)

圖14 氣泡生長過程中不同時刻加熱表面過熱度分布的模擬結果及氣泡界面位置的實驗結果Fig.14 The superheat distributions of the heated surface in the simulation and the positions of bubble interface of the experiment at different moments

需要說明的是，本文研究僅針對較低熱通量(約5.5 kW/m2)下的孤立氣泡沸騰區(qū)域，在此區(qū)域內(nèi)，微液層的蒸發(fā)在沸騰換熱過程中起主導作用。但隨著熱通量增大，微液層蒸發(fā)作用在氣泡生長過程中所占份額會逐漸減少，而大液層以及主流體中兩相界面部分的蒸發(fā)所占的份額會逐漸增加。在氣泡水平合并頻繁的高熱通量區(qū)域(1.0 MW/m2以上)大液層的蒸發(fā)將成為換熱的主要機理[9,38]。不同熱通量區(qū)域的沸騰換熱機理和氣泡動力學特性還需通過實驗和多物理場耦合分析進一步考察。

3結 論

本文首先通過單個氣泡沸騰實驗獲得了氣泡生長過程時間序列圖像，進一步對氣泡生長過程圖像進行圖像處理，獲得了氣泡生長過程中的接觸角、基圓半徑、體積、高度和直徑等動態(tài)參數(shù)變化。相關性分析結果表明，氣泡的生長速率與氣泡底層基圓半徑的變化顯著相關，由此可認為基圓半徑內(nèi)側的微液層蒸發(fā)是孤立氣泡沸騰過程中最主要的相變機理。在此基礎上，結合實驗測得的氣泡動態(tài)參數(shù)建立了加熱表面的換熱分析模型，通過迭代匹配實驗測得的氣泡生長速率和加熱表面過熱度最終得到，當加熱表面過熱度為4.82 K時，氣泡底層微液層的厚度約為3.43μm。液層厚度預測結果與相關實驗測量值吻合，表明低熱通量孤立氣泡沸騰過程中微液層蒸發(fā)是氣泡生長的主要機理，這為進一步進行沸騰過程流熱耦合分析和沸騰換熱表面設計奠定了理論基礎。

符號說明

Cdry——干斑區(qū)域演變經(jīng)驗參數(shù)，m

Cslope——微液層初始斜率

cp,l——液體比定壓熱容，J/(kg·K)

D——氣泡橫向尺寸，m

Db——氣泡直徑，m

Eresi——迭代殘差

Fb——浮力，N

Fσ——壁面附著力，N

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

Hb——氣泡高度，m

Hfg——汽化潛熱，J/kg

I——通過加熱棒的電流，A

Lmic——基圓半徑收縮引起的微液層移動的距離，m

m——經(jīng)驗參數(shù)，取值范圍0.6～1.2

Rb——基圓半徑，m

Rd——干斑區(qū)域半徑，m

Rmic——微液層區(qū)域半徑，m

Rt——熱邊界層與氣泡界面接觸點位置，m

q——平均熱通量，W/m2

qboi——加熱表面的局部熱通量，W/m2

qmic——微液層蒸發(fā)的熱通量，W/m2

qnc——自然對流的熱通量，W/m2

r——距離氣泡中心的距離，m

T——溫度，K

Ts——加熱表面平均過熱度，K

Tsat——工質(zhì)飽和溫度，K

TTC——熱電偶位置過熱度，K

ΔT——加熱表面過熱度，K

t——時間，s

tshrink——基圓半徑開始迅速收縮的時刻，s

U——加熱棒兩端電壓，V

V——氣泡體積，m3

V′——氣泡生長速率，m3/s

z——高度，m

α——熱擴散系數(shù)，m2/s

β——體積膨脹系數(shù)，1/K

γevap——蒸汽產(chǎn)生率，m3/s

δt——熱邊界層厚度，m

δ0——初始微液層厚度，m

ε——迭代相對誤差

θ——接觸角，(°)

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

μ——動力黏度，N·s/m2

ν——運動黏度，m2/s

ρ——密度，kg/m3

ρxy——相關系數(shù)

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

τd——氣泡脫離時刻，s

下角標

l——液體

v——蒸汽

w——固體表面

0——上一時刻參數(shù)