基于VOSET方法模擬并排氣泡的上升過程

李隆鍵，張 磊，朱文冰，申憲文

（重慶大學動力工程學院，重慶 400030）

文章編號：1001?246X（2015）05?0545?08

基于VOSET方法模擬并排氣泡的上升過程

李隆鍵，張 磊，朱文冰，申憲文

（重慶大學動力工程學院，重慶 400030）

采用VOSET方法和耦合表面張力模型的N?S方程，模擬豎直通道內并排氣泡對的上升過程，模擬與實驗結果吻合較好.重點研究表面張力系數對并排氣泡上升軌跡和速度的影響.結果表明，隨著表面張力系數的變化，并排氣泡對的上升過程出現三種類型：兩氣泡融合，兩氣泡反復靠近、遠離但未融合，兩氣泡碰撞反彈后逐漸遠離.在未融合的情況下，并排氣泡對的上升軌跡關于通道中心線對稱，左右兩個氣泡的上升速度基本一致，水平速度大小相同，方向相反.

并排氣泡；數值模擬；VOSET

0 引言

多相流廣泛存在于能源、動力和環境等工程設備和工藝中，如核電站的蒸汽發生器、火箭發動機、鼓泡反應器等工業設備.這些過程往往涉及到多個氣泡的融合、破碎等復雜變化.氣泡之間的相互作用顯著影響氣泡的形狀和運動特性，因此，準確研究氣泡之間的相互作用，如各個氣泡速度的波動，氣泡的融合和反彈等，對研究多相流尤其是泡狀流的運動特性具有重要意義.

Sanada等［1］實驗研究了靜止流體中一對輕微變形的并排氣泡的反彈和融合特性，使用高速照相機對氣泡的運動軌跡和尾跡進行了可視化研究.Duineveld等［2］實驗研究了高Re數下氣泡對在純水中和含有活性劑的水中的相互作用，分析了Weber數對氣泡對運動的影響.

隨著數值模擬技術的發展，越來越多的學者采用數值方法來研究氣泡間的運動行為.張淑君等［3］采用結合PLIC界面重構的VOF方法，模擬了不同放置方式下氣泡之間的相互作用，重點研究了周圍流體黏性及氣泡間距對氣泡融合的影響.趙知辛等［4］采用Level Set方法對水下兩個直徑為8mm的氣泡的運動過程進行數值模擬.Chen等人［5］使用MPS方法（移動粒子半隱式法）對靜止液體中氣泡對的上升進行了數值模擬，研究了水平并排氣泡對和豎直同軸兩個氣泡的融合過程.王太等［6］采用三維的VOF方法模擬研究了豎直同軸兩個氣泡的融合過程，分析了液體粘度和表面張力對同軸氣泡融合的影響.

直徑為2mm左右的氣泡廣泛存在于泡狀流動等實際過程中［2］，Wu等［7］實驗研究了直徑為1mm～2mm的氣泡在水中上升的形狀和路徑，實驗發現直徑小于1.5mm的氣泡呈直線上升，更大的球形氣泡則呈“之”字形或螺旋狀上升.Sone等［8］實驗研究了靜止流體中一對“之”字形氣泡碰撞過程中氣泡和周圍流體的運動.在以往的模擬中，主要研究的是較大氣泡間的相互作用，而小氣泡間的相互作用研究較少.本文采用數值模擬方法，研究小氣泡的相互作用過程，該范圍內氣泡的Re數較大，且呈“之”字形上升.

氣泡的融合過程涉及到復雜的相界面變化，對數值方法的要求較高.VOF方法缺點是難以計算曲率及處理界面附近突變的物理量，Level Set方法不能保證質量守恒，這兩種方法都有一定的局限性.結合二者方法的優點，Sun［9］等人提出VOSET方法，該方法被證明能夠模擬復雜的界面拓撲變化，且能保證質量守恒.Sun分別通過三個二維的不相溶，不可壓的兩相問題對VOF，Level Set和VOSET方法進行計算分析，VOSET的界面曲率、精度、整場的密度和粘度的光順性，質量守恒特性均優于其它兩種方法.葉政欽［10］采用VOSET方法針對兩個同軸氣泡及多個氣泡的上升和融合過程進行了模擬，分析了復雜兩相流動問題中VOSET方法的質量守恒特性、幾何計算方法的可靠性及求解過程的穩定性.

