徐玲君，薛 陽，張 松

（1.中國水電顧問集團華東勘測設(shè)計研究院，杭州 310014；2.四川水利職業(yè)技術(shù)學(xué)院，成都 611231）

0 引 言

氣泡廣泛存在于自然界和工業(yè)過程中，由于其實際應(yīng)用的重要性，氣泡形成的研究已經(jīng)引起了很多科學(xué)家關(guān)注，并且在上個世紀(jì)就開始進行了大量的研究。氣泡在銳孔上生成和脫離時，伴隨著形狀變化和體積變化，并且時刻影響著周圍的流場變化。

在過去的幾十年中，形成了很多種模擬氣泡自由上升的方法，如VOF 法［1］，Level set法［2］，格子玻爾茲曼法［3］，但是，這些方法很少用來模擬氣泡的形成，或者都采用單一模擬方法進行研究，Gerlach［4］研究了水中銳孔處準(zhǔn)靜態(tài)氣泡的形成過程，并且闡述了數(shù)值模擬中VOF 法和表面張力模型結(jié)合的優(yōu)點。Higuera［5］通過Level set法模擬了在高粘性液體中氣泡的形成和在上升過程中氣泡的合并。Zhang［6］對氣泡形狀進行VOF 數(shù)值模擬后與實驗數(shù)據(jù)進行了對比，二者形狀吻合。Martin［7］給出了一個計算模型，這個模型對氣泡在形成過程中脫離銳孔時可能受到的影響因素給予了分類說明。Davidson，Jamialahmadi和Kumar等［8-10］通過求解氣泡受力平衡方程式來預(yù)測氣泡形成和脫離時刻的體積。梁剛濤等［11］采用CLSVOF 法對單液滴撞擊平面液膜初期的飛濺過程進行了數(shù)值模擬，LIU［12］等通過CLS VOF法對低裝載液箱晃蕩抨擊荷載進行了研究，獲得了滿意效果。廖斌等［13］運用CLSVOF 方法數(shù)值模擬了三維單氣泡在液體中的上升和變形過程，分別考察了液體的表面張力和浮力對上升氣泡變形的影響。宋云超等［14］根據(jù)CLSVOF 方法的基本思想，將VOF和Level Set兩種方法結(jié)合起來，用于追蹤不可壓縮兩相流界面的計算，獲得了更加精確的結(jié)果。

采用CLSVOF 法，引入表面張力模型［15］，并指定不同計算條件下的壁面粘附角對銳孔氣泡的形成與脫離進行了數(shù)值模擬，研究了進氣流量和銳孔孔徑對氣泡形成的影響。VOF 法可以避免Level set法中的物理量不守恒問題，Level set法可以彌補VOF法中求解交界面法向和曲率的精度問題，因此，CLSVOF法結(jié)合了二者的優(yōu)點，為不可壓兩相流界面的模擬計算注入了新的活力。

1 數(shù)值解法

將水氣組成的兩相流體認(rèn)為是一種變密度的單流體，二者共享一套連續(xù)方程和動量方程。引入體積分?jǐn)?shù)F，F(xiàn)=0表示單元內(nèi)為純氣相，0＜F＜1表示單元內(nèi)為兩相混合，F(xiàn)=1表示單元內(nèi)為純水。對整個計算區(qū)域進行求解的控制方程如下：

N-S方程：

上面兩個方程中，v是速度矢量，P為壓強，F(xiàn)sv為每個單元的表面張力，ρ和μ分別為流體的密度和粘度，分別定義如下：

式中，σ是表面張力系數(shù)，κ是平均表面曲率，＾n是表面單元法向量，δs是表面delta函數(shù)。ρl和ρb是水和氣的密度，μl和μb是水和氣的粘度。

通過VOF法計算出體積分?jǐn)?shù)F平流方程可以得出運動界面：

在通過VOF法得出重構(gòu)后的交界面之后，在計算區(qū)域中，定義一個標(biāo)量函數(shù)φ，使得任意時刻水氣交界面剛好是φ的零等值面。新的交界面就可以通過求解Level set方程得到：

通過Level set方程，得到表面單元法向量：

平均表面曲率κ：

