氣泡上升過(guò)程中尾流演變的VOF數(shù)值模擬

顧英杰，楊偉棟，劉志遠(yuǎn)，羅志國(guó)，2，鄒宗樹(shù)，2

（1 東北大學(xué)冶金學(xué)院，遼寧沈陽(yáng)110819； 2 東北大學(xué)冶金學(xué)院多金屬共生礦生態(tài)化冶金教育部重點(diǎn)實(shí)驗(yàn)室，遼寧沈陽(yáng)110819）

引 言

在石油、冶金、礦物加工等過(guò)程中都存在泡狀流動(dòng)[1-3]。在這些流動(dòng)過(guò)程中，氣泡的形狀、上升軌跡以及上升速度都在隨時(shí)發(fā)生變化。當(dāng)氣泡的形狀和上升速度發(fā)生變化后其尾流也會(huì)發(fā)生變化，而尾流的變化可能會(huì)使氣泡周?chē)牧鲌?chǎng)發(fā)生變化，進(jìn)而影響氣泡上升過(guò)程中的運(yùn)動(dòng)行為[4-9]。初始直徑不同的氣泡上升過(guò)程中的運(yùn)動(dòng)行為不同[10-13]，其上升過(guò)程中尾流的演變也可能不相同。尾流運(yùn)動(dòng)狀態(tài)的不同會(huì)導(dǎo)致流場(chǎng)的擾動(dòng)程度不同，不同的流場(chǎng)對(duì)工業(yè)生產(chǎn)起著不同的作用。因此，深入研究初始直徑不同的氣泡自由上升過(guò)程中其尾流的演變具有重要意義。

當(dāng)流體以不同的流速流過(guò)不同形狀的物體時(shí)，其尾流的形態(tài)以及脫落頻率都不相同[10,14-15]。為了探究氣泡尾流的特性，國(guó)內(nèi)外學(xué)者做了很多實(shí)驗(yàn)和數(shù)值模擬的研究。目前大多數(shù)研究主要集中在液體的黏度、液體的運(yùn)動(dòng)方式和氣泡界面受污染程度對(duì)氣泡尾流的影響[6,16-19]。實(shí)驗(yàn)研究方面，de Vries等[16]和Brücker[17]使用PIV 技術(shù)研究了氣泡自由上升過(guò)程中其尾流的變化，實(shí)驗(yàn)結(jié)果表明氣泡的形狀為球形或橢球形，實(shí)驗(yàn)結(jié)果還表明自由上升氣泡的尾流由做順時(shí)針運(yùn)動(dòng)和逆時(shí)針運(yùn)動(dòng)的兩個(gè)渦旋構(gòu)成。數(shù)值模擬研究方面，前人將氣泡的位置或形狀固定來(lái)研究其尾流的演化[11,18-21]。倪明玖[22]使用VOF 方法研究了氣泡尾流形成的機(jī)理，但其沒(méi)有研究氣泡尾流的脫落過(guò)程。Xu 等[23]使用VOF 方法對(duì)金屬溶液中自由上升的氬氣泡尾流的運(yùn)動(dòng)狀態(tài)進(jìn)行了模擬研究，結(jié)果表明氣泡的尾流是從其下表面開(kāi)始脫落，且尾流脫落的過(guò)程會(huì)對(duì)氣泡的形狀產(chǎn)生影響。龐明軍等[18-19]數(shù)值模擬研究了氣泡界面受污染程度對(duì)其尾流的影響，研究中他們將氣泡視為位置固定的球。Antepara 等[24]使用數(shù)值模擬方法研究了形狀不同的氣泡上升過(guò)程中其尾流的演變過(guò)程，但其沒(méi)有探究氣泡長(zhǎng)軸與水平線存在夾角情況下氣泡尾流的運(yùn)動(dòng)狀態(tài)，并且其模擬研究中幾何模型的尺寸僅為氣泡初始直徑的2～2.75 倍，這無(wú)法完全消除壁面效應(yīng)對(duì)氣泡尾流的影響。自由上升氣泡的尾流和氣泡的其他運(yùn)動(dòng)行為之間的關(guān)系也是學(xué)者研究的重點(diǎn)。Lee 等[25]和Tripathi 等[26]研究結(jié)果表明氣泡尾流的脫落會(huì)影響其上升的運(yùn)動(dòng)軌跡。Cano-Lozano 等[21,27]通過(guò)數(shù)值模擬方法對(duì)氣泡的上升軌跡和尾流進(jìn)行了研究，結(jié)果表明氣泡的不穩(wěn)定上升是由氣泡尾流的不對(duì)稱脫落造成的。Gaudlitz 等[9]使用直接數(shù)值模擬方法研究了當(dāng)量容積直徑為5.2 mm 的氣泡自由上升過(guò)程中其尾流的變化，研究結(jié)果表明氣泡尾流脫落時(shí)的頻率和氣泡上升軌跡發(fā)生振蕩時(shí)的頻率相同。

