翟建國(guó)，魯天龍，劉慶茶，黃筱云，胡勇虎

(1.中交武漢港灣工程設(shè)計(jì)研究院有限公司，武漢 430040；2.長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，長(zhǎng)沙 410114)

氣泡帷幕被廣泛應(yīng)用于各種實(shí)際工程中，如異齒裂腹魚(yú)的阻攔[1]、河面破冰工程[2]、異重流防護(hù)[3]、水下噪音隔離[4]、水下爆破防護(hù)[5]、氣幕式防波堤[6]、海上浮油攔截[7]等。在這些實(shí)際工程應(yīng)用中，主要以現(xiàn)場(chǎng)試驗(yàn)用于確定氣泡帷幕所需的孔徑、空氣流量以及孔口間距；但由于氣泡本身及其運(yùn)動(dòng)引起的流場(chǎng)的復(fù)雜性，目前對(duì)于氣泡帷幕的研究尚不完善，工程應(yīng)用時(shí)對(duì)于孔徑、空氣流量以及孔口間距的確定變得相對(duì)困難；為減少實(shí)際應(yīng)用中確定孔徑、空氣流量以及孔口間距的工作量，對(duì)于氣泡帷幕的研究是十分必要的。目前國(guó)內(nèi)外學(xué)者就氣泡帷幕的研究相對(duì)較少，但對(duì)單個(gè)及多個(gè)連續(xù)氣泡運(yùn)動(dòng)及其周?chē)鲌?chǎng)運(yùn)動(dòng)特征進(jìn)行了大量理論、實(shí)驗(yàn)及數(shù)值研究。如不可壓縮理想流體的氣泡輸運(yùn)方程在1917年就被Rayleigh[8]提出，為氣泡研究奠定了基礎(chǔ)理論；在實(shí)驗(yàn)方面主要采用粒子圖像測(cè)速儀(PIV)對(duì)氣泡及其運(yùn)動(dòng)流場(chǎng)進(jìn)行測(cè)量。如Liu等[9]采用粒子圖像測(cè)速儀(PIV)對(duì)矩形氣泡柱中氣泡鏈誘導(dǎo)的流動(dòng)結(jié)構(gòu)進(jìn)行了實(shí)驗(yàn)研究，得出隨著液體粘度的降低，氣泡上升軌跡由一維向三維變化，并在實(shí)驗(yàn)中觀察到自由渦旋，錯(cuò)流和不規(guī)則的循環(huán)流；Fan等[10]采用粒子圖像測(cè)速儀(PIV)對(duì)不同質(zhì)量濃度的非牛頓流體中兩個(gè)平行運(yùn)動(dòng)氣泡周?chē)牧鲌?chǎng)進(jìn)行了實(shí)驗(yàn)測(cè)量，并對(duì)速度矢量、速度等值線和速度分量進(jìn)行了分析。為了對(duì)氣泡及其周?chē)鲌?chǎng)特性進(jìn)行更好地分析，學(xué)者們對(duì)數(shù)值模擬也進(jìn)行了廣泛地應(yīng)用；如徐玲君等[11]對(duì)直徑為2.6 mm和4 mm的單個(gè)氣泡進(jìn)行了數(shù)值模擬研究；Bhusare等[12]采用OpenFoam中的歐拉-歐拉雙流體模型進(jìn)行了氣液數(shù)值模擬，對(duì)液體中氣泡的速度分布、含氣率分布和環(huán)流分布情況進(jìn)行了分析；吳恒等[13]應(yīng)用VOF方法，分析了不同狀態(tài)下的氣泡脫離時(shí)間及直徑；Saleh[14]通過(guò)使用k-ε閉合湍流模型對(duì)三維氣液流動(dòng)進(jìn)行了模擬，得到小氣泡由于沒(méi)有足夠的動(dòng)量離開(kāi)壁面，將在壁面跟隨液體循環(huán)運(yùn)動(dòng)。

