單孔出流氣泡周圍流場特性研究

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

(1.中交武漢港灣工程設計研究院有限公司，武漢 430040；2.長沙理工大學 水利工程學院，長沙 410114)

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

根據目前的研究，國內外學者對氣泡及其周圍流場雖然已進行了大量的實驗和數值研究，但并未對氣泡運動引起的周圍流場在斷面的分布進行深入地分析。本研究采用OpenFoam中的VOF和標準k-ε模型對單孔出流氣泡周圍流場特性進行了三維模擬，在此基礎上通過改變孔口氣體流速、孔口直徑進一步對單孔出流氣泡周圍流場在斷面的分布進行深入地分析，為氣泡帷幕周圍流場運動特性研究打下基礎。

1 數值模擬

1.1 控制方程

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

(1)

(2)

相間邊界追蹤采用VOF模型，對于第q相，其體積分數方程為

(3)

(4)

式中：α表示體積分數；q表示為第q相。

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

(5)

(6)

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

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

(7)

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

(8)

(9)

1.2 模型驗證

表1 計算流體物理參數

1.2.1 氣泡形態驗證

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

圖1 幾何模型及網格示意圖

表2 數值模擬結果與文獻[15]實驗結果對比表

1.2.2 氣泡周圍流速驗證

選用文獻[9]實驗中編號S-1下20 ml/min時的實驗工況對氣泡周圍流速進行驗證；根據具體的實驗工況，數值模擬驗證采用底部開孔，孔口直徑d為2 mm。通過提取模擬計算得到的點(0，0，33 mm)在x方向速度隨時間的變化情況，并與文獻此點實驗結果進行對比，對比結果如圖3所示，由圖可見數值計算模擬結果與文獻物理實驗結果數據趨勢基本一致。

1.3 模擬計算工況

本文計算區域幾何尺度為1 100 mm×150 mm×250 mm(長×寬×高)，具體幾何模型如圖4所示，建立直角坐標系，原點(0，0，0)位于底部開孔中心。

圖3 數值結果與實測結果比較

表3 模擬試驗參數值

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

2 計算結果與分析

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

為了清楚地描述不同孔口直徑、不同孔口氣體流速對周圍水體流場的影響，繪制了流線圖及不同孔徑及孔口氣體流速下氣泡上升運動過程中的速度等值線分布圖。

圖5表示孔口氣體流速Vg=1.0 m/s，孔口直徑d=1.5 mm時x=0斷面流線分布情況。由圖可見氣泡上升運動過程中，氣泡首先帶動其周圍的水體向上運動，在水體到達水面后形成環流擴散開來。

圖6表示不同孔口氣體流速、孔口直徑下x=0斷面合速度等值線圖，其中圖6-a～6-e的孔口直徑d均為1.5 mm；圖6-f～6-j的孔口氣體流速Vg均為1.0 m/s。由圖可見氣泡上升運動帶動的周圍流場運動速度主要集中在垂直孔口上方附近兩側；且隨著孔口氣體流速、孔口直徑的增大，連續上升氣泡帶動的周圍水體運動速度相對增大，帶動的周圍水體運動范圍也增大。

2.2 斷面垂線平均合速度

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

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

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

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

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

(10)

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

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

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

2.3 相對速度預測

相對速度D主要受表面張力、粘性力、慣性力的作用，因此可用We、Re等參數對氣泡運動引起的周圍流場特性進行描述，D、xo、We、Re數的計算公式為

(11)

(12)

(13)

(14)

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

圖10 D與xo的關系圖

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

(15)

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

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

圖11表示yo與Re、We和H/d的關系圖，其中圖11-a表示yo與Re數的關系圖，圖11-b表示yo與We數的關系圖。由圖11-a、圖11-b可見yo與Re和We數的關系不明顯，故不使用Re和We數對yo進行預測；由圖11-c可以看出yo總體隨H/d增大而減小。

從圖11可以看出：yo與H/d幾乎呈線性相關，因此可構建σ與H/d之間的關聯式，通過對計算數據進行擬合，找出yo的最佳預測模型

(16)

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

A與Re、We和H/d的關系如圖12所示。由圖12-a、圖12-b可見A整體隨Re數及We數的增大而減少；由圖12-c可以看出A與H/d關系不明顯，因此不使用H/d對A進行預測。

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

從圖12-a、圖12-b可以看出：A隨著We數和Re數的增大呈現減小的趨勢，因此可構建A與We數和Re數之間的冪函數關聯式對數據進行擬合

(17)

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

表4 相關參數

圖13表示σ與Re數、We數及H/d的關系圖，由圖13-a、圖13-b可見σ與Re和We數的關系不明顯，圖中數據點較分散，因此不采用Re數和We數對σ進行預測；由圖13-c可以看出σ隨H/d的增大而增大。

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

圖13表示σ與Re數、We數及H/d的關系圖，由圖13-a、圖13-b可見σ與Re和We數的關系不明顯，圖中數據點較分散，因此不采用Re數和We數對σ進行預測；由圖13-c可以看出σ隨H/d的增大而增大。

(18)

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

綜合式(16)～式(18)及表4，可以得到yo、σ與A最佳預測模型，具體擬合公式如下

(19)

(20)

(21)

圖14 實測數據與預測曲線D與xo關系圖

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

3 結論

本文應用OpenFoam軟件對單孔出流氣泡進行了三維模擬，通過改變孔口氣體流速、孔口直徑對氣泡運動引起的周圍流場特性進行了分析研究，得到以下結論：(1)連續上升氣泡會帶動周圍的水體在上升氣泡兩側形成環流；(2)隨著孔口氣體流速、孔口直徑的增加，氣泡運動引起的周圍流場速度逐漸增大，影響范圍也逐漸增大；(3)斷面垂線平均合速度隨孔口氣體流速、孔口直徑的增大而增大；(4)擬合出基于高斯分布函數的相對速度預測公式，且與實驗數據吻合較好。