三維雙氣泡融合的格子Boltzmann模擬

呂雅琪，聶德明，林建忠

（中國計量學(xué)院計量測試工程學(xué)院，浙江杭州 310018）

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

三維雙氣泡融合的格子Boltzmann模擬

呂雅琪，聶德明?，林建忠

（中國計量學(xué)院計量測試工程學(xué)院，浙江杭州 310018）

針對基于自由能模型的格子Boltzmann方法，推導(dǎo)D3Q15格子模型對應(yīng)的平衡態(tài)分布函數(shù)；采用該模型模擬雙氣泡的融合過程.結(jié)果表明，氣泡的融合不僅與它們之間的初始間距有關(guān)，還與表面張力有關(guān).表面張力越大，氣泡融合的臨界距離也越大.此外，研究氣泡融合速度與初始間距、表面張力及粘性系數(shù)的關(guān)系.

氣泡融合；格子Boltzmann方法；表面張力；氣液界面

0 引言

氣泡運動廣泛存在于能源、動力和環(huán)境工程設(shè)備與工藝中，氣泡的運動特性一直以來都是研究的熱點.與單氣泡相比，多個氣泡的運動存在相互作用，如碰撞、融合及破裂等，會出現(xiàn)豐富而復(fù)雜的運動現(xiàn)象.此外，由于氣泡在運動中可能會變形，因此比一般的顆粒兩相流問題更加復(fù)雜.

氣泡的運動屬于典型的帶有相界面的氣液兩相流動.目前已經(jīng)發(fā)展出多種數(shù)值方法用于模擬氣泡的運動，如Level Set方法［1－3］、VOF（Volume of Fluid）方法［4－6］和格子Boltzmann方法［7－9］等.這些方法在氣泡運動的研究中獲得了許多成果.格子Boltzmann方法（LBM）是20世紀(jì)80年代中期出現(xiàn)的流體計算和物理建模的一種新方法.與傳統(tǒng)數(shù)值方法不同，LBM并非以連續(xù)介質(zhì)理論為基礎(chǔ)，而是基于微觀尺度上的統(tǒng)計力學(xué)Boltzmann方程.在計算中LBM只涉及一系列的碰撞和流過程步驟，并不需要求解復(fù)雜的非線性方程，對于復(fù)雜結(jié)構(gòu)的無滑移邊界條件很容易實現(xiàn)，且能更加準(zhǔn)確合理地描述物理系統(tǒng)的微觀動力學(xué)機制.因此非常適合涉及界面運動的模擬.近二十幾年來已經(jīng)發(fā)展出幾種適用于多相流的格子Boltzmann模型，如顏色模型［10］、偽勢模型［11］和自由能模型［12］.顏色模型［10］通過顏色梯度來表示不同相之間的界面，而偽勢模型［11］則是通過偽勢函數(shù)來描述不同相之間的相互作用.然而這兩種模型在涉及熱力學(xué)方面計算時遇到很大的限制.自由能模型［12］引入對流擴散方程來描述氣液界面，該方程求解簡單，易于實現(xiàn)，且可以較好地與其他熱模型結(jié)合，因此更加適合氣泡運動的模擬和研究.

基于自由能模型，Zheng等人［13］提出一種適用于較大氣液密度比的格子Boltzmann模型，他們采用雙分布函數(shù)，分別用來描述流體運動方程和界面運動方程，并針對D2Q9格子模型給出了相應(yīng)的平衡態(tài)的分布函數(shù).通過模擬計算，他們認(rèn)為雙氣泡的融合與氣泡的初始間距有關(guān)，且給出了以下結(jié)論，即只有當(dāng)雙氣泡初始間距小于兩倍界面厚度時，這兩個氣泡才能融合在一起并形成一個更大的氣泡.但這一結(jié)論并不全面，首先氣泡之間的融合條件應(yīng)該還涉及其它參數(shù)如表面張力和流體粘性等的影響.其次這是基于二維的模擬結(jié)果，三維氣泡的融合過程可能存在較大差異，需要進一步的研究.另外，從物理機制上來說，氣泡的融合現(xiàn)象在多氣泡運動的情況中可能會頻繁地發(fā)生，當(dāng)氣泡相互融合成更大的氣泡時，會顯著影響氣泡群的運動特性，進而影響到氣泡群的運動速度以及平均阻力系數(shù)等［15］，因此值得深入的研究.然而，以往對于氣泡融合條件的研究較缺乏，尤其在三維方面更缺乏系統(tǒng)的研究成果.本文采用Zheng等人［13－14］提出的基于自由能模型的格子Boltzmann方法，推導(dǎo)D3Q15格子模型對應(yīng)的平衡態(tài)分布函數(shù)，并對雙氣泡的融合條件進行模擬研究，考慮表面張力、粘性系數(shù)以及氣泡間距的影響.此外，論文還對氣泡的融合速度進行細(xì)致地研究.

