田　輝，房　媛，王文成 ，鄒克武，葉陽輝

(1.承德石油高等?？茖W(xué)校 機(jī)械工程系，承德　067000；2.承德石油高等專科學(xué)校 學(xué)生處，承德　067000；3.西安交通大學(xué) 能源與動力工程學(xué)院，西安　710049)

?

基于Level Set方法的水中氣泡上升過程數(shù)值模擬

田輝1，房媛2，王文成1，鄒克武1，葉陽輝3

(1.承德石油高等專科學(xué)校 機(jī)械工程系，承德067000；2.承德石油高等?？茖W(xué)校 學(xué)生處，承德067000；3.西安交通大學(xué) 能源與動力工程學(xué)院，西安710049)

摘要：基于Level Set方法發(fā)展了一種高分辨率的求解氣液相界面遷移特性的數(shù)值方法。通過三步Runge-Kutta Crank-Nicholson投影方法求解流場，Level Set方法捕捉氣液交界面位置。在驗(yàn)證算法的捕捉性能的基礎(chǔ)上，研究了氣泡在水中上升過程的變形以及運(yùn)動特性。分析結(jié)果顯示氣泡上升過程中的變形是氣泡頂部/底部壓差、兩側(cè)渦及表面張力綜合作用的結(jié)果；氣泡在接近壁面時(shí)由于兩側(cè)不對稱渦及壁面排擠作用的影響，上升過程不斷向中心區(qū)域靠近并伴隨有逆時(shí)針旋轉(zhuǎn)；氣泡初始位置距壁面越近，壁面的影響越顯著。

關(guān)鍵詞：Level Set;投影法;Zalesak問題;氣泡上升

氣泡運(yùn)動大量存在于能源、動力、石化、航空航天等各個(gè)領(lǐng)域，氣液相界面遷移對流體系統(tǒng)的宏觀流動、傳熱及傳質(zhì)特性有著顯著的影響。掌握氣泡的運(yùn)動規(guī)律和變形機(jī)理，可為設(shè)計(jì)高效節(jié)能并可以持續(xù)穩(wěn)定運(yùn)行的兩相流設(shè)備提供理論依據(jù)[1]。為此國內(nèi)外學(xué)者發(fā)展了各種數(shù)值模擬方法用于研究多相流相界面遷移特性，如當(dāng)前被普遍應(yīng)用的：Front Tracking Method[2]、VOF[3]、Level Set[4]等系列算法。

Front Tracking Method通過拉格朗日思想跟蹤表征界面位置的虛擬標(biāo)記粒子來獲取相界面的運(yùn)動規(guī)律。其計(jì)算精度較高，一些典型算例已將此方法的求解結(jié)果作為基準(zhǔn)解。但此方法會消耗大量的計(jì)算資源，并且在求解拓?fù)浣Y(jié)構(gòu)發(fā)生變化的問題時(shí)虛擬粒子的增減機(jī)制尚需完善。因而也制約了Front Tracking Method在工程問題上的應(yīng)用。VOF方法通過歐拉思想定義計(jì)算網(wǎng)格內(nèi)某相的體積含率來表征流場內(nèi)相的分布，體積含率F∈(0,1)的區(qū)域即為界面位置。VOF方法可以很好的保證系統(tǒng)的質(zhì)量守恒，但不能準(zhǔn)確地求解當(dāng)?shù)厍?，進(jìn)而不能準(zhǔn)確描述表面張力的影響。同時(shí)VOF方法還必須根據(jù)特定的界面重構(gòu)方法來獲得界面的位置，限制了界面的捕捉精度。同樣基于歐拉思想的Level Set方法通過定義計(jì)算域內(nèi)各點(diǎn)到界面的距離為相函數(shù)，定義距離函數(shù)φ=0描述相界面。由于距離函數(shù)為連續(xù)函數(shù)便于進(jìn)行微分、求導(dǎo)等操作，所以可以較準(zhǔn)確地考慮表面張力對相界面遷移的影響，同時(shí)可以獲得更銳利的相界面分布。但Level Set方法隨著計(jì)算的進(jìn)行存在質(zhì)量虧損問題，Sussman[5]提出了一種距離函數(shù)重新初始化方法在獲得銳利的界面分布的同時(shí)有效地保證系統(tǒng)質(zhì)量守恒。陳斌[6-8]基于Front Tracking Method研究了壁面對氣泡上升的影響，計(jì)算結(jié)果顯示壁面對氣泡的作用主要體現(xiàn)在對氣泡周圍流場的抑制上，由此而形成的非對稱流場導(dǎo)致氣泡逐漸偏離壁面。然而文中并未給出氣泡由靜止達(dá)到穩(wěn)定狀態(tài)的過程內(nèi)氣泡變形與其運(yùn)動軌跡的相互影響。

