氣泡群動力特性模擬分析

張淑君,吳錘結(jié)

(1.河海大學(xué)淺水湖泊綜合治理與資源開發(fā)教育部重點實驗室,江蘇南京 210098;2.大連理工大學(xué)機(jī)械工程系,遼寧大連 116014)

氣泡群運(yùn)動的研究在環(huán)境工程、水利工程和航運(yùn)工程等方面具有重要的理論和現(xiàn)實意義[1-3].如污水處理廠的曝氣工藝,河道水環(huán)境治理中的人工受氧過程以及水源地水質(zhì)改善等都與氣泡群的運(yùn)動密切相關(guān).隨著氣泡數(shù)量的增加,氣泡之間的影響和作用變得更加復(fù)雜,特別是當(dāng)粒子濃度很大時對粒子的空間分布也有不同程度的影響.早期工業(yè)氣泡流特性研究只集中于體特性,如平均壓降、平均流量或平均剪切應(yīng)力等,且僅限于一維[4],而更細(xì)致的氣泡研究應(yīng)著重于氣泡本身和其引起的周圍流場的三維運(yùn)動特性,這要借助于各種界面模擬手段.Semereka等[5-6]采用邊界積分方法模擬了大量氣泡在勢流中的運(yùn)動,結(jié)果表明,氣泡上升速度在經(jīng)過一段時間后基本相同且在上升過程中氣泡水平聚集,這一結(jié)論與試驗結(jié)果不盡相同的原因是模擬過程中沒有考慮氣泡尾部漩渦的作用.Spelet等[7-8]模擬分析了湍流對氣泡群平均上升速度的影響,發(fā)現(xiàn)氣泡群在湍流中的上升速度只有其在靜水中上升速度值的一半左右.由于文中把氣泡假設(shè)成點顆粒,且沒有考慮氣泡間的相互作用,因此該結(jié)論只適用于體積分?jǐn)?shù)很小的氣泡運(yùn)動.在考慮了更多因素如慣性、黏性等,Hu等[9-10]采用非結(jié)構(gòu)網(wǎng)格模擬了100個固體顆粒的運(yùn)動,但由于在每一時間步均須重構(gòu)流場網(wǎng)格,與目前廣泛采用的單流體模型相比有很大的不足.目前以Tryggvason[11-13]為代表的鋒面跟蹤法在模擬氣泡群方面取得了一定進(jìn)展,但多數(shù)假設(shè)氣泡在上升過程中不發(fā)生碰撞、融合等現(xiàn)象.本文采用VOF(volume-of-fluid)中的PLIC(piecewise linear interface calculation)界面重構(gòu)方法研究多個三維氣泡的動力學(xué)特性及其對周圍流場的擾動.

1 控制方程及其離散

1.1 控制方程

考慮表面張力的動量方程為

式中:σ——表面張力 ;p——壓強(qiáng);μ——動力黏度;κ——界面的曲率;δs——與界面有關(guān)的 Dirac分布;n——界面單位法向矢量;D——應(yīng)力張量,滿足

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

采用VOF方法模擬氣泡界面的體積函數(shù)C滿足

對于兩相流,方程(1)中的 μ和ρ均由體積函數(shù)C決定

式中 ρ1,ρ2,μ1,μ2分別為2種流體的密度和動力黏度.

1.2 界面跟蹤

本文采用Youngs的PLIC方法在單個網(wǎng)格內(nèi)用直線段近似界面方法進(jìn)行界面重構(gòu).

首先確定界面法向量.界面法向量m(不一定是單位矢量)定義在網(wǎng)格中心,步長為h,則

以二維為例,則網(wǎng)格中心法向量m可由4個角點處的m值確定

界面的位置由界面法向量和同樣定義在網(wǎng)格中心的體積函數(shù)C確定,并采用Lagrangin方法跟蹤隨流動傳播的界面.

物理量定義采用MAC交錯網(wǎng)格,并采用經(jīng)典的中心差分格式.另外,由于在界面處密度、黏度和壓強(qiáng)突變,計算過程中采用多重網(wǎng)格技術(shù),每個節(jié)點處的密度和黏度采用簡單的網(wǎng)格體積平均進(jìn)行計算,壓強(qiáng)通過Possion方程進(jìn)行求解.

為了驗證模型,文獻(xiàn)[14]模擬了氣泡在水中的連續(xù)變形過程,并與Walters等[15]的試驗結(jié)果進(jìn)行了比較.

2 模擬結(jié)果分析

本文模擬了三維圓形氣泡群在靜止無限流場中重力的作用下上升和變形的過程.氣泡初始速度為零,在浮力的作用下在靜止的流場中開始上升,邊界采用無滑移條件.模擬區(qū)域為長方體,計算網(wǎng)格為66×66×128,時間步長和空間步長分別為2.5×10-6和1.95×10-4,周圍流體與氣泡密度和黏度比分別為ρf/ρb=40,μf/μb=80,Eo=ρbgd2/σ描述氣泡的尺寸,其中d為氣泡直徑.

