單氣泡上升對近壁面?zhèn)鳠嵊绊懙臄?shù)值研究

張文峰，龐明軍

(常州大學(xué)機(jī)械與軌道交通學(xué)院，江蘇 常州 213164)

泡狀流動(或氣泡在液相中的運(yùn)動)廣泛存在于在自然界和工業(yè)過程中，如巖漿流動、金屬冶煉、石油輸送等過程.氣泡的存在有助于液相內(nèi)(或固壁)熱量的傳輸和擴(kuò)散，均勻溫度場的分布，起到強(qiáng)化傳熱的作用[1-3]，從而實現(xiàn)節(jié)能減排的目的.文獻(xiàn)研究表明，當(dāng)有小氣泡出現(xiàn)在加熱壁面時，局部傳熱系數(shù)可提高2倍～3倍[4].因此，為了理解氣泡對不同泡狀設(shè)施內(nèi)固壁和液相間的傳熱影響規(guī)律和影響機(jī)理，國外學(xué)者相繼開展了相關(guān)研究.

理論上，Deckwer[1]推導(dǎo)出鼓泡塔反應(yīng)器的理論傳熱方程，揭示了傳熱系數(shù)對氣泡的依賴關(guān)系.Khabeev[5]提出了一種描述氣泡與液相間的傳熱方法，并導(dǎo)出氣泡邊界處的熱流表達(dá)式.數(shù)值上，Larimi和Ramiar[6]用流體體積法(VOF)研究了管道內(nèi)靜止和非靜止液體內(nèi)單氣泡上升對努塞爾數(shù)的影響.他們發(fā)現(xiàn)，當(dāng)液相的雷諾數(shù)較高時，氣泡上升對努塞爾數(shù)沒有明顯的影響；而當(dāng)液相的雷諾數(shù)較低時，氣泡上升對努塞爾數(shù)有顯著的增強(qiáng)作用.Deen和Kuipers[7]使用VOF方法研究了氣泡在熱壁附近上升對固壁與液相傳熱系數(shù)的影響，發(fā)現(xiàn)由于氣泡對流體的擾動，特別是在熱壁附近有氣泡合并時，壁面與液體間的傳熱系數(shù)顯著提高.Bhuvankar和Dabiri[8]使用前沿跟蹤法(FT)計算了單氣泡在剪切流中上升時對固壁面和液相間傳熱的影響；結(jié)果表明，在氣泡的下游和上游區(qū)域，傳熱表現(xiàn)出強(qiáng)化和減弱的現(xiàn)象.實驗上，Donoghue和Delaure[9]研究了單個氣泡撞擊熱水平壁面的傳熱現(xiàn)象，發(fā)現(xiàn)與單純的自然對流相比，反彈氣泡顯著增強(qiáng)了水平壁面的傳熱系數(shù).在隨后的研究中，他們還發(fā)現(xiàn)，單個氣泡撞擊加熱水平壁面時，在氣泡影響區(qū)域和其尾流區(qū)域內(nèi)，對流傳熱的變化在空間上是對稱的，并且氣泡附近的熱通量呈現(xiàn)出大幅、快速波動的現(xiàn)象[10].Qiu和Dhir[11]研究了氣泡沿傾斜加熱壁面上滑時的流動模式和傳熱現(xiàn)象，發(fā)現(xiàn)上滑氣泡尾流區(qū)的傳熱得到了顯著強(qiáng)化.

雖然國內(nèi)外學(xué)者對此類課題已經(jīng)展開了一定的研究，但是氣泡傳熱強(qiáng)化的機(jī)理尚不明確.因此，探究氣泡上升對傳熱的影響規(guī)律和強(qiáng)化機(jī)理對于泡狀設(shè)施的設(shè)計和控制具有重要意義.本文利用VOF法，詳細(xì)研究了表面張力、重力水平和氣泡距固壁距離對單氣泡上升固壁與液相間傳熱的影響規(guī)律，以獲得一些潛在的機(jī)理解釋.