根據(jù)目前的研究，國(guó)內(nèi)外學(xué)者對(duì)氣泡及其周?chē)鲌?chǎng)雖然已進(jìn)行了大量的實(shí)驗(yàn)和數(shù)值研究，但并未對(duì)氣泡運(yùn)動(dòng)引起的周?chē)鲌?chǎng)在斷面的分布進(jìn)行深入地分析。本研究采用OpenFoam中的VOF和標(biāo)準(zhǔn)k-ε模型對(duì)單孔出流氣泡周?chē)鲌?chǎng)特性進(jìn)行了三維模擬，在此基礎(chǔ)上通過(guò)改變孔口氣體流速、孔口直徑進(jìn)一步對(duì)單孔出流氣泡周?chē)鲌?chǎng)在斷面的分布進(jìn)行深入地分析，為氣泡帷幕周?chē)鲌?chǎng)運(yùn)動(dòng)特性研究打下基礎(chǔ)。

1 數(shù)值模擬

1.1 控制方程

本研究中，將氣體與水體均視為不可壓縮相，不可壓縮牛頓流體的連續(xù)性方程和動(dòng)量方程如下

(1)

(2)

相間邊界追蹤采用VOF模型，對(duì)于第q相，其體積分?jǐn)?shù)方程為

(3)

(4)

式中：α表示體積分?jǐn)?shù)；q表示為第q相。

計(jì)算單元中混合流體的密度和粘度的計(jì)算方程如下

(5)

(6)

式中：μ表示粘度，Pa·s；g表示氣體相。

在連續(xù)表面張力(CSF)模型中，氣體和水體兩相的體積力由式(7)給出

(7)

k-epsilon(k-ε)湍流模型是計(jì)算流體動(dòng)力學(xué)中最常用的模型，用于模擬湍流條件下的平均流動(dòng)特性。它是一個(gè)雙方程模型，通過(guò)兩個(gè)傳輸方程(PDE)給出湍流的一般描述。對(duì)于湍流動(dòng)能k和耗散率ε由下式定義

(8)

(9)

1.2 模型驗(yàn)證

表1 計(jì)算流體物理參數(shù)

1.2.1 氣泡形態(tài)驗(yàn)證

根據(jù)文獻(xiàn)[15]物理實(shí)驗(yàn)，選用其孔口氣體流速Vg=3.0 m/s，區(qū)域尺寸為50 mm×50 mm×100 mm，孔口直徑為1 mm的實(shí)驗(yàn)條件進(jìn)行驗(yàn)證。數(shù)值計(jì)算幾何模型及網(wǎng)格繪制如圖1所示；計(jì)算流體物理參數(shù)如表1所示；邊界條件設(shè)置底部孔口為氣體相速度入口邊界。數(shù)值模擬結(jié)果與文獻(xiàn)[15]實(shí)驗(yàn)結(jié)果對(duì)比如圖2所示，同時(shí)給出數(shù)值模擬與文獻(xiàn)物理實(shí)驗(yàn)的初始?xì)馀菝撾x時(shí)間與脫離直徑和上升速度對(duì)比如表2所示；根據(jù)對(duì)比結(jié)果，可以看出數(shù)值模擬得到的氣泡運(yùn)動(dòng)特性與文獻(xiàn)試驗(yàn)結(jié)果吻合較好。

圖1 幾何模型及網(wǎng)格示意圖

表2 數(shù)值模擬結(jié)果與文獻(xiàn)[15]實(shí)驗(yàn)結(jié)果對(duì)比表

1.2.2 氣泡周?chē)魉衮?yàn)證

選用文獻(xiàn)[9]實(shí)驗(yàn)中編號(hào)S-1下20 ml/min時(shí)的實(shí)驗(yàn)工況對(duì)氣泡周?chē)魉龠M(jìn)行驗(yàn)證；根據(jù)具體的實(shí)驗(yàn)工況，數(shù)值模擬驗(yàn)證采用底部開(kāi)孔，孔口直徑d為2 mm。通過(guò)提取模擬計(jì)算得到的點(diǎn)(0，0，33 mm)在x方向速度隨時(shí)間的變化情況，并與文獻(xiàn)此點(diǎn)實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果如圖3所示，由圖可見(jiàn)數(shù)值計(jì)算模擬結(jié)果與文獻(xiàn)物理實(shí)驗(yàn)結(jié)果數(shù)據(jù)趨勢(shì)基本一致。

1.3 模擬計(jì)算工況

本文計(jì)算區(qū)域幾何尺度為1 100 mm×150 mm×250 mm(長(zhǎng)×寬×高)，具體幾何模型如圖4所示，建立直角坐標(biāo)系，原點(diǎn)(0，0，0)位于底部開(kāi)孔中心。