2.1 氣泡群的上升

氣泡群在上升的過程中通常伴隨著碰撞、變形和融合的過程,雷諾數(shù)及氣泡數(shù)量對其亦有不同的影響.

圖1給出了三維球形氣泡群在上升過程中的變化情況,氣泡數(shù)量為144個,氣泡沿z上升方向排成4層,每層36個,且3個方向氣泡間距相同,氣泡在上升過程中保持球形形狀.從圖中可以看出,在無量綱時間τ=45時,氣泡由最初的整齊排列開始變化,占據(jù)中心位置的氣泡上升速度要高于周圍氣泡,這使得每層由平面變?yōu)橄蛏仙较蛲怀龅那?靠近底部的氣泡由于受到前部氣泡的阻擋只得朝著與上升方向垂直的方向運(yùn)動,整個群體呈現(xiàn)扇形體的形狀.隨著時間的推移,外圍的氣泡開始發(fā)生碰撞并且開始融合,從τ=75時刻對應(yīng)的圖中可以觀察到第1個融合后體積變大的氣泡隨著群體繼續(xù)上升,隨后大量氣泡的融合過程開始,其中包括左右相鄰和前后排列的氣泡之間的融合.在τ=111時刻氣泡的第2次融合過程亦已發(fā)生,此時氣泡的體積相差顯著,氣泡之間的上升速度差別亦增大,因此氣泡群體在上升方向距離加大,并在隨機(jī)的位置上繼續(xù)上升、碰撞、融合.這一現(xiàn)象與氣泡數(shù)量減至Nb=64個時變化過程類似,但由于氣泡數(shù)量少的緣故,氣泡融合時間開始得稍晚一些.

隨著雷諾數(shù)的增加,例如Re=74,Eo=6.9,氣泡變形顯著,為橢圓形,且其上升速度與小雷諾數(shù)時相比高出很多,變化過程也大不相同.氣泡在初期也表現(xiàn)出向兩側(cè)擴(kuò)展的扇形體,但程度不太明顯.值得注意的是,與表面相鄰的第2層氣泡幾乎整體與其發(fā)生碰撞、融合形成大氣泡,留下3排體積不等的氣泡繼續(xù)上升,此時氣泡間速度差別開始加大,但接下來的第2次融合過程與第1次類似,仍為整體融合,氣泡數(shù)量再次大量減少.

泡群運(yùn)動引起的速度分布見圖2.從圖2可以清楚地看到與少量氣泡相比,多個氣泡的同時出現(xiàn)使得氣泡間的相互作用、相互影響范圍擴(kuò)大,因此對周圍流體的擾動也大大增強(qiáng).但當(dāng)Re=74時,流體黏性降低,氣泡上升速度加快,在其影響下,外圍流體速度相對Re=0.8要大得多,氣泡對流場的擾動也進(jìn)一步增強(qiáng).

圖1 氣泡群的上升過程(Nb=144,Eo=0.36,Re=0.8)Fig.1 Rising processes of bubble clusters(Nb=144,Eo=0.36,Re=0.8)

圖2 氣泡群運(yùn)動引起的流場速度分布(Nb=144)Fig.2 Distribution of velocity fields induced by motion of bubble clusters(Nb=144)

圖3 氣泡群平均上升速度(Nb=144)Fig.3 Average rising velocities of bubble clusters(Nb=144)

與上升過程相對應(yīng)可以分析氣泡群整體上升速度的變化.圖3給出了氣泡群的平均上升速度,其中u*為速度.當(dāng)Re=0.8,Eo=0.36時,初期速度由零迅速增大至0.8左右,開始振蕩上升,此時由于前排氣泡的阻擋,氣泡群整體形狀向外擴(kuò)展成扇狀體,上升速度趨緩.隨著氣泡之間距離的減小,少量氣泡間開始發(fā)生碰撞與融合.從融合開始到第1個融合后大氣泡的生成時段,氣泡群速度經(jīng)歷了第1個融合后的速度降低與反彈[10],例如 τ=75時速度為1.15.隨著大量體積變大氣泡的出現(xiàn)及碰撞、融合的頻繁發(fā)生,一部分氣泡體積變大,上升速度加快,進(jìn)而帶動整體上升速度,表現(xiàn)為斜率變大、振幅增強(qiáng)和多個峰值的出現(xiàn).圖3(b)考察了氣泡形狀為橢圓形時平均速度的變化過程,顯然其上升速度要遠(yuǎn)遠(yuǎn)高于小Re時的對應(yīng)值.總體來說,氣泡群速度在上升到9.0左右以后變化不大,這是由于氣泡群體擴(kuò)散現(xiàn)象不明顯的原因.隨后在第1次2排氣泡的融合過程中,速度出現(xiàn)了明顯的先降后升現(xiàn)象,并在第2次的大量融合時段出現(xiàn)第2個整體速度的降低與反彈現(xiàn)象.

2.2 氣泡群引起的速度脈動

