基于混合熱格子Boltzmann方法有限熱源加熱氣泡動力學模擬

金志浩， 章港， 戰(zhàn)洪仁， 李帥， 王立鵬， 張志剛

(沈陽化工大學機械與動力工程學院， 沈陽 110142)

沸騰是一種高效的傳熱方式，近些年來，對于沸騰傳熱機理以及強化沸騰傳熱的研究，已經(jīng)成為國際上沸騰領域的研究熱點。此研究對于提升能源利用效率、減少工程投資有著重要意義。隨著電子元器件的集成化越來越高，普通的風冷已經(jīng)無法滿足電子元件的散熱需求，采用沸騰散熱已經(jīng)成為必然的趨勢，近些年，通過調控加熱表面的潤濕性來強化沸騰傳熱已成為研究熱點之一，很多學者通過實驗對親疏水表面沸騰傳熱的影響進行了研究。Takata等[1]通過在銅表面鍍聚四氟乙烯(polytetrafluoroethylene，PTEE)顆粒獲得疏水表面，發(fā)現(xiàn)疏水表面在加熱面溫度低于飽和溫度時，就能形成穩(wěn)定膜態(tài)沸騰。Theofanous等[2-3]通過紅外相機觀察壁面潤濕性對氣泡脫離直徑有較大的影響。Fan等[4]通過實驗發(fā)現(xiàn)隨著壁面潤濕性的增加，氣泡脫離頻率降低。Jo等[5]通過實驗研究了接觸角對換熱特性的影響。姜洪鵬等[6]用激光燒蝕方法在拋光后的銅上制備出超疏水/親水性的規(guī)則微陣列結構表面，發(fā)現(xiàn)疏水表面尾端氣泡容易匯聚，生長周期較長；而親水表面沒有發(fā)生明顯的氣泡匯聚行為，氣泡生長周期較短。但是沸騰傳熱過程的影響因素眾多，在實驗研究中難以區(qū)分，使得實驗結果的準確性受到質疑，而且通過實驗很難觀測氣泡生長、脫離過程的微觀流動細節(jié)，無法詳細描述出氣泡生長的內(nèi)在機理。

隨著計算機技術的飛速發(fā)展，流體體積法(volume of fluid，VOF)[7-9]、Level Set[10]、界面追蹤法(front tracking method，F(xiàn)TM)方法[11]及格子Boltzmann方法等數(shù)值模擬方法得到了廣泛的發(fā)展，并且已應用于研究沸騰和冷凝相變過程。但VOF、Level Set及FTM等數(shù)值模擬方法需要精確追蹤汽液界面，在氣泡和液滴的生長與脫離的過程中，追蹤汽液界面演化的方法并不能很好地滿足要求。此外，在應用這些數(shù)值方法模擬相變時，需要預先設置汽化/液化核。所以這些都意味著用傳統(tǒng)的計算流體力學方法模擬受限空間內(nèi)沸騰相變的傳熱過程仍然是一個挑戰(zhàn)。近些年，飛速發(fā)展起來的格子Boltzmann模型，在并行性和自動追蹤相界面方面具有固有的優(yōu)勢，已成為探索相變現(xiàn)象的強大的數(shù)值模擬工具。Boltzmann模型主要有自由能模型[12-14]、偽勢模型[15-17]、顏色模型[18-20]三種模型形式。偽勢模型[15]有著極高的穩(wěn)定性，可應用于大規(guī)模計算仿真、多相多組分仿真以及特定復雜流體情況和邊界條件時的仿真計算。但該模型存在密度比小，且模擬方法的穩(wěn)定性較差，嚴重限制了格子Boltzmann方法的實際應用，為解決此問題，Zhang等[21]首次將格子Boltzmann偽勢模型與有限差分法結合，并模擬大空間氣泡脫離的過程。Dong等[22]提出一種新的自由能模型，對熱格子Boltzmann模型進行了擴展與發(fā)展，形成混合格子Boltzmann模型。Gong等[23]在偽勢模型基礎上提出一個新的源相，提高數(shù)值穩(wěn)定性，并把P-R狀態(tài)方程代入格子Boltzmann方程，驗證格子Boltzmann結合P-R狀態(tài)方程穩(wěn)定性。Gong等[24]又提出一種引入流固張力的模型，使得壁面接觸角可以在0～180°調整，研究了池沸騰親疏水壁面氣泡生長傳熱過程。Gong等[25]基于P-R狀態(tài)方程的格子Boltzmann模型，用該模型模擬了大空間多氣泡生長過程。Li等[26]提出一種模擬多相流的混合熱格子流密度，研究了臨界熱流密度以及池沸騰氣泡的成核過程。汪鵬軍等[27]通過將密度分布函數(shù)與溫度分布函數(shù)進行耦合求解的方式建立了氣液相變模型，研究了氣泡脫離直徑與重力加速度的關系。綜上所述，運用格子Boltzmann方法對沸騰的研究，主要集中于大空間多氣泡核態(tài)沸騰。而多氣泡沸騰傳熱包括一系列氣泡生長、脫離、多氣泡的相互作用及表面粗糙度狀況等較多的影響因素，不利于從表面潤濕、過熱度等單一因素研究氣泡的生長機理及動力學特征，因此應用格子Boltzmann數(shù)值模擬研究討論單氣泡的生長過程，可有效避免多氣泡相互作用而產(chǎn)生非線性現(xiàn)象的影響，對多氣泡的運動特性的研究有重要意義。