目前，大多數(shù)的數(shù)值模擬研究氣泡尾流時(shí)將其視為一個(gè)形狀和位置固定的球或橢球，而氣泡的形狀和速度在其上升過(guò)程中是不斷變化的，因此研究氣泡尾流的演變時(shí)將其形狀和位置固定是不夠完善的。本文采用VOF 的方法對(duì)靜止水中氣泡上升過(guò)程中尾流的演變進(jìn)行了數(shù)值模擬，并分析了不同尺寸氣泡的尾流演變過(guò)程。

1 數(shù)值模擬

氣泡在自由上升的過(guò)程中其形狀是不斷變化的[28]。在眾多的多相流模型中，VOF 方法不僅能夠準(zhǔn)確地模擬出氣泡上升過(guò)程中的形狀變化、不穩(wěn)定上升等行為[6,29]，還能夠模擬出氣泡的尾流[22]。因此本文利用VOF 模型在商業(yè)計(jì)算流體力學(xué)（CFD）軟件Fluent中對(duì)不同初始直徑的氣泡自由上升的過(guò)程進(jìn)行了數(shù)值模擬。假設(shè)模擬的氣液兩相流是等溫系統(tǒng)，氣體和液體均為黏度固定的不可壓縮牛頓流體，同時(shí)假設(shè)液體初始狀態(tài)為靜止。

1.1 控制方程

質(zhì)量守恒方程和動(dòng)量守恒方程為

式中，t 為時(shí)間，s；u 為速度，m·s-1；ρ 為密度，kg·m-3；p 為壓強(qiáng)，Pa；μ 為動(dòng)力黏性系數(shù)，kg·(m·s)-1；Fs為表面張力，N·m-3；g為重力加速度，m·s-2。

VOF 模型是流體體積函數(shù)法，對(duì)于單相或多相體積分?jǐn)?shù)的連續(xù)方程如下。

對(duì)第q相

式中，αg為氣相體積分?jǐn)?shù)。本文中氣體為空氣[ρ=1.225 kg·m-3，μ=1.7894×10-5kg·(m·s)-1]。

對(duì)主相

式中，αl液相體積分?jǐn)?shù)。本文中液相為水[ρ=998.2 kg·m-3，μ=1.003×10-3kg·(m·s)-1]。

出現(xiàn)混合流體的計(jì)算單元中密度與黏度的計(jì)算方法為

式中，μl為液相黏度；μg為氣相黏度；ρl為液相密度；ρg為氣相密度。

對(duì)于水氣交界面上的表面張力，應(yīng)用Brackbill等[30]提出的表面張力模型（CFS）。該模型中把表面張力作為一個(gè)體積力Fs，對(duì)于兩相系統(tǒng)其表達(dá)式為

式中，σ 為表面張力系數(shù)，N·m-1；k 為界面的曲率；αl為液相體積分?jǐn)?shù)；ρm為混合相的密度，利用式(5)計(jì)算。

1.2 幾何模型和邊界條件

氣泡的上升行為會(huì)受到壁面的影響[25]，當(dāng)壁面到氣泡的距離大于氣泡直徑的3 倍時(shí)，可以忽略壁面效應(yīng)的作用[31]。數(shù)值模擬所用的幾何模型如圖1所示，尺寸為80 mm×150 mm。網(wǎng)格為結(jié)構(gòu)化網(wǎng)格，其總數(shù)量共計(jì)約61萬(wàn)個(gè)。出口設(shè)置為壓力出口，空氣在出口處的體積分?jǐn)?shù)為1。壁面采用無(wú)滑移邊界。模擬中氣泡初始的直徑為2.4、2.7、3.0、3.4、3.7 mm，氣泡的初始位置距離幾何模型底部15 mm。

