張一夫，苑宇，李維重

(1.遼寧工程技術大學 安全科學與工程學院，遼寧 葫蘆島 125105;2.大連交通大學 機械工程學院，遼寧 大連 116020;3，大連理工大學 能源與動力學院，遼寧 大連 116023)

0 引言

氣泡現(xiàn)象普遍存在于自然界和工程實際中［1-2］，廣泛用于環(huán)境、能源、化工及水利工程等領域，如液－氣燃料燃燒問題中的液－氣界面、水輪機和水泵的空化空蝕、船舶螺旋槳水流、氣液化學反應、廢水處理、破碎的波浪在水面下產(chǎn)生的大量氣泡等.氣泡在液體中的上升運動過程是一個非線性、復雜、不穩(wěn)定的動力過程.氣泡在變形過程中還常常伴隨著融合、破碎、翻轉、擺動等一系列的復雜運動，并且氣泡的運動、變形以及尾渦還會對周圍流場的速度、壓力場產(chǎn)生影響.因此，對豎直流道內氣泡上升過程的研究是一項具有挑戰(zhàn)性和現(xiàn)實意義的研究課題.

目前，常用的追蹤氣泡的方法主要有:VOF、Level Set、Front Tracking、Phase Field 等方法［3-8］.相比較而言，Level Set方法不僅能夠比較準確地追蹤運動界面，而且無須進行復雜繁瑣的界面重構技術，容易編程，具有較大的通用性，另一方面，可以采用高精度高分辨率的差分格式，對物質界面定位更加準確.但由于Level Set方法是一個非守恒的算法，將Level Set方法應用于不可壓縮兩相流問題的模擬中，由于數(shù)值方法引起的數(shù)值耗散問題，必然會引起氣泡質量的變化(質量損失或增加)，使相界面逐漸偏離正確的位置，進而使流場失真;另一方面，考慮到Level Set方法在全場進行重新初始化時，會由于邊界初網(wǎng)格非正交性引起的數(shù)值耗散.因此，本文作者引入體積修正法和局部重新初始化方法對臨近并排雙氣泡的上升運動過程進行了數(shù)值模擬，成功地解決了上述問題.

1 數(shù)學模型

1.1 基本假設

本模型基于如下假設:

(1)流體模型為低雷諾數(shù)層流模型;

(2)用均一的單流體模型計算變密度、變黏度系數(shù)的不相混兩相流;

(3)液相和氣相都視為不可壓縮流體.

1.2 控制方程組

本文采用的無量綱控制方程組包括:

(1)不可壓縮流體連續(xù)性方程

(2)考慮表面張力和浮力的動量方程

(3)采用Level Set方法追蹤兩相界面的輸運方程

式中，u為速度矢量;ρ為密度;μ為動力黏度;p為壓力;B為浮力項，B=(1－ρ)/Fr;F為表面張力項，F(xiàn)=－ kδs(φ)n/We;n 為界面單位法向矢量，n=－▽φ/|▽φ|;k為界面的曲率k=▽·n=▽·(－▽φ/|▽φ|);δs為與界面有關的Dirac分布;φ為Level Set函數(shù).無量綱數(shù)

1.3 界面追蹤方法及兩相物理量參數(shù)處理

氣-液兩相界面的區(qū)分是采用一個光滑的Level Set函數(shù)φ(x，t)來實現(xiàn)的.當光滑函數(shù)φ等于零時，該光滑函數(shù)的變量x就代表了分界面上的點，當φ值大于零為界面外部，小于零為界面內部.

為了避免在界面附近由于大密度比及大粘度系數(shù)比引起的數(shù)值不穩(wěn)定，采用修正的Heaviside函數(shù)來光滑兩相參數(shù).因此，無量綱密度和黏度系數(shù)由如下公式得出

式中，λρ= ρg/ρl和λμ= μg/μl分別為密度比和黏度系數(shù)比.

光滑的Heaviside函數(shù)H(φ)定義如下

式中，ε為數(shù)值彌散界面的可調帶寬.

由于數(shù)值方法的內在因素，隨著計算的進行，Level Set函數(shù)往往無法保持符號距離函數(shù)(Signed Distance Function，SDF)的特性.因此，為了使φ(x，t)保持這種特性，需要采用一種方法來將Level Set函數(shù)重置為符號距離函數(shù).利用符號距離函數(shù)的如下性質:

若φ為符號距離函數(shù)，則有:

重新初始化的偏微分方程為