1 物理模型和工況設(shè)置

為了研究氣泡上升對固壁與液相間傳熱的影響，幾何模型如圖1所示.為了簡化計算，將計算模型簡化為二維平面模型.為使氣泡的運(yùn)動在計算域內(nèi)能達(dá)到穩(wěn)定運(yùn)動狀態(tài)，且避免左壁AC對氣泡變形和運(yùn)動的影響，故計算區(qū)域的尺寸設(shè)置為，高H=50d、寬L=(40/3)d.氣泡圓心距右壁BD的距離為l=1.5d和2.5d，距CD為h=(17/6)d.d為氣泡的初始直徑，d=6 mm.

圖1 計算區(qū)域

計算時，對于流動而言，將AB、AC、CD和BD均設(shè)置為無滑移壁面，重力沿著y軸負(fù)方向；對于傳熱而言，分別在BD和AC施加393 K和293 K的壁面溫度，AB、CD設(shè)置為絕熱邊界條件.另外初始狀態(tài)下的氣相和液相都是靜止的，溫度均為293 K.

氣泡在靜止液體內(nèi)上升時，主要受重力、黏性力和表面張力的影響，這三種力可由無量綱參數(shù)E?tv?s(Eo)數(shù)和Gallilei(Ga)數(shù)表示.Eo數(shù)是重力和表面張力的比值，Ga數(shù)是重力和粘性力的比值，其定義式見公式(1)和公式(2).通過改變表面張力、重力水平和氣泡距固壁的距離，分別控制Eo數(shù)、Ga數(shù)和l，來改變氣泡運(yùn)動狀態(tài)，從而研究氣泡對壁面對流傳熱的影響.因此，設(shè)計了三組工況：(1)Eo=1、5、20和80時，Ga=14.5(1g)，l=1.5d；(2)Ga=14.5(1g)、12.6(0.75g)和10.3(0.5g)時，Eo=5、20和80，l=1.5d；(3)l=1.5d和2.5d時，Eo=5和80，Ga=14.5(1g).

(1)

(2)

公式中：ρl為液相密度，kg/m3；g為重力加速度，m/s2；d為氣泡直徑，m；μg為液相黏度，Pa·s；σ為表面張力，N/m.本文計算工況的物性參數(shù)，如表1所示.

表1 氣液物性參數(shù)

2 控制方程和數(shù)值方法

2.1 控制方程

計算時候假設(shè)液相和氣相都是不可壓縮流體，且流動近似為層流，液相和氣相的控制方程如下.

連續(xù)性方程：

?·u=0，

(3)

動量方程：

(4)

公式中：u為速度矢量，m/s；ρ為流體局部平均密度，kg/m3；μ為流體局部平均黏度，Pa·s；p為壓力，Pa；Fσ為由表面張力引起的體積力，N；?為哈密頓算子.

流體的局部平均密度和粘度分別按下式進(jìn)行計算：

μ(f)=μ1(f)+μg(1-f)，

(5)

ρ(f)=ρ1(f)+ρg(1-f)，

(6)

公式中：下標(biāo)g和l分別為氣相和液相，后文與此相同.

為了得到溫度場，忽略黏性耗散的影響，能量守恒方程可寫為

(7)

體積熱容ρCp和局部平均熱導(dǎo)率λ分別用下面的式子計算：

ρCp(f)=ρlCp，l(f)+(1-f)ρgCp，g，

(8)

λ(f)=λl(f)+λg(1-f)

.

(9)

2.2 VOF方法

為了精確捕捉氣液兩相的界面變化，故采用計算精度高計算速度快的VOF法.VOF法的相函數(shù)輸送方程為

(10)

相函數(shù)f可定義為

(11)

得知相函數(shù)，可利用公式(12)計算氣泡的局部體積分?jǐn)?shù)：

(12)