用基于壓力的方法分離和求解模型方程。用壓力隱式的PISO方法求解壓力-速度耦合。動(dòng)量離散方法采用二階迎風(fēng)格式，氣液界面自由表面形狀采用QUICK 方案。殘差設(shè)定收斂的標(biāo)準(zhǔn)為1×10-5，時(shí)間步長(zhǎng)為1×10-5。

圖1 計(jì)算域和網(wǎng)格劃分Fig.1 Schematic diagram of computational domain and mesh

1.3 模型驗(yàn)證

本文采用如圖2所示的實(shí)驗(yàn)設(shè)置來(lái)驗(yàn)證模型的準(zhǔn)確性。實(shí)驗(yàn)所用的超白玻璃水箱尺寸為150 mm×200 mm×200 mm。實(shí)驗(yàn)氣源為空氣，其經(jīng)氣閥調(diào)節(jié)氣流量以產(chǎn)生單個(gè)氣泡。液體為超純水，且實(shí)驗(yàn)與模擬中保持液面高度一致。 高速攝影機(jī)（PASTCAM SA-Z）的拍攝頻率為4000 fps，采用背投光進(jìn)行補(bǔ)光。實(shí)驗(yàn)驗(yàn)證拍攝的氣泡初始直徑為2.7 mm 和3.7 mm，且每個(gè)氣泡直徑下至少有三次重復(fù)拍攝的結(jié)果。下文將模擬所得的結(jié)果與實(shí)驗(yàn)進(jìn)行了對(duì)比，如圖3 和圖4 所示。從圖3 中可以看出，氣泡上升過(guò)程的實(shí)驗(yàn)和模擬結(jié)果中其形狀都是由球形轉(zhuǎn)變?yōu)闄E球形，且橢球形的氣泡上升一段時(shí)間后，氣泡的長(zhǎng)軸與水平方向都會(huì)產(chǎn)生一個(gè)夾角。圖4 中的速度是多次實(shí)驗(yàn)的平均值，且氣泡的速度是由Image-Pro Plus 軟件計(jì)算得到。從圖4 中可以看出，氣泡的上升速度隨氣泡上升高度的增加都呈現(xiàn)為增大、減小、增大后趨于穩(wěn)定，且模擬值與實(shí)驗(yàn)值誤差較小。這表明了數(shù)值模擬結(jié)果的有效性和準(zhǔn)確性。

圖2 實(shí)驗(yàn)設(shè)置Fig.2 Experimental setup

2 結(jié)果與分析

2.1 尾流的形成及其演變

由文獻(xiàn)[10-11]可知，當(dāng)流體流經(jīng)障礙物時(shí)，靠近障礙物的表面會(huì)形成邊界層，邊界層外流體的黏性可以忽略不計(jì)，將其視為理想流體，邊界層內(nèi)流體的黏性不可以忽略。當(dāng)流體的流速增加到一定數(shù)值時(shí)，會(huì)產(chǎn)生邊界層分離現(xiàn)象，形成尾流。流體Reynolds 數(shù)的變化導(dǎo)致尾流的運(yùn)動(dòng)狀態(tài)發(fā)生變化。

圖3 氣泡上升過(guò)程的實(shí)驗(yàn)和數(shù)值模擬結(jié)果Fig.3 Experimental and simulated result of bubble rise process

圖4 氣泡速度隨氣泡上升高度的變化Fig.4 Evolution of bubble velocity with rising height of bubble

