含多個(gè)矩形加熱器通道內(nèi)流動(dòng)沸騰傳熱性能的介觀數(shù)值方法研究*

李迎雪， 王浩原， 婁 欽

（1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

引 言

隨著現(xiàn)代科技的發(fā)展，設(shè)備的散熱問(wèn)題愈發(fā)受到關(guān)注，例如核反應(yīng)堆冷卻和芯片散熱等[1-3].沸騰作為一種有效的散熱方式，在較低溫度下即可實(shí)現(xiàn)較高的傳熱效率，因此得到了廣泛的應(yīng)用.沸騰可以分為池沸騰和流動(dòng)沸騰兩種，和池沸騰相比，流動(dòng)沸騰中常伴隨著更劇烈的對(duì)流傳熱及相變速率，從而可以達(dá)到更高的換熱效率，因此流動(dòng)沸騰受到了越來(lái)越多的關(guān)注[4].

在流動(dòng)沸騰中，熱源對(duì)相變過(guò)程及傳熱效率有著至關(guān)重要的影響[5-7]，因此許多研究者研究了不同加熱器對(duì)流動(dòng)沸騰的影響.例如，Gong和Cheng[8]在水平管道內(nèi)放置一微型加熱器，研究了不同入口速度、接觸角和加熱器尺寸對(duì)飽和流動(dòng)沸騰的影響，他們發(fā)現(xiàn)隨著入口速度的增加或加熱器長(zhǎng)度的減小，流型從環(huán)狀流變?yōu)槎稳?，并且在段塞流階段，通道中的氣泡段塞長(zhǎng)度減小.Tan等[9]研究了水平管內(nèi)壁面潤(rùn)濕性對(duì)流動(dòng)沸騰性能的影響，比較了親水性表面、疏水性表面和梯度潤(rùn)濕性三種情況，研究發(fā)現(xiàn)梯度潤(rùn)濕性可以結(jié)合接觸角的連續(xù)變化效應(yīng)以及親水和疏水表面的優(yōu)點(diǎn)，限制段塞流的出現(xiàn)，使通道內(nèi)更多的區(qū)域表現(xiàn)為氣泡流狀態(tài).Yin等[10]使用流體體積模型和氣液相變模型模擬了分離型熱管蒸發(fā)器中垂直管內(nèi)的流動(dòng)沸騰，研究了低熱流密度下管道內(nèi)不同高度的氣泡形態(tài)變化，發(fā)現(xiàn)沸騰的流動(dòng)模式隨管的高度而變化.可以看出，加熱器垂直或水平的放置方式對(duì)流動(dòng)沸騰中氣泡運(yùn)動(dòng)和壁面換熱有著明顯的影響.因此，Azizifar等[11]通過(guò)實(shí)驗(yàn)研究了水的過(guò)冷流動(dòng)沸騰過(guò)程，結(jié)果表明，與普通管相比，金屬泡沫導(dǎo)致沸騰傳熱系數(shù)（heat transfer coefficient，HTC）顯著增加，金屬泡沫管放置方式不同時(shí)，水平放置的傳熱效率比垂直放置高10%.

上述流動(dòng)沸騰的研究中均采用了單個(gè)加熱器加熱方式，最近有研究指出[12]，多個(gè)加熱器加熱時(shí)，會(huì)產(chǎn)生多個(gè)氣泡成核位點(diǎn)，氣泡間的相互作用會(huì)導(dǎo)致部分成核位點(diǎn)失活從而降低傳熱效率.因此，多個(gè)加熱器引起的流動(dòng)沸騰逐漸受到關(guān)注.Vontas等[13]模擬了多個(gè)成核位點(diǎn)和同一成核位置發(fā)生多次重復(fù)成核的情況,觀察到表面潤(rùn)濕性對(duì)換熱系數(shù)具有顯著影響，發(fā)現(xiàn)親水和疏水表面的主要傳熱機(jī)制分別是液膜蒸發(fā)和接觸線蒸發(fā).Sun等[14]研究了由多個(gè)加熱器引起的流動(dòng)沸騰，通過(guò)氣泡成核、生長(zhǎng)、離開、聚結(jié)的行為，討論了單相流、氣泡流、段塞流和核態(tài)沸騰轉(zhuǎn)折點(diǎn)（departure from nucleate boiling， DNB）多種流動(dòng)模式的行為.同時(shí)研究了熱通量、Reynolds數(shù)、流動(dòng)通道寬度、成核點(diǎn)寬度等因素對(duì)流動(dòng)沸騰的影響，特別是對(duì)DNB點(diǎn)的影響.Yadav等[15]對(duì)多個(gè)加熱器引起的流動(dòng)沸騰進(jìn)行了空隙率的局部測(cè)量，觀察到低入口過(guò)冷度比高入口過(guò)冷度時(shí)的局部空隙率更高.此外，改變加熱器的位置可以改變觀察孔隙率的徑向位置，增加加熱器的數(shù)量可以增加孔隙率.但是，他們工作中并沒有考慮多個(gè)加熱器之間的相互影響.Baltis等[16]研究了垂直排列的加熱器之間的相互影響，發(fā)現(xiàn)上游的加熱器通過(guò)額外對(duì)流可以增加下游加熱器的氣泡成核頻率和氣泡脫離直徑，隨著加熱器間距的增加，這種促進(jìn)作用越來(lái)越小.此外，上游加熱器產(chǎn)生的成核氣泡會(huì)抑制下游的氣泡成核，這種抑制趨勢(shì)的影響隨著下游氣泡成核頻率的增加而減弱.從上述研究可以看出，關(guān)于多加熱器引起的流動(dòng)沸騰研究還存在一些不足，不同加熱器表面氣泡去除的動(dòng)力學(xué)行為差異和不同加熱器表面的傳熱性能差異尚不明確，例如，當(dāng)加熱器的長(zhǎng)度不同或相鄰兩個(gè)加熱器之間距離不同的情況下，氣泡動(dòng)力學(xué)行為差異和加熱器表面的熱流密度差異.