公式中：s為網(wǎng)格編號；AS編號為s網(wǎng)格單元的面積.界面的重構(gòu)利用分段線性法(PLIC)，該方法使用直線段近似表示單個網(wǎng)格內(nèi)氣液兩相的界面.界面法向量通過求解網(wǎng)格單元中氣泡的局部體積分?jǐn)?shù)梯度來獲得，界面的法向量n表示為

(13)

2.3 連續(xù)表面張力模型

Brackbill等[12]所提出的連續(xù)表面張力模型被用來計算表面張力.在計算過程中，表面張力作為源項體積力Fσ出現(xiàn)在動量方程的右邊，即

Fσ=σκδn，

(14)

公式中：δ為狄拉克分布函數(shù)；κ為相界面曲率，即

(15)

2.4 網(wǎng)格劃分

因結(jié)構(gòu)性網(wǎng)格生成速度快，網(wǎng)格質(zhì)量好，數(shù)據(jù)結(jié)構(gòu)簡單[13]，所以網(wǎng)格劃分采用正方形結(jié)構(gòu)性網(wǎng)格.為了獲得較高精確的計算結(jié)果、并且降低計算成本，同時在氣泡周圍采用了自適應(yīng)網(wǎng)格法(AMR).自適應(yīng)網(wǎng)格生成時，把氣相體積分?jǐn)?shù)作為變量，對氣液交界面處的網(wǎng)格進(jìn)行加密，即在氣相體積分?jǐn)?shù)大于0且小于1時的網(wǎng)格，進(jìn)行局部加密，網(wǎng)格劃分如圖2所示.

圖2 網(wǎng)格劃分

2.5 數(shù)值方法

計算采用非穩(wěn)態(tài)的層流模型，速度和壓力場耦合采用PISO算法，選用最小二乘法離散單元中心的變量梯度，壓力項選用PRESTO格式離散，動量和能量方程選用二階迎風(fēng)格式進(jìn)行離散.壓力、動量、體積力和能量的亞松弛因子分別設(shè)為 0.3、0.7、1和1，時間步長設(shè)置為5×10-5s，以保證較好的收斂性和穩(wěn)定性.計算時先進(jìn)行初始化，使計算域內(nèi)充滿靜止液體，然后再把直徑為d的氣泡放到指定位置，進(jìn)行兩相流計算.

3 網(wǎng)格無關(guān)性和結(jié)果準(zhǔn)確性驗證

3.1 網(wǎng)格無關(guān)性驗證

為了保證網(wǎng)格尺寸不對計算結(jié)果產(chǎn)生影響同時又能降低計算量，采用三種不同的網(wǎng)格對最復(fù)雜工況(Eo=80和Ga=14.5)的計算結(jié)果進(jìn)行了檢驗.自適應(yīng)網(wǎng)格的初始尺寸都是d/6，初始網(wǎng)格數(shù)均為96000.為了減少計算量，采用了自適應(yīng)網(wǎng)格劃分技術(shù).即在計算過程中，對氣液界面所臨近的網(wǎng)格進(jìn)行了局部加密.分別嘗試加密了4次(Grid1)、3次(Grid2)和2次(Grid3)，得到氣液界面位置處最小的網(wǎng)格尺寸分別是0.375/d、0.75/d、1.5/d.三種網(wǎng)格下氣泡速度隨時間的變化如圖3所示，可以看出Grid1和Grid2的計算結(jié)果基本一致.結(jié)合計算量和準(zhǔn)確性的考慮，最終選擇Grid2網(wǎng)格系統(tǒng)進(jìn)行計算.

圖3 三種網(wǎng)格系統(tǒng)的氣泡速度對比

3.2 結(jié)果準(zhǔn)確性驗證

為了描述壁面?zhèn)鳠嵝Ч暮脡?，給出了對流傳熱系數(shù)的定義：

(16)