這樣，相交界面的法向向量和曲率則不再用不連續(xù)的F函數(shù)計算，而是通過連續(xù)的φ函數(shù)計算。上述方程采用有限體積法在計算網(wǎng)格上離散，速度、壓力、F函數(shù)和φ函數(shù)均定義在網(wǎng)格中心。

在CLSVOF法的基礎(chǔ)上根據(jù)不同的計算條件指定了不同的壁面粘附角進行了數(shù)值模擬工作，并對模擬結(jié)果采用自編的圖像分析算法進行處理，迅速準(zhǔn)確地獲取了氣泡的運動信息。

為了對氣泡形成特性有更好的了解，采用Fr數(shù)和Bo數(shù)來描述氣泡的形成特性。

公式（10）中U0是銳孔處的進氣速度，d0是銳孔直徑。Fr數(shù)和銳孔處的進氣速度平方成正比，Bo數(shù)和銳孔直徑的平方成正比。

2 計算結(jié)果分析與討論

2.1 模型驗證

實驗采用的水體密度和粘度分別是1000kg／m3和1×10-3kg／（m·s），氣體的密度和粘度分別是1kg／m3和1.8×10-3kg／（m·s），表面張力系數(shù)是0.0728N／m，模擬計算區(qū)域為4cm×4cm×12cm 的方形柱體。方形柱體底部設(shè)置不同直徑的銳孔（0.5，0.8和1mm）用于形成氣泡。模擬條件與實驗條件相同。四個側(cè)邊壁采用無滑移邊界條件，頂部為大氣壓出口，底部銳孔采用入流邊界條件。整個計算區(qū)域覆蓋結(jié)構(gòu)網(wǎng)格，網(wǎng)格間距0.05mm。這樣的精度足夠滿足計算需要。

氣泡運動的實驗在一個橫斷面面積為100mm×100mm，高度500mm 的垂直水槽中完成，實驗裝置見圖1。實驗系統(tǒng)主要由氣液裝置、照明設(shè)備、成像裝置、以及圖像處理部分組成，照明設(shè)備為日本NAC公司生產(chǎn)的MEMRECAM GX-1型高速攝像機，實驗中視窗范圍為4.5cm×3.6cm，拍攝分辨率為28.2pixel／mm，拍攝速度為1000幀／s，采用50W 鹵素?zé)粼诒趁婧蛡?cè)面補光。拍攝的圖像與計算機連接直接輸入到計算機上保存，保存的mcf文件通過與高速攝像機配套的GXLink 控制軟件可直接轉(zhuǎn)化為BMP格式的序列圖像文件。實驗氣源為壓縮空氣，經(jīng)穩(wěn)壓裝置后通過控制閥門開度調(diào)節(jié)進氣流量以產(chǎn)生單個氣泡。模型由透明的有機玻璃板制成，具有良好的透光性能。

圖1 實驗裝置示意圖Fig.1 Schematic diagram of the experimental apparatus

圖2（a）給出了數(shù)值模擬得出的氣泡膨脹拉伸過程和脫離過程，氣泡在形成初期，以球面的形狀鼓出銳孔平面，隨著氣體的注入，氣泡體積不斷增大，呈規(guī)則的球形，氣泡不斷擠壓周圍的水體，在氣泡尾部附近形成渦環(huán)，使氣泡尾部開始收縮，伴隨氣泡體積的進一步增大，氣泡尾部開始出現(xiàn)細長的頸的雛形，尾部渦環(huán)強度逐漸增大，氣泡變形加劇，氣泡在浮力的作用下整體開始向遠離銳孔的方向移動，只有一個很細的頸保持和銳孔的接觸，最后氣泡克服表面張力的作用，從細頸處整體脫離銳孔，以震蕩狀態(tài)開始上浮，同時銳孔處重復(fù)開始出現(xiàn)鼓出銳孔平面的氣泡，一個新的氣泡又開始生成。

模擬結(jié)果與模型實驗圖片非常吻合，在第四步中，氣泡尾部的氣流模擬得很清楚，模型實驗中的氣泡尾部氣流無法看得如此清晰，這在一定程度上彌補了模型實驗的不足。模擬結(jié)果和模型實驗的吻合表明：采用的CLSVOF 方法可以準(zhǔn)確地模擬出氣泡邊界。