圖3 數(shù)值結(jié)果與實(shí)測(cè)結(jié)果比較

表3 模擬試驗(yàn)參數(shù)值

計(jì)算模擬水深h=200 mm，孔徑及孔口氣體流速設(shè)計(jì)如表3所示。

2 計(jì)算結(jié)果與分析

2.1 速度等值線圖及流線分布

為了清楚地描述不同孔口直徑、不同孔口氣體流速對(duì)周?chē)w流場(chǎng)的影響，繪制了流線圖及不同孔徑及孔口氣體流速下氣泡上升運(yùn)動(dòng)過(guò)程中的速度等值線分布圖。

圖5表示孔口氣體流速Vg=1.0 m/s，孔口直徑d=1.5 mm時(shí)x=0斷面流線分布情況。由圖可見(jiàn)氣泡上升運(yùn)動(dòng)過(guò)程中，氣泡首先帶動(dòng)其周?chē)乃w向上運(yùn)動(dòng)，在水體到達(dá)水面后形成環(huán)流擴(kuò)散開(kāi)來(lái)。

圖6表示不同孔口氣體流速、孔口直徑下x=0斷面合速度等值線圖，其中圖6-a～6-e的孔口直徑d均為1.5 mm；圖6-f～6-j的孔口氣體流速Vg均為1.0 m/s。由圖可見(jiàn)氣泡上升運(yùn)動(dòng)帶動(dòng)的周?chē)鲌?chǎng)運(yùn)動(dòng)速度主要集中在垂直孔口上方附近兩側(cè)；且隨著孔口氣體流速、孔口直徑的增大，連續(xù)上升氣泡帶動(dòng)的周?chē)w運(yùn)動(dòng)速度相對(duì)增大，帶動(dòng)的周?chē)w運(yùn)動(dòng)范圍也增大。

2.2 斷面垂線平均合速度

圖5 x=0斷面流線分布圖

為了得到x=0斷面垂線平均合速度的分布情況，截取x=0斷面，參考方崇等[16]對(duì)垂向分布流速的研究，采用三點(diǎn)法對(duì)相對(duì)水深為0.2、0.6及0.8時(shí)的速度進(jìn)行提取，并進(jìn)行相應(yīng)計(jì)算。

水深h=0.2 m，孔口氣體流速Vg=1.25 m/s時(shí)，不同孔口直徑下相對(duì)水深分別為0.2、0.6及0.8時(shí)，合速度在y軸上的分布情況如圖7表示；其中圖7-a和圖7-b分別表示孔口直徑d=1.0 mm和d=1.5 mm；圖8表示水深h=0.2 m，孔口直徑d=1.0 mm，不同孔口氣體流速下相對(duì)水深分別為0.2、0.6及0.8時(shí)，合速度在y軸上的分布情況；其中圖8-a、圖8-b分別表示孔口氣體流速Vg=0.75 m/s和Vg=1.25 m/s；由圖7和圖8可見(jiàn)，水深h=0.2 m，不同孔口直徑、孔口氣體流速下相對(duì)水深為0.2、0.6及0.8時(shí)，合速度在y軸上的分布規(guī)律基本相同；同組工況、不同相對(duì)水深時(shí)，合速度在y軸上的分布差別較小，大體上均呈對(duì)稱分布，因此可以應(yīng)用三點(diǎn)法求斷面垂線平均合速度。

6-a Vg=0.5 m/s6-b Vg=0.75 m/s6-c Vg=1.0 m/s6-d Vg=1.25 m/s6-e Vg=1.5 m/s

7-a d=1.0 mm7-b d=1.5 mm8-a Vg=0.75 m/s8-b Vg=1.25 m/s

三點(diǎn)法求斷面垂線平均合速度的計(jì)算公式如式(10)所示，根據(jù)提取的相對(duì)水深為0.2、0.6及0.8時(shí)的速度，可以通過(guò)三點(diǎn)法進(jìn)行計(jì)算得到x=0斷面垂線平均合速度在y軸上的分布情況。

(10)

式中：Vmag(y)表示y點(diǎn)對(duì)應(yīng)的垂線平均合速度，m/s；Vh(y)表示y點(diǎn)對(duì)應(yīng)的相對(duì)水深為0.2、0.6及0.8時(shí)的合速度，m/s。