公式中：αw為存在氣泡上升工況的局部傳熱系數(shù)；T1為熱壁面溫度；T0為流體的平均溫度.

因努塞爾數(shù)可以表征對流傳熱的強(qiáng)度，下面給出有氣泡上升時壁面局部瞬時努塞爾數(shù)Nu的定義式：

(17)

公式中：L為模型的寬度，m；λ1為液相的熱導(dǎo)率，W/(m·K).

為了驗證計算方法的可靠性，把相似工況的計算結(jié)果和文獻(xiàn)[14]的結(jié)果進(jìn)行了對比，如圖4所示.可見，橢圓形氣泡在固壁附近上升時的Nu隨計算時間的變化為一條平滑的曲線，努塞爾數(shù)隨計算時間的增加逐漸減小、最后趨于一個穩(wěn)定的值，這與文獻(xiàn)[14]得到的結(jié)果一致.

圖4 本文計算結(jié)果和文獻(xiàn)[14]對比

4 結(jié)果分析和討論

為了便于對比分析，定義了局部瞬時無量綱壁面努塞爾數(shù)

(18)

公式中：Nuz為沒有氣泡上升時壁面的瞬時努塞爾數(shù).

4.1 氣泡上升對壁面?zhèn)鳠岬挠绊?/h3>

4.1.1 不同Eo數(shù)下氣泡對壁面?zhèn)鳠岬挠绊?/p>

首先研究了給定重力水平Ga=14.5(g)和壁面距離l=1.5d下，表面張力對傳熱的影響.根據(jù)公式(1)可知，表面張力的影響可用Eo數(shù)表征.在給定重力水平Ga=14.5(g)下，不同Eo數(shù)所對應(yīng)的氣泡形狀如圖5(a)所示.可見，Eo數(shù)越大、氣泡形變越厲害，和文獻(xiàn)[15]的計算結(jié)果基本吻合.

圖5 不同Eo數(shù)下氣泡形狀及對應(yīng)的

4.1.2 不同Ga數(shù)下氣泡對壁面?zhèn)鳠岬挠绊?/p>

其次，在給定的Eo數(shù)和壁面距離l=1.5d下，研究了重力水平(即Ga數(shù))對壁面?zhèn)鳠崮芰Φ挠绊?常重力1 g對應(yīng)的Ga=14.5，0.75 g對應(yīng)的Ga=12.6，0.5 g對應(yīng)的Ga=10.3.

圖6 不同Ga數(shù)下氣泡形狀及對應(yīng)的

對于相同的Eo數(shù)，盡管重力水平(即Ga數(shù))對氣泡的形狀影響不大，但會對氣泡的上升速度產(chǎn)生重要影響.重力水平越高(即Ga數(shù)越大)，相同形狀氣泡的上升速度越大，其對氣泡上升經(jīng)過的流場擾動越大，導(dǎo)致傳熱強(qiáng)化區(qū)的范圍越寬，第一峰值也越大.第二峰值位置的偏移，是因為相同時間下，不同Ga數(shù)下氣泡上升的位置不同造成的；其小基本沒發(fā)生變化，是因為氣泡的形狀和橫向尺寸變化輕微，不足以引起氣泡熱壁面一側(cè)的流場發(fā)生明顯的變化.

4.1.3 氣泡不同初始位置對傳熱的影響

圖7 氣泡不同初始位置l對壁面的影響(Ga=14.5)

圖8 氣泡上升過程對壁面的影響

另外，從圖8也可以看出，第二峰值所對應(yīng)的氣泡方位是變化的.當(dāng)計算時間為0.1 s時，第二峰值對應(yīng)氣泡的赤道位置；在0.2 s時，對應(yīng)氣泡的底部；隨著氣泡上升距壁面距離的增大，其逐漸下移到氣泡的下游位置；當(dāng)計算時間大于0.6 s時，氣泡上升過程趨于穩(wěn)定、距壁面的距離不再增大，此時第二峰值所對應(yīng)的氣泡方位也趨于穩(wěn)定，第二峰值的大小也不再變化.這是因為氣泡對流場的擾動需要更長的時間才能對熱壁面附近的對流傳熱產(chǎn)生影響，即流動和傳熱之間存在滯后效應(yīng)[14].