圖4給出了在融合發(fā)生之前氣泡表面的速度脈動,所跟蹤氣泡分別位于各層的中心,采樣無量綱時間間隔為0.15,采樣點數(shù)為400個.可以看出在氣泡體積增大和上升速度同時增大的 τ=0～30時段,氣泡表面x,y,z3個方向的速度u′,v′,w′脈動振幅都很大,這一現(xiàn)象在氣泡數(shù)量較多Nb=144時更加明顯.這是由于氣泡數(shù)量的增加使得相互之間作用和影響增強(qiáng)的結(jié)果.但在數(shù)量較少時不同層面之間的氣泡表面脈動程度相差較大,例如Nb=64時第2層速度脈動與位于表面第1層氣泡間的差別較明顯,并在最大速度到來之前出現(xiàn)2個脈動峰值.

圖4 氣泡表面的速度脈動(Eo=0.36,Re=0.8)Fig.4 Fluctuating velocities on bubble surface(Eo=0.36,Re=0.8)

3 結(jié) 語

本文采用VOF(volume-of-fFluid)中的PLIC(piecewise linear interface calculation)界面重構(gòu)方法模擬研究三維氣泡群的動力學(xué)特性及其對周圍流場的擾動.結(jié)果表明,在氣泡群上升過程中,伴隨著隨機(jī)的碰撞與融合等現(xiàn)象,氣泡分布經(jīng)歷了由規(guī)則到扇形體再到隨機(jī)的過程;在大雷諾數(shù)時融合過程表現(xiàn)為整體行為.與此相對應(yīng),氣泡群平均上升速度在小雷諾數(shù)時隨著時間的推移逐漸加快,并在氣泡發(fā)生融合過程中出現(xiàn)震蕩.氣泡數(shù)量的增多對氣泡本身和外圍流體的速度擾動都有所增強(qiáng).

[1]鄭源,索麗生,屈波,等.有壓輸水管道系統(tǒng)含氣水錘研究[J].河海大學(xué)學(xué)報:自然科學(xué)版,2005,33(3):277-281.(ZHENG Yuan,SUO Li-sheng,QU Bo,et al.Study onwater hammer with gas in pressuriedwater delivery system[J].Journal of Hohai University:Natural Science,2005,33(3):277-281.(in Chinese))

[2]FERRANTEA,ELGHOBASHI S E.On the accuracy of the two-fluid formulation in DNS of bubble-laden turbulent boundary layers[J].Physics of Fluids,2007,19:1-8.

[3]張淑君,吳錘結(jié),王惠民.單個三維氣泡運(yùn)動的數(shù)值模擬[J].河海大學(xué)學(xué)報:自然科學(xué)版,2005,33(5):534-537.(ZHANG Shu-jun,WU Chui-jie,WANG Hui-min.Numerical simulation of the rise and deformation of a three dimenssionalgas bubble[J].Journal of Hohai University:Natural Science,2005,33(5):534-537.(in Chinese))

[4]HETSRONI G.Handbook of multiphase systems[M].Washington:Hemisphere-McGraw Hill,1982:550-556.

[5]SMEREKA S.On the motion of bubbles in a periodic box[J].J Fluid Mech,1993,254:79-112.

[6]SMEREKA S.A vlasov equations for soundwave propagation in bubbly fluids[J].J Fluid Mech,2002,454:287-325.

[7]SANGANI A S,DIDWANIA A K.Dynamic simulationsof flow of bubbly liquids at large Reynolds numbers[J].JFluidMech,1993,250:307-337.

[8]PELET P D,BIESHEUVEI A.On the motion of gas bubbles in homogeneous isotropic turbulence[J].J FluidMech,1997,336:221-244.

[9]HU H H.Direct numerical simulations of flows of solid-liquid mixtures[J].Phys Fluids,1996,22:335-352.

[10]JOHNSON A A,TEZDUYAR T E.3D simulation of fluid-particale interactions with the number of particles reaching 100[J].Comput Meth Appl Mech Eng,1997,145:301-321.

[11]RADL S,TRYGGVASON G,KHINAST J.Flow and mass transfer of fully resolved bubbles in non-Newtonian fluids[J].AIChE Journal,2007,53:1861-1878.

[12]LU J,TRYGGVASON G.Numerical study of turbulent bubbly downflows in a vertical channel[J].Physics of Fluids,2006,18(10):103-302.

[13]LU J,BISWAS S,TRYGGVASON G.Study of laminar bubbly flows in a vertical channel[J].International Journal of Multiphase Flow,2006,32:643-660.

[14]張淑君,吳錘結(jié).氣泡之間相互作用的數(shù)值模擬[J].水動力研究與進(jìn)展:A輯,2008,23(6):681-686.(ZHANG Shu-jun,WU Chui-jie.Numerical simulation of the intera between two three-dimensional deformable bubbles[J].Chinese Journal of Hydrodynamics,2008,23(6):681-686.(in Chinese))

[15]WALTERS J K,DAVIDSION J F.The initialmotion of a gas bubble formed in an inviscid liquid,part 1:the two dimensional bubble[J].J Fluid Mech,1962,12:408-416.