9-a d=1.5 mm 9-b Vg=1.0 m/s

圖9-a表示孔口直徑d=1.5 mm時(shí)，不同孔口氣體流速時(shí)x=0斷面垂線平均合速度在y軸上的分布；圖9-b表示孔口氣體流速Vg=1.0 m/s時(shí)，不同孔口直徑時(shí)x=0斷面垂線平均合速度在y軸上的分布；由圖可見(jiàn)，斷面垂線平均合速度在y軸上基本呈對(duì)稱分布；當(dāng)孔口氣體流速、孔口直徑較大時(shí)，斷面垂線平均合速度也相對(duì)較大。

2.3 相對(duì)速度預(yù)測(cè)

相對(duì)速度D主要受表面張力、粘性力、慣性力的作用，因此可用We、Re等參數(shù)對(duì)氣泡運(yùn)動(dòng)引起的周?chē)鲌?chǎng)特性進(jìn)行描述，D、xo、We、Re數(shù)的計(jì)算公式為

(11)

(12)

(13)

(14)

式中：Vmag表示垂線平均合速度，m/s；Vg為孔口出氣速度，m/s；xo表示相對(duì)孔徑；y表示y軸坐標(biāo)值，m；d為孔口直徑，m；ρl為液體的密度，kg/m3；ρg為氣體的密度，kg/m3；σ為表面張力系數(shù)，N/m；μl為液體黏度，kg/(s2·m)。

圖10 D與xo的關(guān)系圖

根據(jù)圖9-b中孔口直徑d=2.0 mm，孔口氣體速度Vg=1.0 m/s時(shí)x=0斷面垂線平均合速度的數(shù)據(jù)，應(yīng)用式(11)和式(12)進(jìn)行計(jì)算處理得到D與xo，其關(guān)系曲線如圖10所示；結(jié)合分布曲線對(duì)D與xo的關(guān)系曲線進(jìn)行直接高斯分布概率密度函數(shù)擬合，其相對(duì)速度分布預(yù)測(cè)公式如式(15)所示。

(15)

式中：yo、A、σ和μ為待定系數(shù)；由式(15)可知，在給定孔口氣體速度Vg和孔口直徑d的條件下，為了得到相對(duì)速度D的分布情況，需要對(duì)上式中含有的yo、A、σ和μ四個(gè)未知參數(shù)進(jìn)行預(yù)測(cè)，考慮到對(duì)稱性及水深h的影響，取μ=0；因此，采用h/d、We和Re數(shù)對(duì)上述未知參數(shù)進(jìn)行擬合。

11-a yo～Re11-b yo～We11-c yo～H/d

圖11表示yo與Re、We和H/d的關(guān)系圖，其中圖11-a表示yo與Re數(shù)的關(guān)系圖，圖11-b表示yo與We數(shù)的關(guān)系圖。由圖11-a、圖11-b可見(jiàn)yo與Re和We數(shù)的關(guān)系不明顯，故不使用Re和We數(shù)對(duì)yo進(jìn)行預(yù)測(cè)；由圖11-c可以看出yo總體隨H/d增大而減小。

從圖11可以看出：yo與H/d幾乎呈線性相關(guān)，因此可構(gòu)建σ與H/d之間的關(guān)聯(lián)式，通過(guò)對(duì)計(jì)算數(shù)據(jù)進(jìn)行擬合，找出yo的最佳預(yù)測(cè)模型

(16)

式中：ao、bo為待定系數(shù)；通過(guò)一元非線性方法對(duì)式(16)進(jìn)行計(jì)算數(shù)據(jù)擬合，得到其待定系數(shù)ao=0.011 9、bo=-5.09×10-5以及均方誤差S=0.002 827、平方絕對(duì)誤差Sr=9.138 4×10-7和相關(guān)系數(shù)r=0.891 1。

A與Re、We和H/d的關(guān)系如圖12所示。由圖12-a、圖12-b可見(jiàn)A整體隨Re數(shù)及We數(shù)的增大而減少；由圖12-c可以看出A與H/d關(guān)系不明顯，因此不使用H/d對(duì)A進(jìn)行預(yù)測(cè)。

12-a A～Re12-b A～We12-c A～H/d