式中，需將sign(φ0)光滑化為為重新初始化之前，上一時層的最后計算值.

考慮到Level Set方法在全場進行重新初始化時，會由于邊界處網(wǎng)格非正交性引起的數(shù)值耗散，因此，將Level Set重新初始化方法局限在了界面附近一個很小的范圍內，可以大大提高Level Set計算的精度.因此，本文采用類似于Peng等［9］提出的窄帶計算的局部重新初始化方法，將兩相界面定義在一個有限的區(qū)間(－β≤φ≤β)內.界面仍采用符號距離函數(shù)的定義方法(如圖1所示)，對φ值進行重新初始化計算.

圖1 界面函數(shù)初始值定義(當β=0.15時)

為保證氣泡在計算過程中具有很好的質量守恒性，采用體積校正方法［10］對氣泡體積進行修正.這種方法的主要思想是通過不連續(xù)相(氣泡或液滴)體積的盈虧來調整界面位置:

假設氣泡是球形的，半徑的增量為

其中，C1為模擬體系的維數(shù).由于分散相體積的增加(或減少)對應于從界面到分散相中心距離的增加(或減少)，所以可取Level Set函數(shù)的校正量正比于Δr.因此，Level Set函數(shù)的校正量可表達為

其中，C2是一個經(jīng)驗參數(shù)，一般在0.01～0.1之間的范圍內取值，C2較大可能導致重新初始化方程不熟練，較小會降低計算效率.本文中C2取為0.01 .

這樣，Level Set方法的體積校正方程為

1.4 表面張力項處理

表面張力項采用連續(xù)表面張力模型(Contin-uous Surface Force，CSF)［11］模型計算來處理.對于Level Set方法，CSF模型的表面張力表達式為

其中，δ函數(shù)表示如下:

2 控制方程的離散和求解

分別采用三階精度的QUICK格式以及二階精度的中心差分格式，對N-S方程中對流項與擴散項進行離散.主場的壓力與速度耦合采用經(jīng)典的SIMPLE［12-13］算法求解.物理量及參數(shù)定義采用同位網(wǎng)格［14］，所有物理量都定義在網(wǎng)格節(jié)點的中心處.界面上的逆變速度采用非交錯網(wǎng)格上的Rhie-Chow動量插值方法［15］獲得，在計算界面壓力梯度時，可以有效抑制不合理的壓力場產(chǎn)生.

Level Set輸運方程采用二階上風和帶有Superbee限制算子的平均積分型TVD格式進行離散求解［16］，以保證非線性對流項數(shù)值穩(wěn)定性.重新初始化方程采用修正的Godunov方法對五階WENO差分格式進行求解，得到重新初始化值.時間導數(shù)項采用三階的TVD-Runge-Kutta方法，以保證每一時間步Level Set方程都滿足符號距離函數(shù)的特性.

采用C++語言編制代碼求解控制方程.計算流程如下:

(1)初始化函數(shù)φ為到界面的符號距離函數(shù);

(2)采用SIMPLE算法求解動力學方程，從tn時間步的速度和壓力場獲得tn+1的值;

(3)求解Level Set方程和局部重新初始化方程以保證φ為到界面的符號距離;

(4)采用體積校正法修正氣泡體積;

(5)更新tn+1時層上的流體特性參數(shù);

(6)時間遞進一步，重復(2～5)步直至最后時間步.

3 模擬結果與討論

本文用未采用體積校正法的全場重新初始化方法(方法1)和采用體積校正法的局部重新初始化方法(方法2)分別對豎直流道內并排雙氣泡的上升過程進行了數(shù)值模擬.左右邊界和底部邊界均采用無滑移邊界條件，上端為開口邊界，采用外推邊界條件.初始時，氣泡形狀為圓形，流場靜止，兩氣泡圓心間距為Δd=2.4R，氣泡初始圓心位置距底面為 3R0.模擬工況:Re=10，We=20，F(xiàn)r=1.0，ρb/ρl=1/1 000 和 μb/μl=1/100.計算網(wǎng)格為120×240.利用兩種方法分別模擬出的氣泡形態(tài)對比如圖2所示.

圖2 兩種方法模擬出的氣泡形態(tài)演變對比