格子Boltzmann方法（lattice Boltzmann method, LBM）是20世紀(jì)80年代中期建立和發(fā)展起來(lái)的一種模擬方法，它是介于流體的微觀分子動(dòng)力學(xué)模型和宏觀連續(xù)模型之間的介觀數(shù)值模型，它兼具二者的優(yōu)點(diǎn)，可以直觀、方便地描述流體與周圍環(huán)境之間的相互作用.近年來(lái)該方法經(jīng)常被用于研究滲流、多相流、微流體等流動(dòng)問(wèn)題[17-19].因此，本文基于格子Boltzmann方法模擬了通道內(nèi)多加熱器下的流動(dòng)沸騰，詳細(xì)研究了不同矩形加熱器（包括矩形加熱器的間距、矩形加熱器的尺寸及其表面潤(rùn)濕性的不同）對(duì)氣泡動(dòng)力學(xué)行為和傳熱性能的影響.

1 數(shù)值方法

1.1 偽勢(shì)LB兩相流模型

本文采用Gong和Cheng提出的改進(jìn)偽勢(shì)LB兩相流模型來(lái)模擬通道中的流動(dòng)沸騰[4].模型包含兩個(gè)粒子分布函數(shù)的演化方程，一個(gè)是密度分布函數(shù)，另一個(gè)是溫度分布函數(shù).接下來(lái)，給出了偽勢(shì)LB兩相流模型的簡(jiǎn)要描述.

密度分布函數(shù)的演化方程如下：

其中g(shù)s代 表流固之間的相互作用力強(qiáng)度；s(x)是一個(gè)開關(guān)系數(shù)，當(dāng)格點(diǎn)x為流體點(diǎn)時(shí)，其值為0，此時(shí)對(duì)應(yīng)的Fads為 0，即表示不存在流體與固體間的相互作用力；當(dāng)格點(diǎn)x為固體點(diǎn)時(shí)，s(x)的值為1.

其他外部作用力Fext的表達(dá)式為

其中g(shù)表示重力加速度，ρa(bǔ)ve是單個(gè)時(shí)間步長(zhǎng)整個(gè)計(jì)算區(qū)域內(nèi)流體的平均密度.

流體密度ρ和速度U的表達(dá)式為

在本文工作中，采用二維九速（D2Q9）數(shù)值模型進(jìn)行研究.在D2Q9模型中對(duì)應(yīng)的權(quán)重系數(shù)為 ω0=4/9；ωi=1/9 ，i=1, 2, 3, 4；ωi=1/36，i=5, 6, 7, 8.其離散速度模型為

1.2 能量方程

和密度分布函數(shù)類似，本文的溫度場(chǎng)由溫度分布函數(shù)得到，其演化方程如下：

其中g(shù)i(x,t)表示粒子在x位置t時(shí)刻i方 向的溫度分布函數(shù)； τT是 溫度的無(wú)量綱松弛時(shí)間，它的值與熱擴(kuò)散率α有關(guān)，為平衡態(tài)溫度分布函數(shù)，其具體表達(dá)式為

式（12）中的? 代表的是源項(xiàng)，可以表示為

此處的T是宏觀溫度，可以通過(guò)溫度分布函數(shù)得到：

2 模型與驗(yàn)證

2.1 計(jì)算模型

