偽勢格子Boltzmann方法的汽化熱阱效應模擬研究

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

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

沸騰傳熱具有低溫差、高熱流密度的特性，被廣泛應用于冶金、芯片制造等領域.強化傳熱有著巨大的經濟效益，因此對沸騰傳熱機理的研究已經成為國際上研究熱點.近些年國內外專家先后提出微液膜蒸發(fā)模型[1-2]、瞬態(tài)導熱模型[3]、三相線模型[4]，但是到目前為止，對于沸騰強化傳熱機理沒有統(tǒng)一解釋.

1999年，沈自求在分析異氣沖擊加熱壁面[5]時，提出汽化界面熱阱效應，即異種氣體沖擊加熱壁面,使氣液界面上的液相快速汽化，吸取大量的汽化潛熱，導致氣液界面溫度降低，并且壁面溫度也隨著汽化熱阱效應的形成而降低.該效應很好地解釋了沸騰傳熱增強原理.沸騰時氣泡核與加熱壁面之間存在液體微層，汽化熱阱效應形成，其極大地加強了液體微層的導熱，為我們研究沸騰傳熱機理提供新的思路.通過實驗很難觀測和測量到汽化熱阱現(xiàn)象的微觀過程，數(shù)值模擬可以有效記錄與檢測流場的速度和溫度變化.

格子Boltzmann技術是一種介觀模擬方法，該方法在處理非線性問題、復雜界面問題、相變問題有著先天優(yōu)勢，因此廣受科研工作者青睞.1993年，Shan等首次提出偽勢模型，即把作用力代入偽速度勢中，該方法可以促使兩相自發(fā)分離，使偽勢模型有著極高的穩(wěn)定性[6].2011年，Zhang用偽勢格子Boltzmann模型與熱力學方程相耦合，形成熱格子Boltzmann模型，首次模擬了相變過程，但是控制流體狀態(tài)方程選用范德瓦爾方程，導致模擬結果數(shù)值穩(wěn)定性差[7].2018年，Gong提出把P-R狀態(tài)方程代入能量方程，通過四階龍格庫塔格式求解能量方程，得出演化后的溫度場，該方法數(shù)值穩(wěn)定性強，模擬結果的熱力學性質更接近真實流體[8].

為了避免多氣泡相互作用而產生的影響，本文選用單氣泡模型來解釋沸騰傳熱機理.運用偽勢格子Boltzmann方法，對單氣泡成核、生長、脫離、重新生成，以及周圍液體微層的流體細觀流動與傳熱特性進行了模擬研究，首次結合界面汽化熱阱效應對單氣泡生長及其傳熱規(guī)律進行分析，結合熱流密度變化曲線，討論其強化傳熱機理，為研究相變、氣泡動力學等行為對沸騰傳熱特性的影響提供了依據.

1 偽勢格子Boltzmann模型

單組分偽勢格子Boltzmann模型中流體粒子分布函數(shù)的演化方程為[9]

fi(x+eiδt,t+δt)-fi(x,t)=

(1)

(2)

(3)

密度ρ與粒子速度u計算公式如下:

(4)

宏觀速度U為

(5)

式中:Δu=Fδt/ρ為速度變量;F為流體受到的合力，其由流體粒子間兩相作用力Fint、流體粒子與固體壁面作用力Fs、重力Fg組成，合力F=Fint+Fs+Fg.

引入流體間粒子作用力Fint是為了促使氣液兩相分離.

(6)

式中[10]：β為權重因子，取值1.16，作用是提高數(shù)值穩(wěn)定性；G為氣液兩相相互作用力常數(shù)，通常取值-1；x′為x的相鄰格子點坐標；偽勢函數(shù)ψ(x)的選取與流體的狀態(tài)方程有關.

(7)

其中：p為流體壓力；g為流體作用系數(shù)；c0為格子結構確定系數(shù).

(8)

其中：sw用來判斷格點是否是壁面的函數(shù)，對于固體壁面取1，對于流體壁面取0；Gw是流體與固體壁面的作用常數(shù)，調節(jié)Gw可以獲得不同接觸角壁面.