圖2 不同時刻的氣泡形成形狀Fig.2 Bubble generated shape at different time

2.2 氣流量影響

圖3給出了模型實驗中在氣流量分別為2.36、3.53和4.71cm3／min，銳孔孔徑d0為0.5mm 情況下氣流量對氣泡形成和脫離以及上升運動的影響。模型實驗中氣泡變形與數(shù)值模擬結(jié)果中周圍流場分布情況見圖3。

從圖3和4中可以看出，模擬結(jié)果和實驗結(jié)果在氣泡形成周期、變形與上升過程都吻合良好，驗證了數(shù)值模擬方法的可靠性。

圖3 氣泡形成和運動過程（實驗結(jié)果）Fig.3 The process of bubble formation and movement（experimental result）

圖3（a）和4（a）中，氣流量Q為2.36cm3／min（Fr=8.16），將最初產(chǎn)生的氣泡稱為引導(dǎo)氣泡，引導(dǎo)氣泡脫離孔口后，經(jīng)過初始時段的形狀震蕩之后，以穩(wěn)定的橢球狀遠離孔口，氣泡尾流對新形成的氣泡影響不大；圖4（a）中氣泡形成周期為190ms，與圖3（a）模型實驗氣泡脫離時間186ms相比要長4ms，模擬得到的氣泡形成過程和上升軌跡與實驗基本吻合。

圖4 氣泡周圍流場分布情況（數(shù)值模擬結(jié)果）Fig.4 The flow field distribution around bubble（numerical result）

氣 流 量 從2.36 增 加 到3.53cm3／min（Fr從8.16增大到18.37），見圖4（b）。引導(dǎo)氣泡在脫離孔口后經(jīng)過形狀震蕩過程上升，氣泡開始對周圍流場產(chǎn)生較為強烈的擾動，并在尾部出現(xiàn)兩個不對稱分布的反向旋轉(zhuǎn)的渦環(huán)，這個不對稱渦環(huán)對孔口處新形成的氣泡開始產(chǎn)生影響，導(dǎo)致跟隨氣泡的脫離時間比引導(dǎo)氣泡短，并且脫離時的氣泡體積也比引導(dǎo)氣泡小，圖3（b）中氣泡在125ms的形成周期中，新產(chǎn)生的氣泡和第一個引導(dǎo)氣泡在形狀上有一定的差距，與圖4（b）吻合。

氣流量增大到4.71cm3／min后（Fr從18.37增大到32.65），實驗中氣泡形成周期縮短到96ms，模擬氣泡形成周期為100ms，比前面兩個氣流量條件下的形成周期分別少了90ms與30ms，見圖3（c）和4（c）。引導(dǎo)氣泡脫離銳孔孔口后，形狀很快發(fā)生改變，出現(xiàn)橢球狀和球帽狀交替的震蕩形狀，并且在上升過程中出現(xiàn)翻轉(zhuǎn)現(xiàn)象，氣泡變形對周圍流場的擾動更加明顯，跟隨氣泡的成長嚴(yán)重受到引導(dǎo)氣泡的尾流作用干擾，不斷左右搖擺，最終以不規(guī)則狀態(tài)脫離孔口。在整個過程中，引導(dǎo)氣泡變形劇烈，并在上升過程中伴隨搖擺與翻轉(zhuǎn)運動。

圖5為模型實驗和數(shù)值模擬氣泡的形成周期比較，圖中可以看到氣泡形成周期趨勢是隨著進氣流量Q的增大而縮短，氣流量越大，氣泡形成周期越短，引導(dǎo)氣泡的尾流對跟隨氣泡的成長影響越大，在上升過程中氣泡間的相互作用也越強烈。

圖5 實驗和模擬氣泡形成周期比較Fig.5 The comparison of measured and predicted bubble formation period

圖5中模擬氣泡形成周期均比實驗氣泡形成周期要長4～10ms不等，一部分原因是實驗條件下的水體雜質(zhì)和實驗溫度造成的誤差，另外在模擬中對孔口接觸角的處理現(xiàn)在還沒有找到很好的解決辦法，期望在后期進行進一步的探索，獲得更好的模擬結(jié)果。

2.3 孔徑影響