本文采用相變LB模型研究通道內(nèi)的流動(dòng)沸騰.所研究物理問(wèn)題的示意圖如圖1所示，通道長(zhǎng)和寬分別為L(zhǎng)x=10 mm和Ly=3 mm，在通道的底部分布著若干個(gè)長(zhǎng)和高分別是L和H的矩形加熱器，加熱器之間的距離為D，加熱器個(gè)數(shù)記為N，加熱器表面接觸角為θ.初始時(shí)通道內(nèi)充滿了密度為ρl=5.426，溫度為Tsat=0.9Tc的飽和水，在計(jì)算域內(nèi)部沿x方向?qū)α黧w施加一個(gè)水平慣性力Fx，y方向上浮力為Fy=Fext.加熱器底部為恒溫邊界條件[21]，加熱器左右兩側(cè)絕熱，底部其他區(qū)域和計(jì)算域頂部均為無(wú)滑移絕熱邊界.計(jì)算域左右兩邊均采用周期性的邊界條件，本文中浮力與慣性力的關(guān)系為Fx=0.9Fy，重力加速度設(shè)置為g=（0,g），g=8 × 10-5.其他一些參數(shù)在模擬中固定如下：a=2/49，b=2/21，R=1.0，cv=6.0，ρl=5.426，ρg=0.811 3，加熱器的高度固定為H=2.另外，如無(wú)特別指明，本文中所有單位均為格子單位.

圖1 計(jì)算域示意圖Fig.1 The computation domain

2.2 模型驗(yàn)證

為了驗(yàn)證模型的正確性，我們計(jì)算了沸騰過(guò)程中不同重力加速度值下的單氣泡脫離直徑d.在初始時(shí)刻，計(jì)算域內(nèi)充滿了Tsat=0.9Tc的飽和水，并且在底壁的中心放置一個(gè)L×H的高溫加熱器.加熱器底部為恒溫邊界條件[21]，加熱器左右兩側(cè)絕熱.計(jì)算域底部其他區(qū)域和計(jì)算域頂部均為無(wú)滑移絕熱邊界，左右邊界都為周期性的邊界條件.對(duì)于該模型驗(yàn)證中，氣泡脫離直徑d與重力加速度相關(guān)[22]：

為了驗(yàn)證這種關(guān)系，在不同的重力加速度值下進(jìn)行了模擬，即模擬中g(shù)=2 × 10-5，4 × 10-5， 6 × 10-5， 8 ×10-5， 1 × 10-4，計(jì)算域?yàn)長(zhǎng)x×Ly=400 × 200.得到的氣泡脫離直徑見圖2，其中方形點(diǎn)表示氣泡離開的直徑，曲線為擬合結(jié)果.圖2中擬合曲線的指數(shù)為-0.507 3，與實(shí)驗(yàn)預(yù)測(cè)的-0.5吻合良好.

在進(jìn)行數(shù)值模擬之前，首先進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證.表1給出了網(wǎng)格大小分別為L(zhǎng)x×Ly=500 × 150，Lx×Ly=1 000 × 300，Lx×Ly=2 000 × 600三種情況下流動(dòng)沸騰過(guò)程中加熱器表面的空間-時(shí)間平均熱流密度Qave.其中加熱器表面接觸角θ=90°，其他參數(shù)設(shè)置和2.1小節(jié)相同.從表1中可以看出，網(wǎng)格數(shù)為L(zhǎng)x×Ly=500 × 150時(shí)的平均熱流密度低于其他兩種情況，此外，網(wǎng)格大小Lx×Ly=1 000 × 300和Lx×Ly=2 000 × 600兩種情況下的平均熱流密度變化相對(duì)誤差為4.054%.因此，考慮到計(jì)算效率和計(jì)算精度，本文網(wǎng)格大小選取為L(zhǎng)x×Ly=1 000 ×300，對(duì)應(yīng)的dy=0.01 mm，即最小壁面法向網(wǎng)格尺度為0.01 mm.

表1 網(wǎng)格無(wú)關(guān)性Table 1 Grid independence

3 數(shù)值結(jié)果與分析

3.1 加熱器間距的影響

本小節(jié)研究加熱器間距對(duì)流動(dòng)沸騰的影響，這里主要考慮了加熱器間距D=1 000，500，333，250四種情況下的流動(dòng)沸騰，模擬中加熱器長(zhǎng)度L=22，每個(gè)加熱器的表面接觸角均為90°.

圖3所示為中性表面情況下加熱器間距D=1 000時(shí)氣泡從生長(zhǎng)到離開的過(guò)程.如圖所示，氣泡在t=4 000時(shí)成核，此時(shí)氣泡形狀為一個(gè)半橢圓形，隨著沸騰過(guò)程的進(jìn)行，在t=10 000～21 000時(shí)刻可以很明顯看出氣泡由于慣性力作用向下游偏離.在t=10 000，15 000，17 000時(shí)刻，隨著時(shí)間的推移，氣泡在浮力和慣性力下和加熱器的接觸長(zhǎng)度越來(lái)越小，同時(shí)氣泡面積越來(lái)越大，最終氣泡在t=19 000時(shí)刻脫離加熱器表面.隨后在t=21 000時(shí)刻，加熱器表面產(chǎn)生成核氣泡，代表第二個(gè)沸騰循環(huán)的開始.