現(xiàn)基于混合熱格子Boltzmann模型，建立單氣泡生長過程的數(shù)值模型，從微觀角度在不同潤濕性表面、不同有限熱源長度、不同過熱度等單一影響因素條件下，對單氣泡產(chǎn)生、長大、脫離過程的動力學特點進行模擬研究，并通過流場和溫度場對獲得的氣泡形態(tài)及氣泡生長過程的影響因素進行分析和比較，為研究強化表面沸騰傳熱及表面潤濕性的選擇提供理論依據(jù)。

1 基本模型

采用單組多相偽勢格子Boltzmann模型[18]耦合能量方程模擬氣液相變過程，其中單組多相偽勢格子Boltzmann模型密度演化方程具體表達式為

fi(x+eiδt，t+δt)-fi(x+t)

(1)

(2)

格子點的密度ρ與速度u計算公式為

(3)

所以格子點的密度ρ為單位格子上各個方向粒子分布函數(shù)之和。

在多相流的模擬中，體積力是十分重要的，諸多學者[12]用精確差分方法求解體積力項，這樣可以保證方程穩(wěn)定性，關系式為

(4)

為了使氣液兩相分離，引入流體間粒子作用力Fint，表達式為

(5)

式(5)中：β為提高計算穩(wěn)定的權重因子，取值1.16；G為相互作用常數(shù)，取值-1；x′為x的相鄰格子點坐標；ψ(x)為有效密度，與流體狀態(tài)方程有關，表達式為

(6)

式(6)中：p為流體壓力；g0為流體作用系數(shù);c0為格子系數(shù)。

通過引入Fs可獲得不同的三相接觸角，F(xiàn)s表達式為

(7)

式(7)中：sw為用來判斷格點是否是壁面的函數(shù)，對于固體壁面，可以取1，對于流體取0；Gw為流體與固體壁面的作用常數(shù)。

重力Fg表達式為

Fg=g(ρ-ρave)

(8)

式(8)中：g為重力加速度；ρave為計算域內(nèi)密度平均值。

模擬氣液相變現(xiàn)象，還需要在單組份多相偽勢模型的基礎上，耦合能量方程[28-29]構成相變模型。能量方程為

(9)

宏觀速度U計算公式為

(10)

壓力ρ可運用控制流體狀態(tài)的P-R狀態(tài)方程確定，即

(11)

最后對溫度場采用二階龍格庫塔法離散，獲得區(qū)域溫度場溫度分布，二階龍格庫塔方法表達式為

(12)

2 初始化模型與量綱分析

2.1 初始化模型與程序驗證

寬Lx=200，高Ly=250的封閉空間中充滿飽和溫度Tf=0.86Tc的氣液兩相流體，經(jīng)過網(wǎng)格獨立性驗證，計算域200×250可保證數(shù)值穩(wěn)定性。建立物理模型圖如圖1所示，為保障在容器底部產(chǎn)生單氣泡，在底部中心處設置加熱長度為3個格子長度的有限長度恒溫熱源P，加熱溫度Ts=1.3Tc。容器左右側壁面取周期邊界條件[31]，上下壁面取Zou-He邊界條件[32]。