氣泡自由上升過(guò)程也會(huì)形成尾流，但其與流體流經(jīng)障礙物所形成的尾流不完全相同[5]。障礙物的尾流主要是邊界層分離所形成，而氣泡的尾流還存在氣體內(nèi)部循環(huán)的影響[22]。此外，當(dāng)障礙物形成尾流時(shí)，一般障礙物的位置和形狀都是不變的，而氣泡的位置和形狀在其尾流形成過(guò)程中時(shí)刻都在發(fā)生變化。圖5給出了初始直徑為2.7 mm的氣泡自由上升過(guò)程中其尾流隨時(shí)間的演化。從圖5中可以看出，隨著上升時(shí)間的增加，氣泡的形狀由球形轉(zhuǎn)變?yōu)闄E球形；當(dāng)上升時(shí)間為0.16 s 時(shí)，氣泡的長(zhǎng)軸與水平方向產(chǎn)生角度，且隨上升時(shí)間的增加，這個(gè)傾斜角度的大小會(huì)發(fā)生變化。靜止液體中自由運(yùn)動(dòng)物體傾斜角度的變化往往與物體周?chē)鲌?chǎng)的變化有關(guān)。Ern 等[32-33]研究了軸對(duì)稱平板在水中的自由上升過(guò)程，結(jié)果表明平板沿之字形路徑上升時(shí)，平板的傾斜角度與平板尾流的非對(duì)稱脫落是同時(shí)出現(xiàn)的，且其傾斜角度的變化和尾流的非對(duì)稱脫落都與平板周?chē)鲌?chǎng)的變化有關(guān)，平板和氣泡在水中的上升過(guò)程都是自由運(yùn)動(dòng)物體存在周期性運(yùn)動(dòng)的典型例子，上升過(guò)程中的運(yùn)動(dòng)行為十分相似。

圖5 氣泡上升過(guò)程尾流隨時(shí)間的演變（氣泡初始直徑為2.70 mm）Fig.5 Bubble wake evolution with time during bubble rise（initial diameter of bubble is 2.70 mm）

從圖5 中還可以看到，氣泡尾流的表現(xiàn)形式與de Vries 等[16]和Brücker[17]實(shí)驗(yàn)的結(jié)果一致，尾流由順時(shí)針運(yùn)動(dòng)和逆時(shí)針運(yùn)動(dòng)的兩個(gè)渦旋構(gòu)成。當(dāng)氣泡的上升時(shí)間為0.06 s 時(shí)，氣泡的尾流從其下表面開(kāi)始形成，此時(shí)氣泡的形狀接近標(biāo)準(zhǔn)的橢球形；當(dāng)上升時(shí)間為0.10 s 時(shí)，尾流沿豎直方向的對(duì)稱性好，且尾流的左右兩個(gè)渦旋在同一條水平線上開(kāi)始脫落；隨著上升時(shí)間的增加，尾流的長(zhǎng)度和寬度增加（如0.14 s），且氣泡的形狀在表面自由能的作用下會(huì)趨于球形；當(dāng)上升時(shí)間為0.16 s 時(shí)，尾流左邊渦旋的寬度明顯大于右邊渦旋的寬度，尾流的對(duì)稱性變差，此時(shí)氣泡的形狀接近球帽形，并且隨氣泡上升時(shí)間的進(jìn)一步增加，氣泡的長(zhǎng)軸與水平方向的夾角變得更明顯（圖5 中0.18 s）；當(dāng)上升時(shí)間在0.1～0.22 s 之間，雖然尾流的對(duì)稱性隨上升時(shí)間的增加變差，且左右兩個(gè)渦旋也并不在同一條水平線上脫落，但尾流都是“成對(duì)”脫落的。隨著氣泡上升時(shí)間的繼續(xù)增加，尾流不再“成對(duì)”脫落（圖5 中0.22 s），且當(dāng)氣泡的上升過(guò)程穩(wěn)定以后，尾流會(huì)轉(zhuǎn)變?yōu)橹芷谛悦撀洌▓D5 中的0.42 s）。在氣泡上升過(guò)程中其形狀從球形轉(zhuǎn)變?yōu)闄E球形。氣泡尾流對(duì)稱脫落出現(xiàn)在氣泡形狀為球形的直線上升過(guò)程；當(dāng)氣泡的長(zhǎng)軸與水平方向出現(xiàn)夾角時(shí)，尾流進(jìn)入過(guò)渡態(tài)——非穩(wěn)定脫落狀態(tài)；隨著上升時(shí)間的進(jìn)一步增加，尾流最終表現(xiàn)為周期性脫落的運(yùn)動(dòng)狀態(tài)。對(duì)稱脫落是尾流運(yùn)動(dòng)狀態(tài)的初始態(tài)，尾流的對(duì)稱性好，且尾流的左右兩個(gè)渦旋在同一條水平線上脫落。過(guò)渡態(tài)是尾流運(yùn)動(dòng)狀態(tài)的非穩(wěn)定性脫落時(shí)期，尾流的左右兩個(gè)渦旋的對(duì)稱性隨時(shí)間的增加變差，且尾流的左右兩個(gè)渦旋并不在同一條水平線上脫落。周期性脫落是尾流的最終態(tài)，尾流的左右兩個(gè)旋渦周期性的交替脫落，且當(dāng)氣泡上升過(guò)程穩(wěn)定時(shí)，尾流周期性脫落的頻率穩(wěn)定。