采用VOSET方法，結合考慮表面張力的氣液運動基本方程，數值模擬水平并排放置的一對半徑R＝0.9mm氣泡對在豎直通道內的上升過程，研究氣泡對相互作用規律.

1 數學模型

1.1 控制方程

對粘性不可壓縮兩相流動，其控制方程為

式中，g為重力加速度，u為速度矢量，ρ為流體密度，μ為流體的粘性.

在界面附近，考慮表面張力的作用，對動量方程進行修正，采用CSF模型［11］（連續表面張力模型），將界面上的表面力轉化為相應的控制容積中的體積力.利用Level Set函數［12］?對動量方程進行修正，修正后的動量方程

式中，σ為表面張力系數，κ（?）為曲率，

Hε（?）為光滑的Heaviside函數，定義為

其中ε表示界面處光滑區間單側的寬度，一般ε＝1.5Δ，Δ表示最小的網格的寬度.

為了使物性在界面上連續光滑變化，流體物性可借助于Heaviside函數做插值計算.這樣整場的密度和動力粘度可表示為

在界面追蹤方法中，VOF方法和Level Set方法是目前運動界面問題中應用較為廣泛的數值方法.VOF方法中引入了流體體積函數C的概念，對于兩相流動流體分為目標流體和非目標流體，流體體積函數C等于一個單元內目標流體與該單元體積之比.對于不可壓縮流動，流體體積函數C的輸運方程為

Level Set方法［12］的主要思想是把隨時間運動的物質界面看作某個函數?（x，t）的零等值面.在每個時刻t，我們只要求出函數的值，求出其零等值面的位置，也就能得到運動界面的位置.首先需構造函數?（x，t），使得在任意時刻t運動界面Γ（t）恰是?（x，t）的零等值面．一般可取?（x，t）為x點到界面Γ（0）的符號距離.

Γ（t）表示t時刻的自由界面，其中d（x，Γ（t））表示區域中任一點到Γ（t）的距離．在任意時刻t，對于活動界面Γ（t）上的任意點x，?（x，t）＝0，即符號距離函數的零等值面即為自由界面.

VOSET方法的主要思想是：流體體積函數C的輸運仍然采用VOF中C的輸運方程，并結合PLIC界面重構方法來更新下一時刻的流體體積函數C.PLIC界面重構方法的思路是：首先在單個網格內用直線段來重構界面，逼近真實界面，然后計算流過該網格各邊界的目標流體的體積流量，從而計算出下一時刻的流體體積函數.

與一般的VOF方法中直接采用C來加權計算界面附近的參數不同的是，VOSET方法在每個時刻通過已知的流體體積函數C，采用幾何迭代方法計算單元（i，j）距其計算范圍內網格上所有界面的最小距離后，通過比較可以得出單元（i，j）距離界面的最短距離d，然后通過式（9）計算符號距離函數?，利用此?函數求得界面附近的曲率和物性參數，具體實施方法參見文獻［9］.

1.2 計算區域和邊界條件的選擇

為了滿足計算區域的要求和節省計算資源，根據文獻［5］的結論，當氣泡中心和壁面的距離大于5倍半徑時，壁面的影響可以忽略不計，而且高度要足夠大使氣泡能夠達到終端速度和終端形狀.綜合考慮，本文選擇12mm×24mm的矩形區域，各邊界條件采用無滑移邊界條件.

1.3 計算流程