圖3 D=1 000時(shí)計(jì)算域的密度分布Fig.3 Density distributions for D=1 000

圖4為加熱器間距D=500，333，250時(shí)的氣泡演化過(guò)程.如圖所示，與圖3中的氣泡成核時(shí)間相同，D=500，333，250時(shí)得到的氣泡成核時(shí)間也是t=4 000時(shí)刻，表明加熱器間距不影響氣泡成核時(shí)間.然而不同的加熱器間距對(duì)氣泡脫離加熱器表面的時(shí)間有很大影響.當(dāng)D=500時(shí)，氣泡脫離時(shí)間在t=21 000時(shí)刻；當(dāng)D=333時(shí)，氣泡脫離時(shí)間在t=24 000時(shí)刻；當(dāng)D=250時(shí)，氣泡在t=29 000時(shí)刻脫離加熱器表面，即加熱器間距越小，氣泡越難以脫離加熱器表面，且氣泡的生長(zhǎng)周期越長(zhǎng).這是因?yàn)榧訜崞鏖g距較短時(shí)在加熱器表面產(chǎn)生的氣泡之間的相互作用力越大，該相互作用力影響周圍的流體運(yùn)動(dòng)和傳熱.

圖4 不同加熱器間距時(shí)計(jì)算域的密度分布：（a）D=500；（b）D=333；（c）D=250Fig.4 Density distributions with different heater spacings: (a)D=500; (b) D=333; (c) D=250

從圖4還可以看出不同加熱器間距下得到的氣泡大小也不相同，為了詳細(xì)說(shuō)明這一現(xiàn)象，圖5給出了不同加熱器間距下氣泡脫離前其面積A隨時(shí)間的變化，這里的氣泡面積為計(jì)算域內(nèi)密度小于 0 .5×(ρl+ρg)的流體格點(diǎn)所占的區(qū)域.需要說(shuō)明的是，由于本文主要研究加熱器對(duì)流動(dòng)沸騰的影響，因此這里主要統(tǒng)計(jì)氣泡脫離壁面之前的面積，而氣泡脫離加熱器之后面積變化受很多因素，例如流體溫度、慣性力大小等的影響，這里暫不討論.從圖中可以看出，對(duì)于不同的加熱器間距，氣泡面積幾乎從同一時(shí)刻開始增加，這和圖4所示的不同加熱器間距得到的氣泡成核時(shí)間相同的結(jié)論一致.從圖5可以看出，在氣泡成核之后，四種不同加熱器間距的氣泡面積都隨時(shí)間的推移而增加，氣泡脫離時(shí)加熱器間距越大氣泡面積越大，Wang等同樣也發(fā)現(xiàn)了加熱器間距對(duì)氣泡面積的大小有影響[12].由于較大的氣泡面積可以導(dǎo)致較大的浮力和慣性力，所以，加熱器間距越大，氣泡脫離時(shí)間越早.通常情況下，較大的氣泡可以從加熱器表面帶走更多的熱量.

圖5 不同加熱器間距下的氣泡面積變化Fig.5 The bubble areas with different heater spacings

為了進(jìn)一步說(shuō)明加熱器間距對(duì)流動(dòng)沸騰的影響，圖6給出了不同加熱器間距下加熱器表面平均熱流密度隨時(shí)間的變化，其中圖6（a）給出了加熱器表面空間平均熱流密度隨時(shí)間的變化，圖6（b）給出了加熱器表面空間-時(shí)間平均熱流密度.從圖6（a）中可以看出，氣泡成核之后，間距越大的加熱器表面空間平均熱流密度Q越大.圖中Q值下降到最低點(diǎn)的時(shí)刻對(duì)應(yīng)氣泡脫離時(shí)刻，分別對(duì)應(yīng)圖3中t=19 000時(shí)刻，圖4（a）中t=21 000時(shí)刻，圖4（b）中t=24 000時(shí)刻和圖4（c）中t=29 000時(shí)刻.可以看出，在相同的時(shí)間內(nèi)，沸騰周期隨著加熱器間距的增大而減小.圖6（b）為t=0～60 000時(shí)間內(nèi)加熱器表面的空間-時(shí)間平均熱流密度Qave，不難發(fā)現(xiàn)，隨著加熱器間距的增加，加熱器表面的空間-時(shí)間平均熱流密度值變大.當(dāng)加熱器間距從250增加到1 000時(shí)，其傳熱性能可以提高約12%.