1 數(shù)值方法

流體運動的Navier?Stokes方程和界面演化方程為

式中θM是遷移率，P是壓力，F(xiàn)b是體積力，n是半密度和，?是半密度差（也叫相序參數(shù)），即n＝（ρA＋ρB）/2，?＝（ρA－ρB）/2.ρA和ρB是流體A和流體B的密度.μ?是化學(xué)勢能，根據(jù)Zheng等人［13］，μ?可以表示為

其中，??是與平衡狀態(tài)相關(guān)的常數(shù)，即

ρH和ρL分別代表兩相中的高密度和低密度.系數(shù)A定義為，

式中σ為表面張力，W為界面厚度.κ的表達式為

1.1 連續(xù)方程和動量方程的求解

根據(jù)Zheng等人［13］，連續(xù)性方程（1）和動量方程（2）采用帶有外力項的格子Boltzmann演化方程求解，即

為了模擬三維的氣泡運動，本文采用D3Q15格子模型，其離散速度方向為

采用如下所示的平衡態(tài)分布函數(shù)，

為了得到式中的系數(shù)Ai和權(quán)系wi，引入質(zhì)量、動量和二階速度矩的守恒關(guān)系，即，

其中格子速度cs＝c2/3，將平衡態(tài)分布函數(shù) （10）代入式（11）、（12）及（13），考慮D3Q15的格子速度矢量（9），可以得到對應(yīng)的系數(shù)Ai和權(quán)系數(shù)wi，

1.2 界面捕捉方程的求解

為了求解界面方程（3），采用如下形式的格子Boltzmann演化方程［13］，

式中，碰撞項為Ωi＝［g（i0）（x，t）－gi（x，t）］/τ?，t?為無量綱松弛時間，q與t?有關(guān)即q＝1/（τ?＋0.5）．

根據(jù)Zheng等人［14］，gi的平衡態(tài)分布函數(shù)為

采用D3Q7模型，其離散的速度方向為，

為了得到式中的系數(shù)Ai、Bi和Ci，引入質(zhì)量、動量和二階速度矩的守恒關(guān)系，即，

則可以推得各項系數(shù)

其中Γ與遷移系數(shù)相關(guān).

1.3 模型驗證

根據(jù)Jacqmin［16］，關(guān)于球形氣泡的相序參數(shù)即氣液界面的理論值為，

其中xc、yc和zc是氣泡中心的坐標(biāo).上式可以用來表示氣液的界面曲線.為了驗證本文的方法，采用式（20）進行驗證.

設(shè)置計算區(qū)域為N3＝64×64×64，表面張力σ＝1.0，界面寬度W＝3.5.對以下兩種情況進行模擬：①氣泡半徑R＝10，ρH＝1 000，ρL＝1；②氣泡半徑R＝15，ρH＝500，ρL＝1.結(jié)果如圖1所示，符號表示的是本文D3Q15對應(yīng)的數(shù)值解，實線和虛線對應(yīng)的是理論解即式（20）.可以看到，本文采用的方法準(zhǔn)確地捕捉到了氣液的界面曲線.