本文采用有限體積法對流場進行離散求解，利用VOSET方法捕捉界面.求解步驟為：①對流場進行初始化，給出初始的流體體積函數C的分布；②利用C的分布求得符號距離函數?的分布，按照式（4）～（7）計算相關的參數，使用IDEAL［13］數值算法求解控制式（1）～（3），獲得收斂的速度場；③采用VOSET方法求得下一時刻流體體積函數C的分布；④返回步驟①，進入下一時間步的計算.

2 結果與分析

2.1 模型的驗證

首先對靜止液體中單個氣泡的上升過程進行數值模擬，以驗證模型的正確性.計算區域為12mm×24mm的矩形區域，區域內充滿靜止液體.液體的物性參數：ρl＝1 000 kg·m－3，μl＝0.001 Pa·s；氣體的物性參數：ρg＝1.1 kg·m－3，μg＝1.8×10－5Pa·s，表面張力系數σ0＝0.072 8 N·m－1，E o＝0.5，M o＝2.5× 10－11.直徑d＝0.001 88 m的氣泡從（0.006 m，0.002 m）的位置釋放.經網格無關性驗證，網格數取90× 180，計算得到的單個氣泡上升過程如圖1所示，氣泡由初始的球形連續變化為扁橢球狀，與Fan［14］給出的氣泡形狀圖譜吻合.該氣泡的最終上升的速度為22.4 cm·s－1，與文獻［7］中20 cm·s－1的速度很接近，證明了該模型的正確性.

圖1 單個氣泡上升模擬Fig.1 Simulation of single bubble rising in quiescent liquid

2.2 水平氣泡對的上升

一個12mm×24mm大小的豎直通道，距底部3mm處，兩個直徑為1.8mm的球形氣泡關于通道中心線對稱并排放置，氣液物性和上述的一致.兩氣泡中心的初始水平距離為2.2mm，內邊緣相距0.4mm.隨著氣泡的上升，兩個氣泡相互靠近，之后接觸、融合，變成一個更大的氣泡繼續上升.融合后的氣泡沒有穩定的形狀，在上升的過程中伴隨著形狀的震蕩，模擬結果如圖2所示，數值模擬的結果與Duineveld［2］的實驗結果吻合較好（圖3）.其實驗條件：10cm×10cm×50 cm立方玻璃水箱，水溫20℃，玻璃箱下面有兩個相距一定距離的注射器，磁力閥推動注射器里面的柱塞運動，從而產生氣泡.實驗中采用一系列的措施基本保證兩個氣泡能夠同時釋放，而且初始為球狀且沒有形狀的波動，利用高速攝像機拍出氣泡運動圖.類似的氣泡融合過程也能在Sanada的實驗［1］中觀察到.

圖2 并排氣泡融合數值模擬Fig.2 Simulation of side by side coalescence between two identical bubbles

圖3 Duineveld（1998）實驗結果Fig.3 Experimental observations of side by side coalescence between two identical bubbles

在兩氣泡靠近融合的過程中，兩氣泡融合前對稱上升，上升速度基本相同.氣泡上升的速度隨時間變化如圖4（a）所示.兩氣泡在融合前的上升過程中，中間會形成一個低壓區，從而兩個氣泡會相互靠近，如圖4 （b）所示.可以看出在接觸之前的上升過程中，兩氣泡對稱上升，速度逐漸增大.接觸時，當兩個氣泡開始融合到接觸面積達到最大時，其體積大大增加，在上升方向產生的阻力達到最大，因此氣泡上升速度呈現下降趨勢.而當融合后的形狀開始收縮時，氣泡上升速度亦開始反彈上升，最后呈波動增大.氣泡在融合過程中上升速度出現了上升－下降－上升的波動現象，這個結果與Sanada［1］實驗的結果相吻合.融合后的氣泡形狀波動變化，導致受力不均勻，從而導致速度的波動.這與張淑君等［3］對低Re數下的氣泡對融合上升最終會有一個穩定的形狀和上升速度不同.這是因為高Re數下，氣泡尾跡的作用使得氣泡沒有穩定的形狀.趙知辛［4］模擬水中相同相對中心距下8mm的氣泡，兩氣泡沒有發生融合，可見小氣泡更易融合.