圖6 不同加熱器間距下加熱器表面平均熱流密度的變化：（a）加熱器表面空間平均熱流密度隨時(shí)間的變化；（b）加熱器表面空間-時(shí)間平均熱流密度Fig.6 Heat fluxes on the heater surface with different heater spacings: (a) temporal variations of the spatial average heat flux on the heater surface;(b) the temporal and spatial average heat flux on the heater surface

3.2 加熱器長(zhǎng)度的影響

本小節(jié)研究了加熱器長(zhǎng)度對(duì)流動(dòng)沸騰過(guò)程中氣泡面積和表面平均熱流密度的影響.模擬中，計(jì)算域底部相鄰兩個(gè)加熱器中心的距離D=500，每個(gè)加熱器的表面接觸角均為90°.分別考慮加熱器的長(zhǎng)度L=16，18，20，22四種情況.

圖7所示為4種不同加熱器長(zhǎng)度的氣泡演化過(guò)程.如圖所示，以加熱器長(zhǎng)度L=16的情況為例，在t=12 000時(shí)刻加熱器表面產(chǎn)生成核氣泡，隨著沸騰過(guò)程進(jìn)行，氣泡面積逐漸變大，其受到的浮力和慣性力也隨著變大，在t=32 000時(shí)刻，氣泡在浮力和慣性力作用下向右上滑移，即氣泡與加熱器的接觸長(zhǎng)度逐漸減小，最終氣泡在t=38 000時(shí)刻完全脫離加熱器，隨后在t=45 000時(shí)刻第二個(gè)成核氣泡產(chǎn)生.觀察不同加熱器長(zhǎng)度下的四組氣泡演化過(guò)程發(fā)現(xiàn)，氣泡的成核時(shí)間隨著加熱器長(zhǎng)度的增加而提前，這是因?yàn)殡S著加熱器長(zhǎng)度的增加，更多的飽和液體被加熱，縮短了氣泡的成核時(shí)間.另一方面，越長(zhǎng)的加熱器氣泡脫離時(shí)間越早，這是因?yàn)榧訜崞髟介L(zhǎng)，表面的氣泡越大，意味著加熱器長(zhǎng)度影響氣泡脫離時(shí)間和氣泡生長(zhǎng)周期.

圖7 不同加熱器長(zhǎng)度時(shí)計(jì)算域的密度分布：（a）L=16；（b）L=18；（c）L=20；（d）L=22Fig.7 Density distributions with different heater lengths: (a) L=16; (b) L=18; (c) L=20; (d) L=22

為進(jìn)一步研究加熱器長(zhǎng)度對(duì)氣泡脫離時(shí)間和氣泡生長(zhǎng)速度的影響，圖8給出了不同加熱器長(zhǎng)度下，氣泡脫離加熱器前其面積A隨時(shí)間的變化.如圖所示，加熱器長(zhǎng)度越長(zhǎng)氣泡面積開始增加的時(shí)刻越靠前，雖然在不同加熱器長(zhǎng)度下，氣泡面積都隨時(shí)間而增加，但是氣泡面積的增大速度卻有所不同，在氣泡成核后的10 000時(shí)長(zhǎng)內(nèi)，氣泡面積增加的大小分別為401（L=22），332（L=20），257（L=18）和167（L=16）.顯然，當(dāng)L=22時(shí)氣泡面積增加的速度最快，L=20和L=18時(shí)氣泡面積增加的速度依次變慢，L=16時(shí)，氣泡面積增加速度最慢.圖中的最高點(diǎn)表示氣泡脫離加熱器表面時(shí)的氣泡面積，可以看出，加熱器越長(zhǎng),氣泡脫離時(shí)的面積越大且脫離時(shí)間越早.上述結(jié)果表明，在加熱器間距一定的情況下，加熱器越長(zhǎng)越有利于氣泡生長(zhǎng).

圖8 不同加熱器長(zhǎng)度下的氣泡面積變化Fig.8 The bubble areas with different heater lengths

圖9所示為不同加熱器長(zhǎng)度下加熱器表面空間平均熱流密度以及空間-時(shí)間平均熱流密度的變化.從圖9（a）可以看出，當(dāng)加熱器長(zhǎng)度不同時(shí)，其表面的空間平均熱流密度Q都是先降低.對(duì)于L=22，在t=4 000時(shí)刻產(chǎn)生成核氣泡，隨著氣泡的生長(zhǎng)，加熱器表面的熱流密度逐漸上升至最高點(diǎn)（t=16 000），在此之后，因?yàn)闅馀菖c加熱器的接觸面積逐漸減小，熱流密度Q逐漸下降，這種下降趨勢(shì)在氣泡與加熱器表面脫離時(shí)停止（t=21 000）.不同加熱器長(zhǎng)度下的熱流密度變化趨勢(shì)相同，只是在時(shí)間上有一些差別.比較不同加熱器長(zhǎng)度下的熱流密度曲線發(fā)現(xiàn)，氣泡脫離加熱器的時(shí)間隨著加熱器長(zhǎng)度的增加而提前.這是由于加熱器長(zhǎng)度較大時(shí)，產(chǎn)生的氣泡面積較大，其較大的浮力和慣性力促進(jìn)了氣泡從加熱器表面的脫離.圖9（b）給出了不同加熱器長(zhǎng)度下加熱器表面在前80 000時(shí)刻內(nèi)的空間-時(shí)間平均熱流密度Qave，如圖所示，加熱器越長(zhǎng)，Qave值越大，表面換熱性能越好，當(dāng)加熱器長(zhǎng)度從16增加到22時(shí)，其傳熱性能提高了13% .