本文通過自編程序?qū)崿F(xiàn)了基于投影法及Level Set方法對氣液兩相界面遷移特性的捕捉。在介紹數(shù)值方法的基礎(chǔ)上，通過求解Zalesak問題驗(yàn)證算法在求解曲率變化較大問題中的準(zhǔn)確性、可行性。在此基礎(chǔ)上對水中氣泡自由上升過程的變形及運(yùn)動規(guī)律進(jìn)行了數(shù)值研究。

1數(shù)值方法

1.1控制方程

描述不可壓縮氣液兩相流問題的無量綱Navier-Stokes方程可表示為：

(1)

·U=0

(2)

通過Level Set方法捕捉氣液交界面。相函數(shù)φ被定義為計(jì)算域內(nèi)各點(diǎn)到界面的距離。φ>0表示液體區(qū)域，φ<0表示氣體區(qū)域，而φ=0給出了氣液相界面的位置。Osher[4]等推導(dǎo)了速度場U作用下函數(shù)φ隨時(shí)間變化所滿足的相函數(shù)附加方程(3)，為Level Set方法的提出奠定了基礎(chǔ)。

(3)

采用文獻(xiàn)[5]中Sussman提出的距離函數(shù)重新初始化方法保證系統(tǒng)質(zhì)量守恒。重新初始化方程表示為：

φt=signε(φ0)(1-|φ|)

(4)

(5)

(6)

1.2求解方法

采用三步Runge-Kutta Crank-Nicholson投影方法求解Navier-Stokes方程，在時(shí)間上保證二階精度?？臻g上擴(kuò)散相采用二階中心差分格式離散，對流相采用二階迎風(fēng)格式離散。為增加計(jì)算穩(wěn)定性采用積分平均法求解相函數(shù)附加方程(3)，五階WENO格式離散相函數(shù)重新初始化方程(4)。離散方程可表示為：

(7)

(8)

(9)

(10)

(11)

φm+1=φm-Δt(αm+1Um·φm+βm+1Um-1·φm-1)

(12)

其中：

(13)

(14)

2算法驗(yàn)證

以Zalesak[9]問題為測試算例對本文發(fā)展的界面捕捉算法進(jìn)行有效性驗(yàn)證。計(jì)算域?yàn)?×1的正方形區(qū)域，四個(gè)邊界對應(yīng)為周期性邊界條件。計(jì)算域被一開了矩形槽的圓劃分為內(nèi)、外兩部分。圓心位于(0.5，0.75)，半徑為0.15，矩形槽長0.25，寬0.05，如圖1所示。給定速度場為：

(15)

本文沿用參考文獻(xiàn)[9]中對計(jì)算誤差的定義：

(16)

其中：φe(t)、φc(t)分別為t時(shí)刻精確解及數(shù)值解，L為界面長度。

t=6.28時(shí)刻帶缺口的圓形沿計(jì)算域中心旋轉(zhuǎn)完一周。表1給出了本文算法與文獻(xiàn)[9]中結(jié)果的誤差比較。表中結(jié)果顯示本文發(fā)展算法在較粗網(wǎng)格內(nèi)計(jì)算精度稍差，隨著網(wǎng)格加密計(jì)算精度與ELVIRA方法相當(dāng)。圖2給出了四個(gè)典型時(shí)刻的數(shù)值解及初始位置精確解?？梢娦D(zhuǎn)一周后數(shù)值解除在尖角部分與精確解有所偏差外，其余部分吻合良好。