Zheng等人［13］研究表明，界面厚度W對模擬的結(jié)果有影響，具體而言，對于半徑為R的氣泡，當(dāng)選取的W較小時會出現(xiàn)比較明顯的誤差.為此，我們設(shè)置 R＝20，ρH＝1 000，ρL＝1以及σ＝1.0，分別選取W＝1.5，2.0，2.5，3.5 和4.5，通過計算可以得到氣泡內(nèi)外的壓力差Δp.對于球形氣泡，通過內(nèi)外壓力差和表面張力的平衡方程可以得到Δp ＝2σ/R（二維的表達式為Δp＝σ/R）.根據(jù)這一表達式可以計算表面張力，顯然這個結(jié)果與設(shè)置的值有差別，通過這個差別可以評估相應(yīng)的誤差.結(jié)果如圖2所示，可以看到，隨著界面厚度W的增大，計算所得的表面張力逐漸靠近理論值即σ＝1，當(dāng)W＝4.5時，計算值與理論值的誤差小于5％，因此在后續(xù)的計算中界面厚度均采用W＝4.5.

圖1 氣液界面的數(shù)值解與理論解Fig.1 Analytical solution and numerical results for gas?liquid interface profile

另外，Huang等人［14］曾建立基于D3Q19的格子Boltzmann方法并用來模擬氣液兩相流動.與D3Q19速度離散模型相比，本文基于D3Q15的方法在每個節(jié)點節(jié)省了4個方向的存儲數(shù)據(jù).對于三維計算而言，網(wǎng)格量大，因此既可以節(jié)省內(nèi)存，同時也提高計算效率.為了說明這一點，設(shè)置以下兩個算例：①計算區(qū)域N3＝60×60×60，氣泡半徑R＝10，表面張力σ＝0.5，初始?xì)馀蓍g距d＝0.5W；②計算區(qū)域N3＝120×120×120，氣泡半徑R＝20，表面張力σ＝1.0，初始?xì)馀蓍g距d＝0.1W.為了得到氣泡的運動的速度，需要采用統(tǒng)計平均的方法進行計算，

式中?＜0代表的是氣相區(qū)域.由于兩個氣泡在流場中是水平對稱的，因此根據(jù)（21）得到的氣泡平均速度始終為零.為了更加細(xì)致地研究氣泡的融合過程，僅對計算區(qū)域的左半部分的流場進行統(tǒng)計，這樣得到的是左邊的氣泡在運動和融合過程中的速度.結(jié)果如圖3所示，從圖中可以看出，采用D3Q15和D3Q19速度離散模型計算出的氣泡融合速度曲線幾乎重合，說明兩個模型的計算結(jié)果是一致的.更進一步地，表1給出了二者對應(yīng)的計算時間，當(dāng)計算區(qū)域為N3＝60×60×60，達到20 200迭代步時，D3Q15模型比D3Q19模型所用的時間減少了約18％；當(dāng)計算區(qū)域為N3＝120×120×120，達到22 500迭代步時，D3Q15速度離散模型比D3Q19速度離散模型所用的時間減少了約16％.因此，本文采用的D3Q15模型能較顯著地節(jié)省計算資源和提高運算效率.

圖2 通過界面厚度W計算得到的表面張力σFig.2 Surface tension computed at different thickness of interface layer

圖3 D3Q15和D3Q19格子模型的氣泡融合速度Fig.3Merging velocity in D3Q15 and D3Q19

2 結(jié)果及討論

首先研究氣泡初始間距對氣泡融合的影響.計算區(qū)域設(shè)置為N3＝120×120×120，采用周期性邊界條件，初始時刻流場處于靜止?fàn)顟B(tài).參數(shù)設(shè)置為氣泡半徑R＝20，無量綱時間因子τn＝0.6和τ?＝0.7，ρH＝1 000，ρL＝1，表面張力σ＝0.5，Γ＝800以及界面厚度W＝4.5.

表1 D3Q15和D3Q19速度離散模型效率Table 1 Com putational efficiency in D3Q15 and D3Q19 velocity model