4.2 傳熱強(qiáng)化機(jī)理分析

為了理解氣泡加入后對傳熱的影響機(jī)理，下面分別分析了氣泡加入后對溫度場和流場的影響.

4.2.1 氣泡上升對溫度場的影響

從圖9可以看出，當(dāng)氣泡的Eo數(shù)不同時，氣泡對熱邊界層的影響程度是不同的.為了更好地比較熱邊界層厚度的變化，圖10給出了熱邊界層最薄處沿橫向的溫度分布.從圖10可以看出，隨著Eo的增大，溫度開始變化的位置越來越靠近熱壁面，這說明熱邊界層逐漸變薄；還觀察到Eo=20和Eo=80的曲線十分相似，即它們的熱邊界層厚度十分接近.熱邊界層厚度的關(guān)系和4.1.1節(jié)中得出的結(jié)論一致：Eo數(shù)越大，氣泡引起的熱邊界層越薄，所引起的傳熱強(qiáng)化程度越顯著.

圖9 不同Eo數(shù)下氣泡周圍的溫度場

圖10 最薄熱邊界層處，溫度沿x方向的變化

4.2.2 氣泡上升對速度場的影響

圖11 不同Eo數(shù)下流場的速度矢量圖

圖12 不同Ga數(shù)下流場的速度矢量圖

圖13給出Eo=80和Ga=14.5(1 g)時，不同氣泡初始位置(l)下氣泡的速度矢量圖.氣泡初始位置的改變，和重力水平類似，影響的是氣泡初始位置附近的流動狀況.l=1.5d工況氣泡左側(cè)渦的強(qiáng)度比l=2.5d工況的強(qiáng)；同樣氣泡右側(cè)渦的強(qiáng)度變化不大.這也符合4.1.3節(jié)得出的結(jié)論，氣泡距離熱壁面的距離越小，對流體的擾動越大，氣泡初始位置附近的傳熱強(qiáng)化程度也就越顯著，但對氣泡位置附近的傳熱強(qiáng)化程度影響甚微.

圖13 氣泡不同初始位置對應(yīng)的流場速度矢量圖

綜上分析，氣泡之所以能改變壁面附近流體的對流傳熱效果是因為氣泡在近壁面上升時改變了近壁面區(qū)域的流場分布，形成旋向相反、尺度不同的對渦，改變氣泡附近和其所過區(qū)域流體的流動狀況，引起近壁熱邊界層厚度的變化，從而改變了近壁面流體的對流換熱效果.

5 結(jié) 論

本文采用VOF法研究了靜止液體中單氣泡上升對近壁面?zhèn)鳠岬挠绊懀敿?xì)研究了表面張力、重力水平和氣泡初始位置對熱壁面附近對流傳熱程度的影響.可以總結(jié)出如下幾點結(jié)論：

(1)Eo數(shù)越大，氣泡的橫向尺寸越大，其上升過程對流體的擾動越大，因此對流傳熱強(qiáng)化效果越好；

(2)氣泡形狀一定時，重力水平(即Ga)越大，氣泡上升過程對流體的擾動越大，導(dǎo)致氣泡初始位置附近區(qū)域的對流傳熱強(qiáng)化效果越好；

(3)氣泡的起始位置主要影響初始位置附近處的對流傳熱程度，氣泡距熱壁面越近，初始位置處的傳熱強(qiáng)化效果越好；

(4)氣泡上升過程產(chǎn)生旋向相反的對渦影響了冷熱流體的流動狀態(tài)，從而引起熱邊界層厚度的變化，最終改變了壁面附近流體的對流傳熱效果.