調整重力加速度g的取值分別為-0.000 01、-0.000 02、-0.000 03、-0.000 04、-0.000 05，獲得不同氣重力加速度下泡脫離直徑曲線(圖2)， 橫坐標g為重力加速度，縱坐標D代表氣泡脫離直徑，其中單位為格子長度。氣泡脫離直徑與重力加速度g的擬合函數(shù)關系為D-g-0.5，符合Fritz提出的半經(jīng)驗公式，說明程序可靠。

圖1 物理模型Fig.1 Physical model

圖2 重力加速度g與氣泡脫離直徑D關系Fig.2 The relationship between gravitational acceleration g and detachment diameterD

2.2 無量綱轉換

格子Boltzmann模擬中用格子單位表示物理量，特征長度l0通常用來描述氣液相變過程，相對應的特征時間t0是由特征長度l0除以特征速度u0得到，計算公式[25]為

(13)

式(13)中：ρl為飽和液體密度，取值6.5；ρv為飽和氣相密度，取值0.38；σ為表面張力。因此，無量綱長度l*和無量綱時間t*分別定義為

(14)

因此無量綱直徑計算方法為

D*=D/l0

(15)

無量綱溫度計算方法為

T*=T/Tc

(16)

式(16)中：Tc流體的臨界溫度；T為流體的溫度；T*為無量綱溫度。

過熱度計算方法為

ΔT=Ts-Tf

(17)

式(17)中：Ts為恒溫熱源加熱溫度；Tf為飽和溫度。

3 計算結果與分析

首先初始時刻在容器底部壁面處放置密度為6.5的液滴，周圍充滿密度0.38的氣體，通過調整Gw數(shù)值，改變流體與固體壁面之間作用力，待程序運行穩(wěn)定后，分別獲得不同三相接觸角的壁面，表明不同親疏水表面，如圖3(a)、圖3(b)、圖3(c)、圖3(d)所示，黃色代表放置的液滴，藍色代表空間中的氣體，其中圖3(a)、圖3(b)為親水壁面，圖3(c)、圖3(d)為疏水壁面。

圖3 不同壁面三相接觸角Fig.3 Three phase contact angle on different wall surfaces

3.1 親疏水壁面氣泡生長過程

圖4和圖5分別為在接觸角為60°的親水表面上單個氣泡生成、成長和脫離流線圖和溫度場。初始時刻，在有限熱源附近形成一個過熱邊界層，如圖5(a)所示。隨著加熱時間的推移，由圖5(b)可以觀察到過熱邊界層內(nèi)產(chǎn)生一個氣泡核，氣泡核內(nèi)的蒸汽溫度升高，氣泡內(nèi)的溫差增大，導致氣液壓差加大，氣泡內(nèi)部流體的速度大于氣泡外的速度，依靠慣性作用推動氣泡核快速生長，這時親水表面的氣泡是以圓形狀態(tài)逐漸生長，而這個階段氣泡體積小，其對外部流場擾動可忽略，如圖4(a)、圖4(b)所示。隨后氣泡內(nèi)部慣性作用逐漸減小，氣泡的長大主要依靠過熱液膜提供熱量，液相轉化為氣相，加快周圍冷流體回流，在有限熱源兩側逐漸形成反方向旋轉的渦流，致使底部熱邊界層溫度降低，熱邊界層提供熱量不足以支撐氣泡快速生長，如圖5(c)所示。此時氣泡生長形態(tài)在表面張力、浮力流場剪切力作用下由圓形被拉成橢圓形，并在近壁面處形成頸部，圖4(c)、圖4(d)所示為氣泡生長階段。圖4(e)、圖4(f)為氣泡的脫離階段，可以觀察到氣泡脫離壁面，這是由于在貼壁處氣膜的氣化需要液體補充，冷流體的回流沖擊，強化在壁面處形成兩個反方向的渦流。在流場剪切力、浮力、壁面張力的拉扯下，氣泡頸部緊縮，壁面張力急劇減小，最終導致氣泡脫離。圖5(d)為氣泡脫離溫度場，氣泡脫離后，有限熱源附近又形成一層熱邊界層，氣泡脫離的尾流對熱邊界層帶來巨大的擾動。

圖4 接觸角60°親水壁面氣泡生長流線圖Fig.4 Bubble growth flow diagram on hydrophilic wall with contact angle of 60°

圖5 接觸角60°親水性壁面氣泡生長過程無量綱溫度場Fig.5 Dimensionless temperature field during bubble growth on hydrophilic wall with contact angle of 60°