從圖2中可以看出，利用上述兩種方法模擬出的氣泡形態(tài)演變過程有明顯差異.方法1模擬出的結果來看，隨著氣泡的上升，氣泡明顯變小，氣泡體積守恒性很差.而方法2模擬出的結果，氣泡守恒性很好.另外，兩種方法模擬出的氣泡上升速度也有所不同(如圖3所示)，隨著時間的推移，方法1模擬出的結果中，氣泡在經(jīng)過開始的加速階段后，速度達到一個峰值，但由于氣泡體積的逐漸減小，氣泡的速度開始下降.在整個過程中，氣泡始終都達不到一個穩(wěn)定的狀態(tài).而由方法2模擬的結果來看，當氣泡上升到一定高度后，氣泡形態(tài)不在變化，因此氣泡的速度最終也維持恒定.圖4分別對利用兩種方法模擬并排雙氣泡時，氣泡體積偏差隨時間的變化規(guī)律進行了對比.從圖中可以看出，利用方法1得出的結果中，在初始階段氣泡的體積守恒性保持的很好，但隨著時間的推移，體積偏差逐漸偏離1.0.而利用方法2獲得的結果來看，氣泡的體積偏差始終在1.0附近一個很小的范圍內波動.由此可見，采用方法2對氣泡的模擬更加準確，結果明顯好于方法1.

圖3 兩種方法模擬出的氣泡上升速度隨時間變化

圖4 氣泡體積盈虧隨時間的變化曲線

4 結論

本文通過建立氣液兩相流模型，提出改進的貼體坐標下Level Set方法，并結合考慮表面張力項采用CSF模型，對豎直流道內并排臨近雙氣泡的上升運動過程進行了數(shù)值模擬.速度與壓力的解耦采用基于同位網(wǎng)格下經(jīng)典的SIMPLE算法，準確地追蹤了氣泡上升過程，獲得了氣泡運動過程形狀，位置及速度等特性參數(shù)變化規(guī)律.采用基于體積校正法和局部重新初始化的Level Set方法，解決了Level Set方法中由于數(shù)值耗散引起的氣泡體積的非守恒問題，獲得了滿意的結果.

［1］MONTES F J，GALAN M A，CERRO R L.Mass transfer from oscillating bubbles in bioreactors［J］.Chemical Engineering Science，1999，54:3127-3136.

［2］DARMANA D，HENKET R L B，DEEN N G，et al.Detailed modelling of hydrodynamics，mass transfer and chemical reactions in a bubble column using a discrete bubble model:Chemisorption of CO2into NaOH solution，numerical and experimental study［J］.Chemical Engineering Science，2007，62(9):2556-2575.

［3］HIRT C W，NICHOLS B D.Volume of fluid(VOF)method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39(1):201-225.

［4］YOUNGS D L.Time-dependent multi-material flow with large fluid distortion numerical method for fluid dynamics［M］.New York:Academic Press，1982:273-285.

［5］OSHER S，SETIAN J A.Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations［J］.Journal of Computational Physics，1988，79:12-49.

［6］CHEN S，MERRIMAN B，OSHER S，et al.A simples Level Set method for solving stefan problems［J］.Journal of Computational Physics，1997，135:8-29.

［7］CHARRIER P，TESSIERAS B.On front-tracking methods applied to hyperbolic systems of nonlinear conservation laws［J］.SIAM Journal on Numerical Analysis，1986，23(3):461-472.

［8］WANG S L，SEKERKA R F.Algorithms for phase field computation of the dendritic operating state at large supercoolings［J］.Journal of Computational Physics，1996，127(1):110-117.

［9］PENG D，MERRIMAN B，OSHER S，et al.A PDE-based fast local Level Set method［J］.Journal of Computational Physics，1999，155(2):410-438.

［10］李向陽，王躍發(fā)，禹耕之，等.改進水平集法模擬不可壓縮兩相流動時質量守恒的體積校正法［J］.中國科學 B 輯:(化學)，2008，38(7):636-644.

［11］BRACKBILL J U，KOTHE D B，ZEMACH C.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100:335-354.

［12］PATANKAR S V.Numerical Heat Transfer and Fluid Flow［M］.New York:Hemispher，1980.

［13］陶文銓.數(shù)值傳熱學［M］.西安:西安交通大學出版社，2001.

［14］PERICCM，KESSLER R，SCHEUERER G.Comparison of finite-volume numerical methods with staggered and colocated grids［J］.Computers and Fluids，1988，16(4):389-403.

［15］RHIE C M，CHOW W L.A numerical study of t he turbulent flow past an isolated airfoil wit h t railing edge separation［J］.AIAA Journal，1983，21(11):1525-1532.

［16］劉儒勛，王志峰.數(shù)值模擬方法和運動界面追蹤［M］.合肥:中國科學技術大學出版社，2001.