圖6（a）給出了d0分別為0.5、0.8和1mm，進氣量分別為2.36、6.03和9.42cm3／min情況下（即進口氣速U0=0.2m／s，F(xiàn)r=8.16、5.10、4.08；Bo=0.034、0.087、0.136）引導(dǎo)氣泡和跟隨氣泡在上升過程中的位置變化圖，圖中曲線的斜率即為氣泡的瞬時上升速度。圖6（b）和（c）分別為進口氣速U0=0.3m／s，0.4m／s條件下各孔徑形成氣泡的位置變化圖。

在圖6（a）中，進氣量為2.36cm3／min，d0=0.5mm（Fr=8.16，Bo=0.034）形成的氣泡在上升運動過程中氣泡間的相互影響基本可以忽略不計，引導(dǎo)氣泡與其它跟隨氣泡一樣，上升狀態(tài)穩(wěn)定，氣泡間距離較大。d0增大到0.8和1.0mm 后，引導(dǎo)氣泡上升速度相對其它氣泡慢，跟隨氣泡上升速度在引導(dǎo)氣泡尾流影響下逐步增大，二者最終發(fā)生合并，這是因為引導(dǎo)氣泡脫離銳孔之后，氣泡底部產(chǎn)生的尾流對后面氣泡產(chǎn)生吸附作用，促使跟隨氣泡產(chǎn)生一個加速過程，上升速度比引導(dǎo)氣泡快。通過合并前氣泡的運動軌跡可以看到，兩個氣泡間的間距越小，跟隨氣泡受到引導(dǎo)氣泡的影響就越大，上升速度就越大。圖6（c）中由于進氣量與孔徑d0均較大，氣泡在發(fā)生合并后隨之發(fā)生破裂，在圖中可以看到破裂后氣泡的上升軌跡。

圖6 氣泡位移過程圖Fig.6 Process of bubble displacement

圖6表明在相同進氣速度條件下，銳孔孔徑越大，氣泡間的相互作用越強烈，越容易發(fā)生氣泡合并與破裂現(xiàn)象。

為了進一步分析氣泡之間產(chǎn)生不同作用的形成機理，了解氣泡從未發(fā)生合并到發(fā)生合并甚至在氣泡合并后再次破裂的過渡階段之間的關(guān)系，采用Fr數(shù)與Bo數(shù)構(gòu)建了氣泡運動形態(tài)過渡圖。

從圖7中可以看到，隨著Bo數(shù)的增大，氣泡從不發(fā)生合并到發(fā)生合并最后到氣泡合并后破裂的過渡階段中，F(xiàn)r數(shù)隨之減小。d0為0.5mm，銳孔處的進氣速度為0.4m／s（Bo數(shù)為0.034，F(xiàn)r數(shù)為32.65），此時氣泡之間相互擾動不大，氣泡在600ms才發(fā)生合并現(xiàn)象；d0為1mm，銳孔處的進氣速度為0.3m／s（Bo數(shù)為0.136，F(xiàn)r數(shù)為9.18），氣泡在350ms發(fā)生合并與破裂，表明進氣速度改變了氣泡之間運動形態(tài)的過渡時間，這個時間又隨著孔徑的增大而不斷減小。為了驗證得到的結(jié)論，將模擬數(shù)據(jù)和實驗數(shù)據(jù)與Tufail［16］的實驗數(shù)據(jù)做了比較，Tufail的結(jié)論中，d0為0.78mm，銳孔處的進氣速度為5.25m／s（Bo數(shù)為0.083，F(xiàn)r數(shù)為3602.09），此時氣泡未發(fā)生合并；Zhang［6］的實驗結(jié)果中，d0為2mm，銳孔處的進氣速度為2.43 m／s（Bo數(shù)為0.545，F(xiàn)r數(shù)為300.93），此時氣泡發(fā)生合并，二者運動形態(tài)從未合并到合并的過渡過程中，進氣速度減小的同時伴隨著孔徑的增大，這與模擬得出的氣泡從不發(fā)生相互擾動到氣泡之間發(fā)生合并與破裂的過渡階段中，氣泡銳孔處的進氣速度隨著孔徑的增大而減小的結(jié)論一致。

圖7 氣泡Fr數(shù)和Bo 數(shù)關(guān)系圖Fig.7 The diagram of Fr and Bo