表1　測試問題計(jì)算誤差E(t=6.28)

網(wǎng)格CLSVOFELVIRA本文算法100×1000.005720.005670.00591200×2000.002520.002620.00259400×4000.001060.001210.00117

3氣泡上升過程的數(shù)值模擬

本部分以靜止水中上升的三維氣泡為研究對象。記氣泡初始半徑為R0，計(jì)算域?yàn)?1R0×11R0×22R0。計(jì)算域上邊界設(shè)置自由面邊界條件，?w/?z=0；其余五個(gè)邊界設(shè)為固壁，邊界無滑移條件。計(jì)算采用多重網(wǎng)格技術(shù)，四重均分網(wǎng)格，網(wǎng)格數(shù)為41×41×81。

3.1上升氣泡的穩(wěn)定過程

圖3(a)中每隔單位時(shí)間給出氣泡的變形情況，展示了氣泡在水中自由上升的過程?？梢钥闯鲎畛蹼A段上升速度、變形均較小，由于慣性造成氣泡頂部與底部速度差異使得氣泡有被拉長的趨勢。但隨著氣泡與周圍水的相對運(yùn)動加劇，在氣泡兩側(cè)各形成一個(gè)渦。其表現(xiàn)為推動氣泡底部加速上升并使氣泡不斷向兩側(cè)伸展。當(dāng)上升到一定階段，氣泡頂部與底部的速度差、兩側(cè)渦的強(qiáng)度對氣泡的作用與氣泡表面張力作用相平衡時(shí)，氣泡形狀將不再變化保持穩(wěn)定上升。圖3(b)給出了氣泡達(dá)到穩(wěn)定時(shí)兩側(cè)渦的分布情況。由于周圍水的慣性較大，使氣泡在上升的啟動過程以及達(dá)到穩(wěn)定上升前受到較大的阻力。當(dāng)浮力與阻力達(dá)到平衡時(shí)氣泡便以恒定的速度上升。上升過程中氣泡平均速度變化情況如圖4所示。

3.2壁面對氣泡運(yùn)動的影響

通過考察初始位置距壁面不同距離的氣泡上升過程來分析壁面對氣泡上升過程的影響。圖5中所示分別為初始位置距壁面3.5R0及2R0的氣泡自由上升過程。圖中可見氣泡在上升過程不斷向中心區(qū)域靠近。并在此過程中伴隨有逆時(shí)針旋轉(zhuǎn)，且氣泡初始位置越靠近壁面旋轉(zhuǎn)越明顯。

如圖6(a)所示，上升初期氣泡左側(cè)(遠(yuǎn)離壁面的一側(cè))已經(jīng)形成了較明顯的渦，而右側(cè)渦強(qiáng)度較小。受左側(cè)逆時(shí)針旋轉(zhuǎn)渦的影響氣泡則呈現(xiàn)出明顯的旋轉(zhuǎn)特性。同時(shí)在壁面排擠作用影響下氣泡逐漸向區(qū)域中心運(yùn)動。如圖6(b)所示，隨著氣泡的上升右側(cè)的渦并沒有發(fā)展成熟，而表現(xiàn)為明顯的壁面排擠作用。使氣泡及左側(cè)的渦向區(qū)域中心運(yùn)動。此時(shí)由于氣泡與左側(cè)渦中心相對位置的變化，旋轉(zhuǎn)作用有所減弱。

4結(jié)論