重力Fg表達式如下:

Fg=g(ρ-ρave).

(9)

式中g為重力加速度.

以上為格子Boltzmann單組分偽勢模型，模擬相變過程采用以下溫度場的能量方程[11-12].

(10)

式中:λ為導熱系數(shù)；cv為定壓比熱；PEOS為流體壓力;?為哈密頓算子.

導熱系數(shù)為

1.3.1 面積 種公羊面積需求為每只4m2以上，妊娠母羊面積需求為每只1m2以上，空懷母羊面積需求為每只0.8m2以上，育成羊面積需求為每只0.6m2以上。

λ=0.028ρ.

(11)

壓力PEOS采用控制流體狀態(tài)的P-R狀態(tài)方程來得到.

(12)

然后對溫度場應用二階龍格庫塔法進行離散.

Tt+δt=Tt+(h1+2h2).

(13)

其中:Tt為當前時間步長t下的溫度場;Tt+δt為下一個時間步長的溫度場.

2 模擬與結果

2.1 計算域與物理模型

計算區(qū)域如圖1所示.計算空間內充滿飽和溫度為Tf=0.86Tc的氣液兩相流體，飽和液體密度ρl取值為6.5，飽和氣體密度ρv取值為0.38，液相與氣相導熱系數(shù)比為17，無量綱重力加速度g取為-0.000 05.經過網格獨立性驗證，計算域200×250可保證數(shù)值的穩(wěn)定性.設置加熱長度為3個格子長度的恒溫熱源，其位置在下表面中心處，熱源溫度Ts=1.3Tc，熱源附近壁面設為絕熱邊界條件[13-14]，容器左右側壁面取周期邊界條件[15-16]，上下固體壁面取Zou-He邊界條件[17].

圖1 物理模型

為驗證程序正確性，對不同重力加速度下氣泡脫離進行模擬，此前有數(shù)值模擬表示氣泡脫離直徑D與g-0.5成正比.在程序驗證中，取無量綱重力加速度分別為-0.000 05、-0.000 04、-0.000 03、-0.000 02、-0.000 01，得到氣泡脫離直徑與無量綱重力加速度的關系，如圖2所示，橫坐標g為無量綱重力加速度，縱坐標D為氣泡脫離直徑.由圖2計算可得，氣泡脫離直徑與無量綱重力加速度的擬合函數(shù)關系為D～g-0.5，符合Fritz理論.

圖2 無量綱重力加速度g與脫離直徑D的關系曲線

模擬中引入了特征長度l0、特征速度u0和特征時間t0.

(14)

(15)

t0=l0/u0.

(16)

其中:l0為毛細長度，是描述氣液兩相流的典型特征尺度，表示表面張力與浮力之比，根據拉普拉斯毛細方程獲得對應飽和溫度下的表面張力σ.

無量綱長度和無量綱時間可定義為：

l*=l/l0,

(17)

t*=t/t0.

(18)

式中:l為格子長度；t為格子時間.

2.2 親水壁面單氣泡生長過程

通過改變公式(8)中Gw的數(shù)值，得到接觸角為60°的親水壁面，以其模擬單個氣泡成核、生長、脫離過程，如圖3所示.圖3(a)中液體開始相變形成氣泡核，但是由于親水壁面的親水疏氣性,相變的位置并不是緊貼加熱面,而是在加熱面正上方的飽和液體中進行，因此加熱面與氣泡核之間存在一層薄薄的液體微層.從圖3(a)至圖3(c)，氣泡不斷生長，期間主要依靠液層微層的相變?yōu)闅馀蓍L大提供氣相，至圖3(b)時，氣泡下方的氣液界面開始與壁面接觸，至圖3(c)時液體微層全部汽化，隨后有流體補充至加熱面，氣泡繼續(xù)生長[圖3(d)].在圖3(c)至圖3(d)過程中，氣泡在壁面張力的束縛下由圓形逐漸被拉伸成橢圓形，隨后完全脫離.圖3(e)為氣泡脫離后的形態(tài)圖.