3 結(jié) 論

采用CLSVOF法結(jié)合表面張力模型同時配合實驗驗證對銳孔氣泡的形成過程及上升運動進行了研究，模擬了不同進氣流量Q和銳孔孔徑d0對氣泡形成的影響，數(shù)值模擬結(jié)果與模型實驗一致。表明所采用的CLSVOF 數(shù)值模擬方法是可行的，驗證了數(shù)值模擬方法的可靠性。

模 擬 了0.5mm ≤d0≤1.0mm，2.36 ≤Q≤18.84cm3／min條件下氣泡的形成過程，得出在同孔徑不同進氣量和同進氣速度不同孔徑時，大進氣量和大孔徑條件下形成的引導(dǎo)氣泡尾流作用對跟隨氣泡的形成影響較大，氣泡合并現(xiàn)象也出現(xiàn)得越早。氣泡的形成周期隨著進氣流量Q的增大而縮短，并且隨著進氣流量Q不斷增大，氣泡形成周期縮短的趨勢慢慢變得平緩。

通過Fr數(shù)和Bo數(shù)關(guān)系圖得出：氣泡從不發(fā)生相互擾動到氣泡之間發(fā)生合并與破裂的過渡階段中，氣泡銳孔處的進氣速度隨著孔徑的增大而減小。

［1］ HIRT C，NICHOLS B.Volume of fluid（VOF）method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，3：201-225.

［2］ CHEN C，F(xiàn)AN L.Discrete simulation of gas-liquid bubble columns and gas liquid solid fluidized beds［J］.A.I.Ch.E.Journal，2004，50：288-296.

［3］ SANKARANARAYANAN K，SHAN X，KEVREKIDIS I.Analysis of drag and virtual mass forces in bubbly suspensions using an implicit formulation of the lattice Boltzmann method［J］.Journal of Fluid Mechanics，2002，452：61-96.

［4］ GERLACH D，BISWAS G，DURST F.Quasi-static bubble formation on submerged orifices［J］.International Journal of Heat and Mass Transfer，2005，48：425-438.

［5］ IGUERA F.Injection and coalescence of bubbles in a very viscous liquid［J］.Journal of Fluid Mechanics，2005，530：369-378.

［6］ ZHANG L，SHOJI M.A periodic bubble formation from a submerged orifice［J］.Chemical Engineering Science，2001，5（6）：5371-5381.

［7］ MRTIN M，MONTES F，GALAN M.A Numerical calculation of shapes and detachment times of bubbles generated from a sieve plate［J］.Chemical Engineering Science，2006，6（1）：363-369.

［8］ DAVIDSON F，SCHULER G.Bubble formation at an orifice in a inviscous liquid［J］.Transactions of the Institute of Chemical Engineers，1960，38：335-342.

［9］ JAMIALAHMADI M，ZEHTABAN R，SARRAFI A.Study of bubble formation under constant flow conditions［J］.Chemical Engineering Research and Design，2001，79（A5）：523-532.

［10］KUMAR R，KULOOR R.The formation of bubbles and drops［J］.Advances in Chemical Engineering，1970（8）：256-368.

［11］梁剛濤，沈勝強，楊勇.單液滴撞擊平面液膜飛濺過程的CLSVOF模擬［J］.熱科學(xué)與技術(shù)，2012，11（1）：8-13.

［12］LIU Yong-tao，MA Ning，ZHU Ren-qing.CLSVOF method for violent sloshing with impact load in tanks under shallow water depth［J］.船 舶力學(xué)，2012，16（9）：990-1000.

［13］廖斌，陳善群.基于CLSVOF 方法的三維單個上升氣泡運動的數(shù)值模擬［J］.水動力學(xué)研究與進展，2012，27（3）：275-283.

［14］宋云超，王春海，寧智.追蹤不可壓縮兩相流相界面的CLSVOF方法［J］.農(nóng)業(yè)機械學(xué)報，2011，42（7）：26-33.

［15］BRACKBILL J，KOTHER D，ZEMACH C.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100（2）：335-354.

［16］TUFAILE A，SARTORELLI J.Chaotic behavior in bubble formation dynamics［J］.PHY A，2000，275：336-346.