圖9 不同加熱器長(zhǎng)度下加熱器表面平均熱流密度的變化：（a）加熱器表面空間平均熱流密度隨時(shí)間的變化；（b）加熱器表面空間-時(shí)間平均熱流密度Fig.9 Heat fluxes on the heater surface with different heater lengths: (a) temporal variations of the spatial average heat flux on the heater surface;(b) the temporal and spatial average heat flux on the heater surface

3.3 潤(rùn)濕性的影響

在沸騰過(guò)程中，受熱表面的潤(rùn)濕性也對(duì)傳熱性能有顯著影響.為了測(cè)試潤(rùn)濕性對(duì)流動(dòng)沸騰的影響，本小節(jié)模擬了親水性表面和疏水性表面兩種不同情況.對(duì)于每種情況，加熱器間距設(shè)置為D=500，加熱器長(zhǎng)度設(shè)置為L(zhǎng)=22.

圖10所示為親水性表面下（θ=77°）流動(dòng)沸騰過(guò)程的氣泡演化.氣泡在t=7 000時(shí)刻成核，t=16 000時(shí)刻后，氣泡在浮力和慣性力作用下與加熱器表面的接觸長(zhǎng)度越來(lái)越小，最終在t=21 000時(shí)刻脫離加熱器，下一個(gè)氣泡的成核時(shí)間在t=36 000時(shí)刻.圖11是一個(gè)疏水性表面（θ=120°）的氣泡演化過(guò)程，與親水表面相比，氣泡在疏水表面成核時(shí)間（t=2 000時(shí)刻）更早.另外，在t=26 000時(shí)刻，氣泡從疏水表面脫離后疏水表面有殘余氣泡，這是因?yàn)槭杷鎸?duì)氣泡的吸附力較大，這些殘余氣泡會(huì)作為下一個(gè)氣泡的核繼續(xù)生長(zhǎng)[23-24].

圖10 親水性表面（θ=77°）時(shí)計(jì)算域的密度分布Fig.10 Density distributions on the hydrophilic surface（θ=77°）

圖11 疏水性表面（θ=120°）時(shí)計(jì)算域的密度分布Fig.11 Density distributions on the hydrophobic surface（θ=120°）

不同潤(rùn)濕性下加熱器表面平均熱流密度的變化如圖12所示.可以看出，在氣泡成核初期（t=2 000 ～ 5 000），疏水性表面存在熱流密度下降期，這是因?yàn)槭杷砻嫔蠚馀莩珊怂俾矢?，由于加熱器尺寸的限制，較快的成核速率導(dǎo)致加熱器表面上形成一層來(lái)不及脫離的較厚氣泡，從而影響了加熱區(qū)域向外的熱量傳遞，導(dǎo)致熱流密度降低.在t=5 000時(shí)刻，慣性力的作用表現(xiàn)明顯，氣泡發(fā)生偏離，開始向下游運(yùn)動(dòng)，熱流密度在此時(shí)開始上升，在t=21 000時(shí)刻達(dá)到最大值，對(duì)應(yīng)圖11中t=21 000時(shí)刻，隨后熱流密度下降，t=26 000時(shí)熱流密度下降到最小，在這個(gè)時(shí)刻氣泡脫離，因?yàn)榇藭r(shí)加熱器表面仍然存在氣泡，所以表面平均熱流密度Q的值仍然較大.但是對(duì)于親水性表面，加熱器表面容易被冷的飽和液體潤(rùn)濕，氣泡脫離之后加熱器表面沒有殘余氣泡，熱流密度值明顯低于疏水性表面.同時(shí)，親水性表面氣泡脫離之后加熱器表面熱流密度在t=20 000 ～ 30 000時(shí)刻變化緩慢，因?yàn)榇藭r(shí)加熱器表面被冷液體潤(rùn)濕，沒有相變過(guò)程產(chǎn)生.比較兩種潤(rùn)濕性下的平均熱流密度，疏水性表面的平均熱流密度比親水性高26%.

圖12 不同潤(rùn)濕性表面下加熱器表面平均熱流密度隨時(shí)間的變化Fig.12 Temporal variations of the spatial average heat flux on the heater surface with different wettabilities