圖6、圖7為接觸角為120°的疏水壁面氣泡生長脫離過程流場圖和溫度場。起初在有限熱源附近形成過熱邊界層中，如圖6(a)、圖7(a)所示，隨后在過熱邊界層內(nèi)形成氣泡核，由于氣泡內(nèi)蒸汽溫度較高，氣泡快速向外膨脹，在慣性力和表面張力作用下，氣泡沿著疏水壁面延伸生長，致使氣泡與疏水壁面三相線長度不斷增大，同時可觀察出疏水壁面氣泡兩端有著反方向旋轉的渦流，該渦流破壞邊界層強化換熱，加速氣泡生長，如圖6(b)和圖7(b)所示的氣泡成核階段。而在圖6(c)、圖6(d)氣泡生長階段，氣泡在疏水壁面的三相線長度縮小，主要原因是隨著氣泡生長慣性力逐漸減小，氣泡浮力不斷增大，氣泡垂直向上增長。氣泡向上生長兩側形成較大的渦流，該渦流促進氣泡在中間緊縮，形成頸部。隨著渦流的增強，氣泡內(nèi)部分布的溫度梯度逐漸減弱，熱量不足以支撐氣泡快速生長，氣泡生長速度減慢，如圖7(c)所示。在圖6(e)、圖6(f)氣泡脫離階段中，氣泡在中間斷裂脫離，這是由于冷流體回流較大，在兩側形成較強循環(huán)，在流場剪切力與氣液張力的拉扯下，在中間形成頸部，直至斷裂脫離。然而由于疏水壁面表面張力過強，依然有部分氣泡附著在疏水表面上作為下一個氣泡的汽化核心。如圖7(d)所示，氣泡脫離后帶走一部分熱量，氣泡脫離形成尾流，促使熱邊界層擾動。在氣泡生長、脫離過程中，溫度場與流場相互作用，促使氣泡成核、生長、脫離。

圖6 壁面接觸角120°疏水性壁面氣泡生長流線圖Fig.6 Hydrophobic wall bubble growth flow diagram with contact angle of 120°

圖7 壁面接觸角120°疏水性壁面氣泡生長 無量綱溫度場Fig.7 Dimensionless temperature field of hydrophobic wall bubble growth with contact angle of 120°

3.2 氣泡生長無量綱直徑與無量綱時間的關系

圖8表示在相同的過熱度下，親疏水壁面氣泡生長大小、脫離與時間關系圖，可以看出，親疏水壁面氣泡生長曲線先快速上升，隨后緩慢上升，最后保

圖8 親疏壁面氣泡D*與t*關系Fig.8 The relationship of bubble growth D*with t* on close and thin wall surfaces

持不變。疏水壁面氣泡脫離后，壁面有部分氣泡殘留，因此疏水壁面產(chǎn)生的氣泡大幅度減小。無論是親水壁面還是疏水壁面。分析可知，都存在快速生長階段、穩(wěn)定生長階段、氣泡脫離3個階段，其中快速生長階段氣泡快速膨脹，穩(wěn)定生長階段氣泡生長速度減慢，脫離階段氣泡不再生長，其內(nèi)在原因已在溫度場中論述。還可以觀察到疏水壁面較親水壁面先產(chǎn)生氣泡。疏水壁面曲線增長速率大于親水壁面，說明疏水壁面氣泡生長速度大于親水壁面氣泡生長速度。疏水壁面曲線增長速率大于親水壁面，說明疏水壁面氣泡生長速度大于親水壁面氣泡生長速度。疏水壁面氣泡脫離時刻有2個氣泡無量綱直徑，分別是100和61，說明由于疏水壁面具有親氣疏水的特性，氣泡脫離時被分成了2部分，一部分氣泡脫離壁面，其直徑為100，一部分氣泡殘留在壁面，直徑為61，殘留氣泡作為下一個氣泡生長的氣泡核。

3.3 不同有限熱源長度下壁面接觸角對氣泡生長脫離的影響