2.3 表面張力對氣泡融合的影響

氣液界面的表面張力系數是重要的物性參數，它與諸多工業問題密切相關，而表面張力系數又不是一個固定的數值，隨工質的溫度，壓力等因素變化.本文采用數值模擬方法來研究表面張力系數的變化對氣泡作用的影響.

圖4 兩氣泡上升速度和速度矢量Fig.4 Vertical velocity and velocity vector distribution of two bubbles

2.3.1 改變表面張力，氣泡融合的情況

改變表面張力系數，分別取σ1＝0.1σ0，σ0，5σ0，10σ0，這四種情況下，兩氣泡在上升的過程中發生融合.氣泡上升速度隨時間變化如圖5所示.σ1＝0.1σ0時兩個氣泡融合后速度波動不大，如圖5（a）所示.σ1＝5σ0時，兩個氣泡融合后速度變大，波動也更為劇烈.σ1＝10σ0時的情況和σ1＝5σ0時的類似.兩個氣泡在接觸后融合成一個更大的氣泡，在上升過程中，氣泡的形狀震蕩，速度也在波動，從圖5中可以看出，隨著表面張力系數σ的增大，融合后的氣泡速度震蕩的更加劇烈，幅度也更大，且平均上升速度呈上升趨勢.

圖5 氣泡上升速度隨時間變化Fig.5 Vertical velocity of coalescence bubbles

2.3.2 改變表面張力，氣泡分離的情況

改變表面張力系數，分別取σ2＝0.04σ0，0.02σ0，0.01σ0，0.005σ0，這四種情況下，兩個氣泡獨立上升，均不會融合.在以往的研究中，大都是針對氣泡的形狀變化進行模擬研究，很少對氣泡運動軌跡進行研究，本文采用數值方法，再現了兩氣泡在相互作用過程中質心軌跡的變化.未融合情況下，兩個氣泡質心的運動軌跡如圖6所示.為了研究氣泡之間的相互作用對單個氣泡運動的影響，圖6也給出了右氣泡單獨存在時的運動軌跡.從圖中可以看出，隨著表面張力系數的減少，單個右氣泡的上升軌跡水平波動越來越小，到后面近似呈直線上升，而且氣泡之間的相互作用加劇了右氣泡的波動程度.從圖中還可以看出，單氣泡上升的距離更大，說明另一個氣泡的存在減小了氣泡的平均上升速度.

圖6 氣泡質心軌跡Fig.6 Trajectory of single bubble and a pair of bubbles

并排氣泡對的上升過程出現了兩種運動類型：反復靠近?遠離?靠近但未融合（σ2＝0.04σ0），如圖6（a）所示；碰撞反彈后相距越來越遠（σ2＝0.02σ0，0.01σ0，0.005σ0），這三種情況下運動軌跡類似，如圖6（b）所示.表面張力較小時，兩個氣泡在上升過程中出現靠近－遠離－再靠近的周期性現象，這與李彥鵬模擬水中8mm氣泡對的研究結果［15］相類似.可見，氣泡直徑增大與表面張力系數減小的效果類似，這是因為氣泡直徑的增大同樣會引起所受表面張力的減小.表面張力繼續變小時，兩個氣泡在上升過程中水平距離越來越大，并沒有出現靠近的現象，如圖6（b）所示，在Sone［8］的實驗中也觀察到了同樣的現象.氣泡對之間反復出現靠近－遠離－靠近的現象，主要是因為高Re數下，氣泡尾部出現了渦旋的脫落，它們之間的相互作用導致了氣泡對的運動軌跡出現波動［1］，如圖7所示.而隨著表面張力系數的減小，氣泡尾跡的相互作用表現為排斥力作用，導致兩氣泡逐漸遠離.