當(dāng)加熱器表面的潤(rùn)濕性不同時(shí)，流動(dòng)沸騰過(guò)程中產(chǎn)生的氣泡面積也存在差別，圖13所示為不同潤(rùn)濕性表面下氣泡面積的變化.很明顯，親水表面的成核時(shí)間晚于疏水表面，這與池沸騰的結(jié)果相同[25].對(duì)于氣泡面積大小而言，親水性表面產(chǎn)生的氣泡面積小于疏水性表面產(chǎn)生的氣泡面積.所以疏水性表面可以從加熱器表面帶走更多的熱量，優(yōu)于親水性表面的傳熱性能.

圖13 不同潤(rùn)濕性表面下氣泡面積的變化Fig.13 The bubble areas with different surface wettabilities

3.4 加熱器間距、長(zhǎng)度及潤(rùn)濕性的綜合影響

前文中已經(jīng)研究了加熱器的間距、長(zhǎng)度及潤(rùn)濕性對(duì)流動(dòng)沸騰的影響，本節(jié)將利用正交試驗(yàn)方案對(duì)這三種因素進(jìn)行綜合考慮， 表2給出了一個(gè)三水平三因素的正交試驗(yàn)表及其數(shù)據(jù)分析.表2的A行中三個(gè)值分別對(duì)應(yīng)的是各種因素下①水平試驗(yàn)指標(biāo)的平均值，即（8.463 × 10-3+ 1.043 × 10-2+ 1.264 × 10-2）/3=1.051 × 10-2；（8.463 × 10-3+ 1.032 × 10-2+ 1.327 × 10-2）/3=1.069 × 10-2；（8.463 × 10-3+ 9.856 × 10-3+ 9.087 × 10-3）/3=9.135 ×10-3.同理，B、C行中三個(gè)值分別對(duì)應(yīng)的是各種因素下②水平和③水平試驗(yàn)指標(biāo)的平均值.R行中三個(gè)值為各種因素下3個(gè)水平的極差值.

表2 正交試驗(yàn)表及數(shù)據(jù)分析Table 2 The orthogonal test table and data analysis

從表2可以看出，對(duì)于三種因素，都是③水平的情況下Qave值最大，①水平的情況下Qave值最小，因此可以得出結(jié)論，當(dāng)加熱器長(zhǎng)度L=22，加熱器間距D=1 000，加熱器表面接觸角θ=120°時(shí)，加熱器表面的熱流密度最高.另一方面，通過(guò)觀察R行的值發(fā)現(xiàn)，接觸角對(duì)Qave的影響最大，加熱器長(zhǎng)度次之，加熱器間距對(duì)Qave的影響最小.

4 結(jié) 論

利用偽勢(shì)LB兩相流模型，模擬了含多個(gè)矩形加熱器通道內(nèi)流動(dòng)沸騰過(guò)程中氣泡成核、生長(zhǎng)到脫離的過(guò)程.詳細(xì)分析了加熱器之間的距離、加熱器的長(zhǎng)度和加熱器表面的潤(rùn)濕性對(duì)通道內(nèi)流動(dòng)沸騰過(guò)程中氣泡動(dòng)力學(xué)行為和傳熱的影響.結(jié)論總結(jié)如下：

1）不同加熱器間距影響氣泡的脫離時(shí)間和生長(zhǎng)速率，加熱器間距越大，氣泡脫離時(shí)間越短，生長(zhǎng)速率越快，加熱器間距從250增加到1 000時(shí)其傳熱性能可以提高12%.

2）氣泡成核時(shí)間隨著加熱器長(zhǎng)度的增加而提前，增加加熱器長(zhǎng)度，更多的飽和液體被加熱，成核氣泡的生長(zhǎng)速度加快，使氣泡脫離加熱器表面的時(shí)間提前.當(dāng)加熱器長(zhǎng)度從16增加到22時(shí)，其傳熱性能提高了13%.

3）親水性表面的氣泡成核時(shí)間晚于疏水性表面的氣泡成核時(shí)間，與親水性表面相比，疏水性表面在氣泡脫離加熱器之后存在殘余氣泡.親水性表面的平均熱流密度和產(chǎn)生的氣泡面積小于疏水性表面.加熱器表面接觸角從77°變化到120°，其傳熱性能提高了26%.

4）通過(guò)正交試驗(yàn)方案發(fā)現(xiàn)，矩形加熱器表面的潤(rùn)濕性對(duì)流動(dòng)沸騰的傳熱性能影響最大，加熱器長(zhǎng)度的影響次之，加熱器間距的影響最小.

參考文獻(xiàn)（ References ）：

［1］KANG M.Effect of surface roughness on pool boiling heat transfer[J].International Journal of Heat and Mass Transfer, 2000, 43（22）: 4073-4085.