從圖12-a、圖12-b可以看出：A隨著We數(shù)和Re數(shù)的增大呈現(xiàn)減小的趨勢(shì)，因此可構(gòu)建A與We數(shù)和Re數(shù)之間的冪函數(shù)關(guān)聯(lián)式對(duì)數(shù)據(jù)進(jìn)行擬合

(17)

式中：a1、b1、c1、f1為待定系數(shù)；Qi分別取We數(shù)和Re數(shù)，通過(guò)對(duì)上式應(yīng)用一元非線性方法對(duì)模擬數(shù)據(jù)進(jìn)行擬合，計(jì)算得到其待定系數(shù)a1、b1、c1、f1以及相關(guān)參數(shù)如表4所示。

表4 相關(guān)參數(shù)

圖13表示σ與Re數(shù)、We數(shù)及H/d的關(guān)系圖，由圖13-a、圖13-b可見(jiàn)σ與Re和We數(shù)的關(guān)系不明顯，圖中數(shù)據(jù)點(diǎn)較分散，因此不采用Re數(shù)和We數(shù)對(duì)σ進(jìn)行預(yù)測(cè)；由圖13-c可以看出σ隨H/d的增大而增大。

13-a σ～Re13-b σ～We13-c σ～H/d

圖13表示σ與Re數(shù)、We數(shù)及H/d的關(guān)系圖，由圖13-a、圖13-b可見(jiàn)σ與Re和We數(shù)的關(guān)系不明顯，圖中數(shù)據(jù)點(diǎn)較分散，因此不采用Re數(shù)和We數(shù)對(duì)σ進(jìn)行預(yù)測(cè)；由圖13-c可以看出σ隨H/d的增大而增大。

(18)

式中：a2、b2為待定系數(shù)；通過(guò)一元非線性方法對(duì)上式進(jìn)行計(jì)算數(shù)據(jù)擬合，得到其待定系數(shù)a2=-5.14、b2=0.096 5以及均方誤差S=2.181、平方絕對(duì)誤差Sr=0.790 5和相關(guān)系數(shù)r=0.961 4。

綜合式(16)～式(18)及表4，可以得到y(tǒng)o、σ與A最佳預(yù)測(cè)模型，具體擬合公式如下

(19)

(20)

(21)

圖14 實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)曲線D與xo關(guān)系圖

為了驗(yàn)證以上公式的合理性，選用文獻(xiàn)[9]中實(shí)驗(yàn)編號(hào)為S-4下Q=20 mL/min，孔口直徑d=2 mm工況下的實(shí)驗(yàn)數(shù)據(jù)對(duì)上式進(jìn)行驗(yàn)證；首先根據(jù)實(shí)驗(yàn)條件通過(guò)式(19)～式(21)計(jì)算得到相對(duì)速度D與xo的預(yù)測(cè)曲線；再根據(jù)實(shí)驗(yàn)數(shù)據(jù)得到的直線上速度分布的數(shù)據(jù)，根據(jù)式(11)、(12)計(jì)算得到D與xo，并一同繪制于圖14；由圖14可以看出實(shí)測(cè)數(shù)據(jù)基本分布在預(yù)測(cè)曲線周?chē)?，所以此模型是可以用于斷面垂線平均合速度的預(yù)測(cè)。

3 結(jié)論

本文應(yīng)用OpenFoam軟件對(duì)單孔出流氣泡進(jìn)行了三維模擬，通過(guò)改變孔口氣體流速、孔口直徑對(duì)氣泡運(yùn)動(dòng)引起的周?chē)鲌?chǎng)特性進(jìn)行了分析研究，得到以下結(jié)論：(1)連續(xù)上升氣泡會(huì)帶動(dòng)周?chē)乃w在上升氣泡兩側(cè)形成環(huán)流；(2)隨著孔口氣體流速、孔口直徑的增加，氣泡運(yùn)動(dòng)引起的周?chē)鲌?chǎng)速度逐漸增大，影響范圍也逐漸增大；(3)斷面垂線平均合速度隨孔口氣體流速、孔口直徑的增大而增大；(4)擬合出基于高斯分布函數(shù)的相對(duì)速度預(yù)測(cè)公式，且與實(shí)驗(yàn)數(shù)據(jù)吻合較好。