本文發(fā)展了一種基于Level Set Method求解氣液相界面遷移特性的數(shù)值方法。通過求解Zalesak問題驗(yàn)證了算法在捕捉曲率變化較大問題中的準(zhǔn)確性、可行性。在此基礎(chǔ)上以水中的氣泡為研究對象，模擬了氣泡初始位置在計(jì)算域中央以及靠近壁面情況下的上升過程。結(jié)果表明：1)氣泡受浮力驅(qū)動向上運(yùn)動，在其頂部與底部間的壓差、兩側(cè)渦及表面張力共同作用下發(fā)生變形，當(dāng)以上各種作用達(dá)到平衡時(shí)氣泡達(dá)到穩(wěn)定上升狀態(tài)；2)靠近壁面的氣泡受到其兩側(cè)不等強(qiáng)度渦的影響及壁面排擠作用的影響，在上升的過程中伴隨著逆時(shí)針的旋轉(zhuǎn)及向中心區(qū)域的運(yùn)動。氣泡初始位置距壁面越近，壁面的影響越顯著。

參考文獻(xiàn)：

[1]田輝,房媛,王文成，等. 離心泵內(nèi)變工況流動特性的數(shù)值研究[J]. 承德石油高等?？茖W(xué)校學(xué)報(bào), 2015(2):28-31.

[2]Tryggvason, G. A front-tracking method for the computations of multiphase flow[J]. Journal of Computational Physics, 2001. 169(2): 708-759.

[3]Rider, W.J. and D.B. Kothe. Reconstructing volume tracking[J]. Journal of Computational Physics, 1998. 141(2): 112-152.

[4]Osher, S. and J.A. Sethian. Fronts Propagating with Curvature-Dependent Speed-Algorithms Based on Hamilton-Jacobi Formulations[J]. Journal of Computational Physics, 1988. 79(1): 12-49.

[5]Sussman, P. Smereka, S. Osher. A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[J]. Journal of Computational Physics, 1994. 114(1): 146-159.

[6]陳斌, T.Kawamura, Y.Kodama. 靜止水中單個(gè)上升氣泡的直接數(shù)值模擬[J]. 工程熱物理學(xué)報(bào), 2005(6): 82-84.

[7]陳斌. 傾斜壁面附近上升氣泡的直接數(shù)值模擬[J]. 工程熱物理學(xué)報(bào), 2007(6): 965-967.

[8]陳斌. 高粘度流體中上升氣泡的直接數(shù)值模擬[J]. 工程熱物理學(xué)報(bào), 2006,02: 255-258.

[9]Sussman, E.G. Puckett. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J]. Journal of Computational Physics, 2000. 162(2): 301-337.

Numerical Simulation of Bubble Rising in Water Based on Level Set Method

TIAN Hui1, FANG Yuan2, WANG Wen-cheng1, ZOU Ke-wu1, YE Yang-hui3

(1.Department of Mechanical Engineering, Chengde Petroleum College,Chengde 067000, Hebei, China;2.Department of Student Affairs, Chengde Petroleum College, Chengde 067000, Hebei, China;3.School of Energy and Power Engineering, Xi’an Jiaotong University,Xi’an 710049, Shaanxi, China)

Abstract:Based on level set method, a high resolution numerical method was developed to study the migration characteristics of gas-liquid interface. Three step Runge-Kutta Crank-Nicholson projection method was used to solve the Navier-Stokes equations. Level set method was applied to capturing the phase interface between air and water. Zalesak’s problem was solved to validate the numerical method. The numerical results of bubble rising in water show that the deformation of bubble during the rising process is subject to (1) the pressure difference between the top and bottom of bubble, (2) the vortex, and (3) surface tension. When close to the wall, due to the wall repulsion and the asymmetric vortex, the deformed bubble moves to the centerline with clockwise rotation. The closer to the wall, the wall effect on the bubble will be more significant.

Key words:level set; projection method; Zalesak’s problem; bubble rising

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(兩相非定常空化空泡時(shí)空演化機(jī)理及其對高效離心泵水動力性能影響的研究)：51076126

收稿日期：2015-10-30

作者簡介：田輝(1982-)，男，河北承德人，承德石油高等?？茖W(xué)校機(jī)械工程系講師，博士，主要從事氣液兩相界面遷移特性數(shù)值研究。

中圖分類號：O359

文獻(xiàn)標(biāo)識碼：A

文章編號：1008-9446(2016)02-0028-05