［2］GOEL P, NAYAK A K, KULKARNI P P, et al.Experimental study on bubble departure characteristics in subcooled nucleate pool boiling[J].International Journal of Multiphase Flow, 2017, 89: 163-176.

［3］HE H, LIU Y, LIU L, et al.Numerical simulation of bubble growth on and departure from the heated surface by an improved lattice Boltzmann model[J].Kerntechnik, 2018, 83: 186-192.

［4］GONG S, CHENG P.A lattice Boltzmann method for simulation of liquid-vapor phase-change heat transfer[J].International Journal of Heat and Mass Transfer, 2012, 55（17/18）: 4923-4927.

［5］WANG H, LOU Q, LI L.Mesoscale simulations of saturated flow boiling heat transfer in a horizontal microchannel[J].Numerical Heat Transfer(Part A):Applications, 2020, 78（4）: 107-124.

［6］DONG Z, XU J, JIANG F, et al.Numerical study of vapor bubble effect on flow and heat transfer in microchannel[J].International Journal of Thermal Sciences, 2012, 54: 22-32.

［7］NIE D, GUAN G.Study on boiling heat transfer in a shear flow through the lattice Boltzmann method[J].Physics of Fluids, 2012, 33: 043314.

［8］GONG S, CHENG P.Numerical investigation of saturated flow boiling in microchannels by the lattice Boltzmann method[J].Numerical Heat Transfer(Part A):Applications, 2014, 65: 644-661.

［9］TAN K, HU Y, HE Y.Effect of wettability on flow boiling heat transfer in a microtube[J].Case Studies in Thermal Engineering, 2021, 26: 101018.

［10］YIN X, TIAN Y, ZHOU D, et al.Numerical study of flow boiling in an intermediate-scale vertical tube under low heat flux[J].Applied Thermal Engineering, 2019, 153: 739-747.

［11］AZIZIFAR S, AMERI M, IMAN B.An experimental study of subcooled flow boiling of water in the horizontal and vertical direction of a metal-foam tube[J].Thermal Science and Engineering Progress, 2020, 20（1）: 100748.

［12］WANG J, CHENG Y, LI X B, et al.Experimental and LBM simulation study on the effect of bubbles merging on flow boiling[J].International Journal of Heat and Mass Transfer, 2019, 132: 1053-1061.

［13］VONTAS K, ANDREDAKI M, GEORGOULAS A, et al.The effect of surface wettability on flow boiling characteristics within microchannels[J].International Journal of Heat and Mass Transfer, 2021, 172: 121133.

［14］SUN T, GUI N, YANG X, et al.Numerical study of patterns and influencing factors on flow boiling in vertical tubes by thermal LBM simulation[J].International Communications in Heat and Mass Transfer, 2017, 86: 32-41.

［15］YADAV A, ROY S.Void fraction distribution for convective boiling flows in single and multiple heater rods assembly[J].Chemical Engineering Science, 2022, 247: 117063.

［16］BALTIS C, GELD C.Experimental investigation of the thermal interactions of nucleation sites in flow boiling[J].International Journal of Heat and Mass Transfer, 2014, 78（11）: 1208-1218.

［17］CHENG M, ZHANG B, LOU J.A hybrid LBM for flow with particles and drops[J].Computers & Fluids, 2017,155（20）: 62-67.

［18］SONG F, WANG W, LI J.A lattice Boltzmann method for particle-fluid two-phase flow[J].Chemical Engineering Science, 2013, 102: 442-450.

［19］LI Q, KANG Q, FRANCOIS M, et al.Lattice Boltzmann modeling of boiling heat transfer: the boiling curve and the effects of wettability[J].International Journal of Heat and Mass Transfer, 2015, 85: 787-796.

［20］YUAN P, SCHAEFER L.Equations of state in a lattice Boltzmann model[J].Physics of Fluids, 2006, 18:042101.

［21］ZHANG T, SHI B, GUO Z, et al.General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method[J].Physical Review E, 2012, 85: 016701.

［22］FRITZ W.Berechnung des maximal volumes von Dampfblasen[J].Physik Zeitschr, 1935, 36: 379-384.

［23］FENG Y, CHANG F, HU Z, et al.Investigation of pool boiling heat transfer on hydrophilic-hydrophobic mixed surface with micro-pillars using LBM[J].International Journal of Thermal Sciences, 2021, 163: 106814.

［24］GONG S, CHENG P.Lattice Boltzmann simulations for surface wettability effects in saturated pool boiling heat transfer[J].International Journal of Heat and Mass Transfer, 2015, 85: 635-646.

［25］GONG S, CHENG P.Lattice Boltzmann simulation of periodic bubble nucleation, growth and departure from a heated surface in pool boiling[J].International Journal of Heat and Mass Transfer, 2013, 64: 122-132.