圖4給出氣泡初始間距d＝2.0W時雙氣泡融合的過程.當(dāng)兩氣泡開始接觸時氣泡的界面發(fā)生形變，在表面張力的作用下，兩氣泡有融合成一個氣泡的趨勢，如圖4（a）和（b）所示.氣泡呈橫向橢球形時氣泡的變形最大，如圖4（c）所示，此時在表面張力作用下氣泡有恢復(fù)成球形的趨勢，如圖4（d）所示.當(dāng)氣泡恢復(fù)成球形后，在慣性的作用下，氣泡繼續(xù)朝向中心運動.此時圓形逐漸變成縱向橢球形，如圖4（e）所示.經(jīng)過幾次振蕩變形，伴隨著粘性的耗散，氣泡最終變成穩(wěn)定的球形，如圖4（f）所示.

圖4 兩個氣泡的融合過程（d＝2.0W）Fig.4 Merging process shown by contours of order parameter at different times（d＝2.0W）

若增大氣泡的初始間距如d＝2.5W，則發(fā)現(xiàn)在計算中即使經(jīng)歷很長時間兩個氣泡也沒有相互靠近的趨勢，一直保持靜止不動，不發(fā)生融合，如圖5所示.顯然，氣泡的初始間距是氣泡融合的關(guān)鍵參數(shù)之一.

圖5 兩個氣泡相序參數(shù)隨時間演化（d＝2.5W）Fig.5 Merging process shown by contours of order parameter at different times（d＝2.5W）

通過計算發(fā)現(xiàn)，氣泡的融合不僅與初始間距的大小有關(guān)，還與表面張力有關(guān).為了研究這二者對氣泡的影響，選取三組表面張力即σ＝0.5，1.0和1.5，在不同間距d下進行模擬，觀察氣泡是否能靠近并融合在一起.統(tǒng)計這些結(jié)果可以得到如圖6所示的圖形，可以看到，表面張力越大，氣泡能融合的最大間距及臨界間距dc越大.例如，當(dāng)表面張力σ＝0.5時，兩個氣泡能融合在一起的臨界間距為dc＝2.0W，而當(dāng)σ＝1.0時，臨界間距為dc＝2.5W.顯然，這與Zheng等人［13］在二維情況下研究的結(jié)果不一致.說明圓形氣泡間的相互作用與球形氣泡之間的相互作用并不等同.特別地，當(dāng)表面張力σ＝1.5時，即使在氣泡間距為4W時還能觀察到融合的現(xiàn)象，說明表面張力對氣泡融合的影響非常顯著.究其原因，表面張力σ越大，化學(xué)勢能越大，因此氣泡融合的驅(qū)動力越大.

圖6 氣泡融合與間距和表面張力的關(guān)系Fig.6 Bubblemerging dependence on surface tension and bubble distance

圖7給出不同表面張力對應(yīng)的氣泡融合速度曲線.可以看到，融合過程的速度曲線相似，即經(jīng)歷一正一負(fù)的峰值之后迅速衰減為零.對于同一表面張力，不同的初始間距對氣泡的融合速度有一定的影響，如圖所示，當(dāng)氣泡間距小于1.5W時，速度峰值明顯更大，而當(dāng)間距大于1.5W時，可以看到，所有的速度曲線的幾乎一樣.此外，間距越小氣泡應(yīng)該越早融合，這一點從圖7（a）和（b）可以得到證實，但比較特殊的是，當(dāng)表面張力σ＝1.5時，發(fā)現(xiàn)d＝5W較之d＝2.5W，3.5W和4W更早融合，這與計算中設(shè)置的周期性邊界條件有關(guān).圖8給出了相同間距的情況下不同表面張力對氣泡融合速度的影響，可以看到，表面張力越大，氣泡開始融合的時間越早，而且氣泡融合的速度峰值也越大.

圖7 不同初始間距對應(yīng)的氣泡融合速度Fig.7 Merging velocity at different initial bubble distances