圖6 不同初始直徑的氣泡上升過(guò)程中的尾流Fig.6 Bubble wake during bubble rising with different initial diameters

2.2 不同直徑氣泡的尾流演變

前人研究結(jié)果表明，隨著氣泡的上升，氣泡外部的渦流會(huì)形成尾部駐渦，尾部駐渦隨著氣泡的上升會(huì)形成尾流，且在氣泡上升過(guò)程中其尾流有對(duì)稱脫落和非對(duì)稱脫落[22,26]。圖6 是不同初始直徑的氣泡上升過(guò)程中其尾流的演變過(guò)程。從圖6中可以看出，氣泡的初始直徑在2.4 ～3.7 mm 時(shí)，氣泡上升過(guò)程中其尾流都是對(duì)稱脫落、過(guò)渡態(tài)和周期性脫落這三種運(yùn)動(dòng)狀態(tài)。當(dāng)尾流為對(duì)稱脫落時(shí)，隨著氣泡的初始直徑的增加，尾流的長(zhǎng)度和寬度增加，且尾流的左右兩個(gè)渦旋與氣泡的中心線之間的距離變大。當(dāng)尾流表現(xiàn)為過(guò)渡態(tài)時(shí)，隨著氣泡的初始直徑的增大，尾流渦旋的不對(duì)稱程度增大。

氣泡初始直徑的不同只改變尾流的長(zhǎng)度和寬度，不改變尾流的三種運(yùn)動(dòng)狀態(tài)。剛體尾流的運(yùn)動(dòng)狀態(tài)和渦旋尺寸主要與流體的Reynolds數(shù)和剛體的形狀有關(guān)[15]。在本研究中，氣泡周?chē)囊后w被假設(shè)為靜止的，所以流體的速度可以近似認(rèn)為是氣泡的速度；氣泡的形狀可以表示為繞流物體的形狀。因此本研究中流體的Reynolds 數(shù)可以使用氣泡的Reynolds 數(shù)來(lái)代替，所以本文所指的Reynolds 數(shù)均為氣泡的Reynolds數(shù)。Reynolds數(shù)的計(jì)算公式如下

式中，ρl為液體密度；d 為任意時(shí)刻氣泡當(dāng)量直徑，mm；ug為任意時(shí)刻氣泡速度；μ為動(dòng)力黏性系數(shù)。

從圖7中可以看出：隨著氣泡直徑的增大，氣泡尾流發(fā)生轉(zhuǎn)變的臨界Reynolds數(shù)增大。這主要是因?yàn)?，初始直徑?.4 ～3.7 mm 氣泡的形狀隨其直徑的增加變化較大[29]。氣泡形狀的變化也代表著氣泡當(dāng)量直徑的變化。此外，結(jié)合圖6，還可以發(fā)現(xiàn)，氣泡的初始直徑越大，則尾流發(fā)生轉(zhuǎn)變的臨界Reynolds數(shù)越大，且其尾流的長(zhǎng)度和寬度也越大。尾流的長(zhǎng)度和寬度對(duì)生產(chǎn)過(guò)程有重要意義。例如在利用氣泡的尾流增加尾礦絮凝沉降時(shí)，尾流的長(zhǎng)度和寬度越大，能影響的礦物范圍也越大，提高了生產(chǎn)效率[34]。

圖7 臨界Reynolds數(shù)隨氣泡直徑的變化Fig.7 Evolution of critical Reynolds number with different diameter of bubble