圖7 兩氣泡速度矢量（σ2＝0.04σ0）Fig.7 Velocity vector distribution of two bubbles

并排氣泡和單個右氣泡獨立存在時的上升速度如圖8所示.可以看出，對于單個右氣泡的獨立上升，氣泡的上升速度先增大，后呈波動狀態.對于并排上升氣泡對而言，左右兩個氣泡上升速度基本完全一致均呈波動狀態，但比單個氣泡獨立上升時的速度小.可見氣泡之間的相互作用減小了氣泡的上升速度.

圖8 氣泡上升速度隨時間變化Fig.8 Vertical velocity of single bubble and a pair of bubble

從圖8中可以看出，隨著氣泡的上升，并排氣泡上升速度不斷增大，后面開始震蕩，而且隨著表面張力系數的減少，震蕩的頻率和幅度均下降.上升初期，兩個氣泡不斷靠近，上升速度也在逐漸增大，相互靠近到某一距離時，兩氣泡開始反彈，此時，上升速度顯著下降，之后呈波動狀態.

氣泡水平速度隨時間變化如圖9所示.右氣泡單獨上升時，其水平速度并不為零而且呈波動狀態，從而解釋了其上升軌跡水平方向有波動的原因.從圖中可以看出，并排上升兩氣泡的水平速度大小相等，方向相反，而且同樣也呈現出波動狀態.隨著表面張力系數的減小，單個氣泡水平速度波動的幅度和頻率均減小，兩個并排氣泡的水平速度也是如此.另一個氣泡的存在加劇了右氣泡的水平波動，這也解釋了兩氣泡并排上升時對稱波動上升的軌跡特征.

圖9 氣泡水平速度隨時間變化Fig.9 Horizontal velocity of single bubble and a pair of bubbles

2.3.3 結果分析

從前面兩種情況可以看出，存在一個臨界的表面張力系數σcr，因此選擇更多的表面張力值，分別取σ＝0.09σ0、008σ0、0.007σ0，當 σ＝0.09σ0時氣泡融合，其它兩種情況下兩個氣泡不會融合，可見 σcr約為0.09σ0．在常見的液態工質中，液態金屬和熔融鹽以及水的表面張力比較大，是大于這個范圍的，這些工質中的氣泡更容易融合.汽油、乙醚等有機溶液表面張力較小，在這些工質中氣泡不易融合.

在兩氣泡的相互作用過程中，液橋是一個重要的因素.液橋力主要是由氣－液相界面的表面張力和液橋內外的壓力差引起的［16］.因此，改變表面張力系數，兩氣泡的作用過程會出現不同的結果.

3 結論

1）采用VOSET方法模擬了并排氣泡對在豎直通道內的相互作用過程，模擬結果與相關文獻和實驗吻合較好.

2）改變表面張力系數，氣泡對出現了三種運動類型.較大的表面張力系數下，氣泡發生融合，融合后形成一個更大的氣泡上升，融合后的氣泡沒有穩定的形狀，也沒有穩定的上升速度；中等表面張力系數時，兩氣泡沒有融合，在上升過程中出現反復靠近－分離－再靠近的運動特征；較小表面張力系數下，兩個并排氣泡未發生融合，碰撞反彈后相互遠離，排斥上升，在上升的過程中，二者的水平距離越來越大.

3）氣泡對的相互作用過程中，如果沒有融合，則兩個氣泡的上升速度基本一致，和單個氣泡的上升速度相比較，均有所下降但變化趨勢基本趨勢一致，和水平速度相比差異較大.從而可以看出，兩氣泡的相互作用主要影響氣泡對的水平運動，對豎直方向的運動影響不大.

［1］ Sanada T，Sato A，et al.Motion and coalescence of a pair of bubbles rising side by side［J］.Chemical Engineering Science，2009，64（11）：2659－2671.

［2］ Duineveld PC.Bouncing and coalescence of bubble pairs rising at high Reynolds number in pure water or aqueous surfactantsolutions［J］.Appl Sci Res，1998，58：409－439.