下面研究流體粘性系數(shù)對氣泡融合過程的影響，固定氣泡初始間距 d＝1.5W，分別取無量綱時間因子τn＝0.55，0.57和0.60，則粘性系數(shù)ν分別為1/60，1.4/60 和1/30，對應(yīng)的氣泡融合速度曲線如圖9所示.可以看到，粘性系數(shù)越小氣泡融合速度的峰值越大，且速度曲線的振蕩越明顯.ν＝1/60時出現(xiàn)了5個速度峰值，ν＝1.4/60時出現(xiàn)了4個速度峰值，而ν＝1/30時只出現(xiàn)了3個速度峰值.這說明氣泡融合所具有的化學(xué)勢能是通過粘性進行耗散的，所以粘性越大，耗散越快.更進一步，圖10給出了兩個氣泡融合成一個氣泡的過程中，氣泡呈縱向橢球形時變形最大的形態(tài).可以看到，隨著粘性系數(shù)的增加，氣泡最大變形越來越小.這是因為粘性系數(shù)體現(xiàn)了流體對氣泡運動阻礙的強烈程度.粘性系數(shù)大，氣泡受到的阻力大，運動過程中的粘性耗散也大.所以粘性系數(shù)越大，氣泡融合的速度峰值越小，氣泡形變越小.

圖8 不同表面張力對應(yīng)的氣泡融合速度（d＝2.0W）Fig.8 Merging velocity with different surface tension coefficients（d＝2.0W）

圖9 不同粘性系數(shù)對應(yīng)的氣泡融合速度Fig.9 Merging velocity with different viscosities

圖10 不同粘性系數(shù)下氣泡縱向最大變形形態(tài)Fig.10 Contours of order parameter for bubble ofmaximum vertical deformation with different viscosities

從圖9和圖10可以看到，流體的粘性顯然對氣泡的融合過程具有影響，但應(yīng)該與氣泡是否融合無關(guān).為了說明這一點同時為了節(jié)省計算時間，我們設(shè)置計算區(qū)域為N3＝60×60×60，氣泡半徑R＝10，ρH＝1 000，ρL＝1，表面張力σ＝0.5以及界面厚度W＝2.25，分別取τn＝0.55，0.57和0.60.考察不同粘性下氣泡的融合過程.通過計算發(fā)現(xiàn)，對于上述三種情況，當(dāng)氣泡間距小于兩倍界面厚度時都發(fā)現(xiàn)氣泡融合現(xiàn)象，否則都不會融合.由此可見，粘性的大小不影響氣泡是否融合.

3 結(jié)論

采用基于自由能模型的Boltzmann方法，模擬了三維情況下雙氣泡的融合.首先針對D3Q15格子模型推導(dǎo)了對應(yīng)的平衡態(tài)分布函數(shù)，并通過氣液界面的理論解進行了驗證.計算結(jié)果表明，氣泡的融合不僅與氣泡之間的初始間距有關(guān)，還與表面張力有關(guān)，結(jié)論如下：

1）與D3Q19速度離散模型相比，D3Q15速度離散模型在每個節(jié)點節(jié)省了4個方向的存儲數(shù)據(jù)，且所用的計算時間減少15％以上.因此D3Q15速度離散模型能較顯著地節(jié)省計算資源和提高運算效率.

2）表面張力越大，雙氣泡融合成一個氣泡的臨界間距越大.特別地，當(dāng)表面張力σ＝1.5時，即使間距d＝4.0W時也能觀察到氣泡的融合，這與二維的結(jié)論不一致.另外，表面張力越大，氣泡融合的時間越早，融合的速度也越快.

3）氣泡初始間距越小，氣泡融合的時間越早.而且，當(dāng)氣泡間距d＜1.5W時，間距越小氣泡融合速度的峰值越大，而當(dāng)氣泡間距d＞1.5W時，氣泡間距對融合速度的峰值幾乎沒有影響.

4）流體黏性系數(shù)對氣泡能否融合不起關(guān)鍵作用，但流體黏性會影響氣泡的融合過程.流體粘性系數(shù)越小，氣泡融合的速度峰值越大，氣泡變形越大，融合速度曲線的振蕩越明顯，從定量的角度描述了粘性耗散的過程.

［1］ 李彥鵬，張乾隆，白博峰.豎直通道內(nèi)相鄰氣泡對上升的直接數(shù)值模擬［J］.熱能動力工程，2007，22（4）：375－379.