2.3 尾流周期性脫落的頻率

氣泡的尾流脫落會(huì)影響其周?chē)黧w的流場(chǎng)[34]。氣泡上升過(guò)程中其尾流脫落所引起的力占浮力的20%～30%[5]。氣泡尾流脫落時(shí)的頻率與氣泡上升軌跡發(fā)生振蕩時(shí)的頻率保持一致[9]。氣泡上升過(guò)程的軌跡、速度、變形、受力以及流場(chǎng)的擾動(dòng)等都存在周期性的變化，且這些過(guò)程都與氣泡尾流的周期性脫落密切相關(guān)。葉片尾流脫落時(shí)的頻率（f）可以使用式（9）來(lái)進(jìn)行計(jì)算[35]。所以，下文對(duì)氣泡尾流周期性脫落的頻率進(jìn)行分析，并將數(shù)值模擬和公式計(jì)算的結(jié)果進(jìn)行了對(duì)比，如圖8所示。

式中，u 是來(lái)流速度，模擬中假定氣泡周?chē)后w靜止，采用氣泡速度；d 是垂直于流體方向，流體形成脫流尾跡處的最大寬度，本文為氣泡的當(dāng)量直徑；Sr是Strouhal數(shù)，與Reynolds數(shù)有關(guān)，取值為0.3。

圖8是初始直徑不同的氣泡上升過(guò)程中尾流周期性脫落的頻率。從圖8 中可以發(fā)現(xiàn)，計(jì)算的結(jié)果與數(shù)值模擬的結(jié)果吻合較好。這表明，在數(shù)值模擬氣泡自由上升的過(guò)程中，當(dāng)氣泡的上升過(guò)程穩(wěn)定時(shí)，可以使用式（9）直接計(jì)算其尾流周期性脫落的頻率。從圖8 中還可以看到，隨著氣泡初始直徑的增大，氣泡尾流周期性脫落的頻率減小。這主要是因?yàn)闅馀菸擦鞯拿撀漕l率與氣泡的上升速度呈正比，而與氣泡的當(dāng)量直徑呈反比。直徑為2.4 ～3.7 mm 氣泡，隨著氣泡初始直徑的增加，當(dāng)氣泡上升過(guò)程穩(wěn)定之后，氣泡的上升速度變化較小[28]，但氣泡的初始直徑越大，氣泡上升過(guò)程中的變形程度越大，當(dāng)量直徑的變化也就越大[29]。

圖8 初始直徑不同的氣泡上升過(guò)程中尾流周期性脫落的頻率Fig.8 Frequency of periodic shedding of bubble wake during bubble rise with different initial diameters

3 結(jié) 論

采用VOF 方法模擬了初始直徑不同的氣泡在靜止水中的自由上升過(guò)程，系統(tǒng)地分析了氣泡上升過(guò)程中其尾流的演變。本文的研究結(jié)果可對(duì)氣泡上升過(guò)程中其運(yùn)動(dòng)行為的研究提供參考。主要結(jié)論如下。

（1）本研究范圍內(nèi)，氣泡上升過(guò)程中的尾流主要有對(duì)稱脫落、過(guò)渡態(tài)和周期性脫落三種運(yùn)動(dòng)狀態(tài)。對(duì)稱脫落是尾流的初始態(tài)，其對(duì)稱性好，且尾流的左右兩個(gè)渦旋在同一條水平線上脫落。過(guò)渡態(tài)是尾流的非穩(wěn)定性脫落時(shí)期，其對(duì)稱性隨氣泡的上升變差，且尾流的左右兩個(gè)渦旋并不在同一條水平線上脫落。周期性脫落是尾流的最終態(tài)，其尾流的左右兩個(gè)渦旋隨著氣泡的上升表現(xiàn)為周期性的交替脫落。

（2）對(duì)于2.4 ～3.7 mm 氣泡尾流的三種運(yùn)動(dòng)狀態(tài)不隨其初始直徑的變化而變化，且尾流運(yùn)動(dòng)狀態(tài)發(fā)生轉(zhuǎn)變的臨界Reynolds 數(shù)隨氣泡直徑的增大而增加。