由圖9(a)可以看出，氣泡脫離直徑隨著壁面三相接觸角的增大而增大，即由親水表面轉變成疏水表面，氣泡脫離時間也逐漸加長。主要原因是三相接觸角的增大，導致壁面張力增大，氣泡脫離壁面需要更大的浮力克服壁面張力，因此在壁面需孕育更大直徑氣泡，生長時間變長，氣泡脫離直徑亦增加。故可通過減小壁面三相接觸角，即改變表面親疏性來使氣泡脫離直徑減小，縮短氣泡從生成到脫離的時間，從而可以達到壁面強化傳熱。親水壁面氣泡脫離直徑隨著有限熱源加熱長度增加而增加，疏水壁面氣泡脫離直徑與有限熱源長度無關。親水壁面氣泡脫離時間隨著有限熱源長度增加而減小，疏水壁面氣泡脫離時間隨著有限熱源加熱長度增加而小幅度減小。由圖9(b)可知，親水壁面有限熱源長度增加，單位時間提供的熱量增加，因此脫離直徑增加，脫離時間減小。疏水壁面上覆蓋一層氣膜，有限熱源長度增加對疏水表面氣泡脫離直徑影響不大，脫離時間小幅度減小。由此可以說明，通過調整有限熱源的大小，可以改善壁面的換熱效果。

圖9 不同有限熱源長度下壁面接觸角對氣泡脫離的影響Fig.9 Influence of wall contact angle on bubble detachment under different finite heat source lengths

3.4 有限熱源過熱度對氣泡脫離影響

圖10為有限熱源過熱度對氣泡生長脫離的影響。圖10(a)為過熱度與氣泡脫離直徑關系圖，親水壁面氣泡脫離直徑隨著有限熱源過熱度升高而增大，表明氣泡脫離直徑與過熱度成正比關系。由圖10(a)中曲線可知，有限熱源過熱度升高，疏水壁面氣泡脫離直徑基本不變，這是由于疏水壁面氣泡生長，氣泡在中間斷裂脫離主要受浮力、流場剪切力、氣泡張力影響下，斷裂脫離，故疏水表面過熱度對氣泡脫離直徑影響不大，如圖6(a)～圖6(f)所示。

圖10 有限熱源過熱度對氣泡生長脫離影響Fig.10 Influence of superheat of heat source on bubble growth detachment

圖10(b)為過熱度對氣泡脫離時間影響的曲線圖，親水壁面有限熱源過熱度的增加，導致氣泡增長速度加快，因而氣泡脫離時間減小。對于疏水壁面，氣泡脫離時間隨著過熱度增加而大幅度減小，這是由于疏水表面在氣泡形成初期先形成了一層氣膜，但是隨著過熱度的增大，有限熱源表面氣膜傳熱能力增強，氣泡的脫離時間將持續(xù)大幅度減小。

4 結論

在混和熱Boltzmann方法基礎上模擬了重力加速度對氣泡脫離直徑影響，分析了親水壁面氣泡的生長機理，得出以下結論。

(1)親疏水壁面氣泡生長機理不相同，親水壁面氣泡內(nèi)部溫度高，在慣性力的推動下長大，隨后慣性力減弱，浮力增大，流體回流，氣泡在浮力、流場剪切力，表面張力的作用下形成頸部直至脫離，壁面無氣泡殘留。疏水壁面氣泡初始在慣性力與壁面張力下，沿壁面生長，隨后慣性力減弱，浮力增強，氣泡在流場剪切力與氣液張力的拉扯下，在中間斷裂脫離，壁面有大量氣泡殘留。

(2)在相同過熱度下，親疏水壁面都有快速生長、穩(wěn)定生長、氣泡脫離3個階段，其中快速生長階段氣泡生長最快，穩(wěn)定生長階段氣泡生長速度降低，脫離階段氣泡大小基本不變。疏水壁面較親水壁面先產(chǎn)生氣泡，疏水壁面氣泡生長速度高于親水壁面。

(3)氣泡脫離直徑隨著壁面三相接觸角的增大而增大，脫離時間隨三相接觸角增大而增大。親水壁面氣泡脫離直徑隨著有限熱源加熱長度增加而增加，疏水壁面氣泡脫離直徑幾乎不隨有限熱源長度增加而改變。親水壁面氣泡脫離時間隨著有限熱源長度增加而減小，疏水壁面氣泡脫離時間隨著有限熱源加熱長度增加而小幅度減小。

(4)揭示了親疏水壁面氣泡脫離直徑與過熱度的關系，親水壁面氣泡脫離直徑與過熱度成正比，過熱度的增加會使氣泡生長速度加快，氣泡脫離時間減小。疏水壁面過熱度對氣泡脫離直徑影響不大，加熱溫度升高，傳熱推動力增強，氣泡脫離直徑基本不變，所以氣泡脫離時間大幅度減小，疏水壁面氣泡脫離時間明顯大于親水壁面，表明表面潤濕特性對氣泡脫離起主導作用，疏水性表面會增大汽泡的脫離阻力。