［3］ 張淑君，吳錘結.氣泡之間相互作用的數值模擬［J］.水動力學研究與進展：A輯，2009，（6）：681－686.

［4］ 趙知辛，王煥然，李彥鵬，等.用Level Set方法對水下兩個氣泡運動的數值模擬［J］.西安交通大學學報，2009，43 （7）：11－15.

［5］ Chen R H，Tian W X，et al.Numerical investigation on coalescence of bubble pairs rising in a stagnant liquid［J］.Chemical Engineering Science，2011，66（21）：5055－5063.

［6］ 王太，李會雄，李陽.同軸兩個氣泡融合特性的數值研究［J］.西安交通大學學報，2013，47（7）：1－6.

［7］ Wu M，Gharib M.Experimental studies on the shape and path of small air bubbles rising in clean water［J］.Physics of Fluids，2002，14（7）：L49－L52.

［8］ Sone D，Sakakibara K，et al.Bubblemotion and its surrounding liquid motion through the collision of a pair of bubbles［J］. Journal of Power and Energy Systems，2008，2（1）：306－317.

［9］ Sun D L，TaoW Q.A coupled volume?of?fluid and level set（VOSET）method for computing incompressible two?phase flows ［J］.International Journal of Heat and Mass Transfer，2010，53（4）：645－655.

［10］ 葉政欽，劉啟鵬，李星紅，等.復雜兩相流中界面追蹤方法——VOSET的性能分析［J］.化工學報，2011，62（6）：1524－1530.

［11］ Brackbill JU，Kothe D B，et al.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100（2）：335－354.

［12］ Osher S，Sethain JA.Fronts propagating with curvature dependent speed：Algorithms based on Hamilton Jacobi formulations ［J］.JComput Phys，1988，79：12－49.

［13］ Sun D L，Qu ZG，et al.An efficient segregated algorithm for incompressible fluid flow and heat transfer problems—IDEAL （inner doubly iterative efficient algorithm for linked equations）Part I：Mathematical formulation and solution procedure［J］. Numerical Heat Transfer，Part B：Fundamentals，2008，53（1）：1－17.

［14］ Fan L S，Tsuchiya K.Bubble wake dynamics in liquids and liquid?solid suspensions［M］.Boston，Massachusetts，USA：Butterworth?Heinemann，1990：263－270.

［15］ 李彥鵬，張乾隆，白博峰.豎直通道內相鄰氣泡對上升的直接數值模擬［J］.熱能動力工程，2007，22（4）：375－379. ［16］ 韓笑.氣固流化床中持液顆粒的流化特性及反應器模型研究［D］.杭州：浙江大學，2013.

Simulation of a Pair of Bubbles Rising Side by Side Using VOSET M ethod

LILongjian，ZHANG Lei，ZHUWenbing，SHEN Xianwen
（College of Power Engineering，Chongqing University，Chongqing 400030，China）

A numerical simulation was performed to investigate two bubbles rising side by side in vertical channel using VOSET （Coupled Volume?of?Fluid and Level Set）method and N?Sequations coupled with continuous surface force（CSF）model.Behaviors of coalescence between two identical bubbles predicted were in good agreementwith experimental results reported in literature.Effects of surface tension on trajectories and velocities during rising process are investigated.Three types ofmotion were observed depending on surface tension：Coalescence，two bubbles repeatedly attracted and bounced against each other，bounced and separated.Without coalescence，trajectory shows symmetry about channel center line.Vertical velocities of both bubbles are almost the same while magnitude of horizontal velocity with opposite direction equals.

bubble pairs；numerical simulation；VOSET

O359＋.1

A

2014－09－02；

2014－12－08

國家自然科學基金（51076172）資助項目

李隆鍵（1966－）男，博士，教授，主要從事氣液兩相流的數值模擬研究，E?mail：longjian＠cqu.edu.cn

Received date： 2014－09－02；Revised date： 2014－12－08