圖3 親水壁面單氣泡生長、脫離形態(tài)

圖4模擬了親水壁面單個氣泡成核、生長、脫離過程的溫度場.在氣泡成核階段，隨著熱源不斷地加熱，熱量首先通過導熱傳遞給附近的流體，由圖4(a)可以看出，此時熱源附近液體微層的溫度大于周圍飽和液體的溫度，處于過熱狀態(tài).對比圖4中的(a)、(b)、(c)這3個階段可以看出：近壁處的流體溫度呈下降趨勢，這是氣泡的生長階段，因為這個狀態(tài)下的氣泡生長主要依靠液體微層的相變，而相變的同時會帶走大量的潛熱，所以導致近壁處的流體溫度下降.

圖4 親水壁面單氣泡生長、脫離溫度場

從圖4(c)、圖4(d)可以看出：熱源附近的流體溫度逐漸升高，這是由于氣泡脫離加熱面，液體回流至加熱壁面，由于液相的熱阻遠小于氣相，所以，熱量迅速由熱源傳遞至其附近流體中，使得溫度升高.

2.3 疏水壁面單氣泡的生長過程

通過改變公式(8)中Gw的數(shù)值，得到接觸角為120°的疏水壁面，以其模擬單個氣泡成核、生長、脫離過程，如圖5所示.從如圖5(a)中可以看出：初始階段，由于熱源對流體不斷加熱，此時靠近熱源的液體開始發(fā)生相變，由于疏水壁面的親氣疏水性，相變是在緊貼壁面處進行.在圖5(a)至圖5(b)過程中，因為疏水壁面表面張力作用，氣泡開始橫向生長，三相接觸線逐漸變長，此時主要依靠氣泡周圍的流體相變?yōu)闅馀萆L提供氣相.在圖5(b)至圖5(c)這個階段，隨著氣泡的生長其體積不斷增大，導致其所受的浮力不斷增大，氣泡開始垂直向上移動，三相接觸線開始變短.從圖5(c)至圖5(d)過程中可以看出：在壁面張力與氣泡浮力的拉扯下，氣泡在中間形成頸部并開始脫離.從圖5(d)至圖5(e)，氣泡完成脫離，同時在壁面殘留一部分氣泡，作為下一個氣泡生成的氣泡胚核. 對比親水壁面，疏水壁面氣泡生長周期遠大于親水壁面氣泡生長周期.

圖5 疏水壁面氣泡生長形態(tài)

圖6模擬了疏水壁面單個氣泡成核、生長、脫離過程的溫度場.從圖6(a)中可以看出,此時熱源兩側近壁處的溫度降低.這是由于疏水壁面的相變緊貼壁面進行，在壁面形成三相接觸線，而此時氣泡的長大主要依靠三相線處的液體相變，并且這時的三相線剛剛形成，距離熱源比較近，所以,導致熱源兩側近壁處的溫度降低.從圖6(b)中可以看出,此時熱源附近的流體溫度逐漸升高.這是由于隨著氣泡的生長,三相線逐漸遠離熱源，此時熱源附近被氣體覆蓋,不會發(fā)生相變，所以，隨著加熱過程的進行，熱源附近的流體溫度逐漸升高.圖6(c)至圖6(d)，隨著氣泡的上移，三相線開始縮短逐漸靠近熱源，導致熱源附近的流體溫度再次降低.

圖6 疏水壁面氣泡溫度場

從模擬的結果中發(fā)現(xiàn)：加熱親水表面至其相變所需要時間為t*=1.8，同樣加熱條件下，加熱疏水表面至其相變所需要時間為t*=1.2.由此可見，相同加熱溫度下，相較于親水壁面，疏水壁面率先汽化.這是由于疏水表面具有親氣疏水特性，導致流體不會在加熱面附近形成過熱液體微層，而是直接相變；親水表面由于其親水疏氣特性，熱源處的流體緊貼壁面，在加熱初期，相比于疏水壁面這部分流體很難發(fā)生相變，從而形成過熱液體層，隨著加熱過程的進行，這部分流體積累能量逐漸增大，導致其在內部發(fā)生相變.Gao[18]在實驗中，用激光干涉法測得疏水壁面相變不存在微液層，親水壁面相變存在微液層.模擬結果與實驗結果一致.

