單氣泡上升對近壁面傳熱影響的數值研究

張文峰，龐明軍

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

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

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

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

1 物理模型和工況設置

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

圖1 計算區域

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

氣泡在靜止液體內上升時，主要受重力、黏性力和表面張力的影響，這三種力可由無量綱參數E?tv?s(Eo)數和Gallilei(Ga)數表示.Eo數是重力和表面張力的比值，Ga數是重力和粘性力的比值，其定義式見公式(1)和公式(2).通過改變表面張力、重力水平和氣泡距固壁的距離，分別控制Eo數、Ga數和l，來改變氣泡運動狀態，從而研究氣泡對壁面對流傳熱的影響.因此，設計了三組工況：(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.本文計算工況的物性參數，如表1所示.

表1 氣液物性參數

2 控制方程和數值方法

2.1 控制方程

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

連續性方程：

?·u=0，

(3)

動量方程：

(4)

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

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

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

(5)

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

(6)

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

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

(7)

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

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

(8)

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

.

(9)

2.2 VOF方法

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

(10)

相函數f可定義為

(11)

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

(12)

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

(13)

2.3 連續表面張力模型

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

Fσ=σκδn，

(14)

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

(15)

2.4 網格劃分

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

圖2 網格劃分

2.5 數值方法

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

3 網格無關性和結果準確性驗證

3.1 網格無關性驗證

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

圖3 三種網格系統的氣泡速度對比

3.2 結果準確性驗證

為了描述壁面傳熱效果的好壞，給出了對流傳熱系數的定義：

(16)

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

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

(17)

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

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

圖4 本文計算結果和文獻[14]對比

4 結果分析和討論

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

(18)

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

4.1 氣泡上升對壁面傳熱的影響

4.1.1 不同Eo數下氣泡對壁面傳熱的影響

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

圖5 不同Eo數下氣泡形狀及對應的

4.1.2 不同Ga數下氣泡對壁面傳熱的影響

其次，在給定的Eo數和壁面距離l=1.5d下，研究了重力水平(即Ga數)對壁面傳熱能力的影響.常重力1 g對應的Ga=14.5，0.75 g對應的Ga=12.6，0.5 g對應的Ga=10.3.

圖6 不同Ga數下氣泡形狀及對應的

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

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

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

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

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

4.2 傳熱強化機理分析

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

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

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

圖9 不同Eo數下氣泡周圍的溫度場

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

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

圖11 不同Eo數下流場的速度矢量圖

圖12 不同Ga數下流場的速度矢量圖

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

圖13 氣泡不同初始位置對應的流場速度矢量圖

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

5 結 論

本文采用VOF法研究了靜止液體中單氣泡上升對近壁面傳熱的影響，詳細研究了表面張力、重力水平和氣泡初始位置對熱壁面附近對流傳熱程度的影響.可以總結出如下幾點結論：

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

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

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

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