［2］ 陳菲，張夢萍，徐勝利.運動激波和氣泡串相互作用的初步數(shù)值模擬［J］.計算物理，2004，21（5）：443－448.

［3］ Sussman M，F(xiàn)atemi E，Smereka P，Osher S.An improved level setmethod for incompressible two?phase flows［J］.Comput Fluids，1998，27：663－680.

［4］ Pilliod JE，Puckett EG.Second?order volume?of?fluid algorithms for trackingmaterial interfaces［J］.JComput Phys，2004，199：465－502.

［5］ 李維仲，趙大勇，陳貴軍.豎直流道寬度對氣泡運動行為影響的數(shù)值模擬［J］.計算力學(xué)學(xué)報，2006，23（2）：196－201.

［6］ 張淑君，吳錘結(jié).氣泡之間相互作用的數(shù)值模擬［J］.水動力學(xué)研究與進展A輯，2008，23（6）：681－686.

［7］ 孫濤，李維仲，楊柏丞，祝普慶.氣泡群上升過程相互作用的格子Boltzmann三維數(shù)值模擬［J］.化工學(xué)報，2013，（5）：1586－1591.

［8］ 曾建邦，李隆鍵，蔣方明.氣泡成核過程的格子Boltzmann方法模擬［J］.物理學(xué)報，2013，62（17）：176401.

［9］ Cheng M，Hua J，Lou J.Simulation of bubble?bubble interaction using a lattice Boltzmann method［J］.Comput Fluids，2010，39：260－270.

［10］ Gunstensen A K，Rothman D H，ZaleskiS，ZanettiG.Lattice Boltzmannmodel of immiscible fluids［J］.Phys Rev A，1991，43：4320－4327.

［11］ Shan X，Chen H.Lattice Boltzmannmodel for simulating flowswithmultiple phases and components［J］.Phys Rev E，1993，47（3）：2557－2562.

［12］ SwiftM R，OsborneW，Yeomans JM.Lattice Boltzmann simulation of nonideal fluids［J］.Phys Rev Lett，1995，75：830－833.

［13］ Zheng H W，Shu C，Chew Y T.A lattice Boltzmann model formultiphase flows with large density ratio［J］.JComput Phys，2006，218：353－371.

［14］ Huang H B，Zheng H W，Lu X Y，Shu C.An evaluation of a 3D free?energy?based lattice Boltzmann model for multiphase flowswith large density ratio［J］.Int JNumer Meth Fluids，2010，63：1193－1207.

［15］ Gupta A，Kumar R.Lattice Boltzmann simulation to studymultiple bubble dynamics［J］.Int JHeat Mass Trans，2008，51：5192－5203.

［16］ Jacqmin D.Calculation of two?phase Navier?Stokes flows using phase?fieldmodeling［J］.JComput Phys，1999，155：96－127.

Lattice Boltzmann Simulation of Gas Bubble M erging in Three Dimensions

LV Yaqi，NIE Deming，LIN Jianzhong
（College ofMetrology and Technology Engineering，China Jiliang University，Hangzhou 310018，China）

Equilibrium distribution functions for latticemodel D3Q15 are proposed in the framework of lattice Boltzmann method in free energymodel.Bubblemerging process is studied.It shows that bubblemerging not only depends on initial bubble distance，but also depends on surface tension.The greater the surface tension is，the greater the critical distance for bubblemerging.Furthermore，influence of initial bubble distance，surface tension and viscosity onmerging velocity are investigated.

bubblemerging；lattice Boltzmann method；surface tension；gas?liquid interface

O359

A

2014－09－26；

2014－12－21

國家重點基礎(chǔ)研究發(fā)展計劃（2011CB706501）及國家自然科學(xué)基金（11272302）資助項目

呂雅琪（1991－），女，碩士生，主要從事氣液兩相流的研究，E?mail：lv_yaqi＠126.com

?通訊作者：E?mail：nieinhz＠cjlu.edu.cn

Received date： 2014－09－26；Revised date： 2014－12－21