2.4 親水壁面氣泡生長時氣泡氣液界面對熱源附近流體溫度的影響

圖7為熱源與氣泡氣液界面溫差ΔT和時間t*關系圖.A—B段為氣泡成核階段，由于過熱液體微層的相變帶走大量潛熱,導致氣液界面處的流體溫度下降，使得熱源與氣液界面的溫差逐漸增大.B—C段為氣泡縱向生長階段(此時主要依靠液體微層的相變維持氣泡生長)，隨著氣泡的生長,氣液界面逐漸與壁面接觸,導致其換熱面積逐漸減小，所以,這一階段熱源與汽化界面的溫差增長斜率小于A—B段.C—D段為氣泡在加熱面橫向生長階段(此時主要依靠三相線處的液體相變維持氣泡生長)，從C點開始,近壁處的氣泡氣液界面全部消失,此時熱源附近的流體全為氣體,即C—D段不會發(fā)生相變,所以,熱源會直接加熱其附近的氣體,使其溫度升高,從而導致熱源與附近流體的溫差減小.

圖7 熱源與氣液界面溫差ΔT與無量綱時間t*的變化關系

2.5 親水壁面氣泡生長時氣泡氣液界面對熱流密度曲線的影響

為了研究氣泡氣液界面對沸騰傳熱的影響，需要計算熱源與其附近流體的熱流密度.氣泡氣液界面主要依靠熱源的導熱傳遞熱量，所以應用傅里葉導熱定律進行計算，具體公式如下：

q(x)=-λ(?T/?y)|y=0.

圖8為近壁處熱流密度q與時間t*關系圖，A、B、C、D與圖7中的4個點相對應.由圖7可知A—B流體的溫度升高，導致傳熱溫差降低，因此，熱流密度曲線在A—B段緩慢下降.B—C段的熱流密度包括兩部分，即導熱和熱阱效應，熱阱效應起主導作用，由圖7可知在B—C段溫差逐漸升高，從而導致液體轉化為氣體的量也在增加，而隨后由于氣泡氣液界面的高度降低，液體逐漸減少，導致相變釋放的熱量也隨之減少，所以，熱流密度在B—C段出現(xiàn)最大極值，即氣泡氣液界面緊貼壁面時，熱流密度大幅度上升到C點.C—D段進入脫離階段，氣泡逐漸脫離壁面，氣泡氣液界面高度逐漸升高，源源不斷冷流體補充到壁面，相變減小，因此，C—D段熱流密度大幅度減小.

圖8 熱流密度q與無量綱時間t*的關系

3 結 論

使用混合熱格子Boltzmann方法，對親疏水壁面單氣泡成核、生長、脫離進行模擬，結合氣液界面汽化熱阱效應對傳熱特性的影響進行分析，得到以下結論：

(1) 相同加熱條件下，相比于親水壁面，疏水壁面率先汽化.這是由于疏水壁面具有親氣疏水性，致使其不會形成過熱液體微層，直接發(fā)生相變；而親水壁面由于親水疏氣特性，會先形成過熱液體微層，而后發(fā)生相變.

(2) 疏水壁面氣泡平鋪壁面生長，熱源附近不存在液體微層，所以，疏水壁面不存在汽化熱阱現(xiàn)象；親水壁面氣泡在壁面生長，存在液體微層，因此，親水壁面存在汽化熱阱現(xiàn)象.

(3) 親水壁面氣泡生長脫離周期較疏水壁面短.親水壁面氣泡生長不同階段，熱源與氣液界面間熱流密度不相同：初始階段，熱源與附近流體溫差逐漸降低，因此，熱流密度逐漸減小；氣泡生長階段，由于液相相變產生汽化熱阱效應，帶走大量汽化潛熱，大幅度提升熱流密度，強化傳熱;脫離階段，氣泡逐漸脫離壁面，熱流密度大幅度減小.