靜水中單氣泡運動特性實驗研究

李文強，焦守華，唐 珂，楊宜昂，曲文海，柴 翔

(1.上海交通大學 核科學與工程學院，上海 200240；2.普渡大學 核工程系，美國 西拉法葉 47906)

氣液兩相流動和傳熱廣泛存在于核能工程、石油工程等工業領域以及自然科學研究領域。核反應堆正常工況運行時，定位格架前后的兩相流特性將大大影響反應堆燃料組件的熱工水力性能[1]，當反應堆發生事故，且堆芯熔融物融穿壓力容器下封頭掉進堆坑內后，由于堆芯余熱導致的堆坑水內兩相流動特性將顯著影響放射性物質的釋放速率[2]。板式核燃料組件矩形流道內的流動換熱研究中[3]，矩形通道內的熱工特性與氣泡的生成和分布密切相關，氣泡的生成將極大影響矩形通道內的熱工特性。在海上石油鉆井平臺中，氣態天然氣和水的混合物形成兩相流動，顯著增強了鉆井內流動阻力，增加了海底井口的堵塞風險[4]。對氣泡的運動速度、形態和軌跡特性進一步研究將有助于提高工程的安全等級或生產效率。針對單氣泡上升運動特性，文獻[5-11]進行了廣泛的理論和實驗研究，然而由于靜止水中氣泡運動過程存在著隨機性，氣泡上升軌跡對壁面與氣泡發生點距離、氣泡產生方法等初始條件敏感，氣泡運動軌跡、變形程度和瞬時速度之間存在著復雜的相互作用，因此迄今為止仍未形成完整、通用、準確的末速度預測關系式。

本文采用高速攝像機和圖像識別技術，對不同等效體積直徑下氣泡的運動軌跡、變形程度、瞬時速度和末速度進行研究。

1 實驗裝置

實驗裝置如圖1所示，主要由氣泡發生裝置、實驗段、光源系統、高速攝像機和圖像處理程序組成。氣泡發生裝置由注射泵、針閥、毛細管、橡膠塞組成，注射泵采用雷弗TSD-01雙向推拉型注射泵，注射速率控制范圍為1 μL/min～87 mL/min，通過同步器快速切換針閥達到快速切換注射器正負行程的目的。毛細管內徑分別為0.3、0.7、1、2、3、4和5 mm，采用可控緩慢生長的方式產生氣泡，橡膠塞用于固定可更換的毛細管，用于產生等效體積直徑為2.5～7.11 mm的氣泡。為消除壁面效應影響，實驗段采用200 mm×200 mm×600 mm的透明亞克力水箱組成，為保證方形水箱內具有很低的導熱系數來排除外界溫度變化的影響，水箱厚度達20 mm，方形水箱內為去離子水工質。實驗照明系統為200 mm×200 mm的平行白色光源，背面打光。采用400幀/s的幀率拍攝氣泡上升過程。本文實驗條件接近Liu等[12]的實驗條件，氣泡脫離毛細管150 mm后，速度已充分發展，達到穩定狀態，基于保守性的考慮，本文選取的拍攝中心位于管嘴出口265 mm處，測量面積為60 mm×90 mm的矩形區域。測量得到的運動圖像通過圖像處理程序進行處理，最終得到氣泡運動瞬時坐標、瞬時速度、變形程度等相關信息。

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

2 圖像處理技術和誤差分析

2.1 圖像處理技術

圖形處理流程如圖2所示，經過背景減除、對比度增強、邊緣識別、形態學處理等過程后即可獲得圖3所示的氣泡輪廓圖像，進一步處理即可獲得氣泡速度、氣泡變形程度和氣泡等效直徑等信息。計算氣泡速度前先通過下式計算氣泡質心：

(1)

(2)

式中：xi和zi為氣泡輪廓內的所有像素坐標；X和Z為求得的氣泡質心坐標；N為記錄的氣泡上升過程所用的總幀數。

圖2 圖形處理流程Fig.2 Flow chart of image processing

圖3 氣泡的長軸短軸示意圖Fig.3 Schematic diagram of bubble major axis and minor axis

氣泡速度為：

(3)

(4)

氣泡的變形程度和等效體積直徑均與氣泡圖像長軸a和短軸b有關(圖3)，氣泡的變形程度采用縱橫比E表示：

(5)

氣泡等效體積直徑de采用下式計算：

(6)

不同邊緣識別算法的縱橫比和等效體積直徑對比如圖4所示，經過Canny、Robert、Sobel、Prewitt邊緣算子得到的縱橫比和等效體積直徑差異很小，不同算子的圖像處理過程對后續得到氣泡運動參數敏感性很小，因此理論上4種算子均可采用。鑒于Canny算子通常具有更強的抑制噪聲的能力[13]，本文選擇Canny算子。

2.2 氣泡運動參數不確定度分析

借鑒Celata等[14]和Zhang等[15]氣泡運動參數不確定度處理方法，對氣泡末速度、縱橫比和等效體積直徑進行不確定度分析。圖像識別過程中，氣泡邊緣處某個方向上多識別或少識別1個像素所造成的質心不確定度為±0.5像素，氣泡長軸和短軸的識別不確定度為±1像素。

圖4 不同邊緣識別算法的縱橫比(a)和等效體積直徑(b)的對比Fig.4 Comparison of aspect ratio (a) and equivalent volume diameter (b) for different edge detection algorithms

根據不確定度傳遞公式，氣泡末速度的絕對不確定度εu為：

εu=

(7)

式中：u為氣泡末速度；Z1、Z2分別為氣泡前后兩幀的像素坐標。

氣泡等效體積直徑不確定度εde和縱橫比不確定度εE分別為：

(8)

(9)

根據以上分析方式得到了氣泡末速度、等效體積直徑和縱橫比的相對不確定度分別為0.67%、2.17%和4.63%，所有參數的相對不確定度均小于5%。

3 實驗結果和討論

3.1 氣泡運動軌跡和變形程度

圖5示出不同等效體積直徑氣泡運動的軌跡。由圖5可見，不同等效體積直徑的氣泡在上升過程中受到空間不對稱的湍流渦脫離的影響，運動軌跡呈現出不同振蕩幅度的運動形態，等效體積直徑小于5 mm的氣泡呈現出螺旋形運動，等效體積直徑大于5 mm的氣泡逐漸變為Z字形運動。

氣泡在浮力作用下開始逐漸加速到末速度狀態，受到阻力和浮力影響，氣泡上部與下部存在一較大的壓力差，壓力差作用下氣泡底部形成渦環并引導射流向上運動引發氣泡變形。氣泡直徑越大，Eo越大，氣泡頂部和底部誘導出的壓力差越大[16]，氣泡底部射流越大，而表面張力越弱，表面張力引起的氣泡內部附加壓強越小，因此氣泡變形越大。由圖5不難發現，氣泡直徑越大，變形程度越明顯。

3.2 氣泡振蕩幅度、瞬時速度和瞬時縱橫比

單氣泡上升過程中變形程度、瞬時速度、運動軌跡之間存在復雜的相互耦合關系，氣泡變形程度的瞬時變化程度會影響氣泡上升過程中曳力受力面積，受力面積的瞬時變化進而會影響氣泡在橫向和縱向運動的速度，氣泡橫向運動的瞬時變化速度會影響氣泡上升過程中氣泡橫向遷移的軌跡，氣泡縱向速度的瞬時變化會影響曳力大小，而曳力大小瞬時變化將改變氣泡運動過程中的受力平衡，最終影響氣泡變形程度。

圖6示出不同等效體積直徑下氣泡的橫向振蕩幅度，圖7示出不同等效體積直徑下氣泡瞬時速度與瞬時縱橫比的對比。由圖6、7可見，氣泡等效體積直徑為2.51 mm時，氣泡橫向速度振蕩幅度為150 mm/s，橫向振蕩幅度為6 mm，當氣泡直徑增加到6.86 mm時，氣泡橫向速度振蕩幅度為130 mm/s，橫向振蕩幅度為3 mm，因此氣泡橫向振蕩幅度與橫向速度vx峰值呈正比。

由圖7可見，氣泡等效體積直徑由2.51 mm逐漸增加至6.86 mm時，氣泡瞬時縱橫比振蕩幅度由0.075增大至0.2，氣泡的瞬時縱向速度vz的振蕩幅度由12.5 mm/s增大至50 mm/s，氣泡瞬時橫向速度振蕩幅度由150 mm/s降低至130 mm/s。這表明，小直徑氣泡的瞬時縱橫比振蕩幅度較小，氣泡瞬時橫向速度振蕩幅度較大，瞬時縱向速度振蕩幅度較小；大直徑氣泡瞬時縱橫比振蕩幅度較大，氣泡瞬時橫向速度振蕩幅度較小，而氣泡瞬時縱向振蕩頻率較大。其原因為氣泡直徑增大后，表面張力對氣泡變形程度的影響逐漸變弱，氣泡周圍流場的變化更容易導致氣泡發生大幅度的變形，氣泡瞬時縱橫比的振蕩幅度增大，氣泡的變形程度最終影響了氣泡橫向速度的變化。從不同直徑下氣泡瞬時橫向速度、瞬時縱向速度及瞬時縱橫比的變化頻率來看，氣泡縱向運動頻率和縱橫比變化頻率為氣泡橫向運動頻率的兩倍，不隨直徑的變化而發生變化。同樣不難發現，氣泡瞬時縱橫比的大小與瞬時縱向速度之間呈反比關系，這種反比關系可根據Moore[17]的理論進行解釋：氣泡速度越大，氣泡頂部的動壓更強，因此變形程度更大，瞬時縱橫比越小；反之氣泡速度越小，氣泡頂部的動壓更弱，因此變形程度更小，瞬時縱橫比越大。

圖6 不同等效體積直徑下氣泡的橫向振蕩幅度Fig.6 Bubble axial oscillation amplitude for different equivalent volume diameters

等效體積直徑：a——2.51 mm；b——3.75 mm；c——5.02 mm；d——6.86 mm圖7 不同等效直徑下氣泡瞬時速度與瞬時縱橫比的對比Fig.7 Comparison of bubble instaneous velocity and instaneous aspect ratio under different equivalent volume diameters

圖8示出氣泡橫向位移、橫向速度和縱橫比的對比，圖9示出氣泡運動坐標示意圖。由圖8、9可見：當氣泡速度從0逐漸加速到速度最大時，氣泡的縱橫比逐漸變小，氣泡逐漸接近橢球形，位置越接近中軸線。根據Fan等[18]的研究結果，氣泡在上升過程中下表面附近形成了不對稱的湍流渦，其周期性從氣泡表面脫落，影響了氣泡周圍流場的分布和氣泡受力，從而造成上升速度和氣泡橫縱比周期性變化。

圖8 氣泡橫向位移、橫向速度和縱橫比示意圖Fig.8 Schematic diagram of lateral motion,lateral velocity and aspect ratio of bubble

3.3 氣泡末速度與等效體積直徑的關系

在工程應用領域，更關注的是氣泡的末速度vT，在分析氣泡末速度時，多數學者僅考慮浮力和曳力：

(10)

式中：g為重力；ρL、ρg分別為水和氣體的密度；Cd為曳力系數。

圖9 氣泡運動坐標示意圖Fig.9 Schematic diagram of bubble motion coordinate

氣泡末速度由重力、氣液密度之比、等效體積直徑、曳力系數共同決定。其中：Rybczynski[5]從理論出發推導出了Re小于1的適用于球形氣泡的運動末速度關系式；Wallis[6]將該關系式推廣到Re介于1～100之間，適用于微變形的球形氣泡；考慮到氣泡運動過程中存在的巨大變形，Davies等[7]推導得到了一個半經驗末速度運動關系式；Mendelson[8]從波動理論出發，得到了氣泡直徑介于1.4～6 mm之間的末速度理論關系式；Tomiyama等[10]從勢流理論出發，推導得到了純凈水中氣泡末速度與氣泡縱橫比、氣泡直徑和物性的顯式關系式；Park等[11]總結了文獻中不同直徑下單氣泡運動末速度實驗數據，獲得了氣泡末速度與氣泡直徑、物性之間的通用顯式關系式。

圖10示出不同等效體積直徑下氣泡末速度的對比。由圖10可見，Rybczynski、Davies等的關系式在預測等效體積直徑為2～8 mm的氣泡末速度時，預測值與實驗值存在100%以上的相對誤差，不再適用。Mendelson的關系式在預測小直徑氣泡末速度關系式時存在過高估計，直徑越小，誤差越大。Wallis的關系式在預測等效體積直徑為2～5 mm的氣泡末速度時，存在20%左右的相對誤差。Tomiyama等的關系式低估了單氣泡運動在等效體積直徑為2～8 mm之間的末速度，這與Tang[19]的結論一致。Park等的顯式關系式與實驗值相比偏高，其主要原因為該關系式采用函數構造的方式獲得，通過抑制因子fsc來調節水中雜質對球形氣泡區域阻力系數的影響，而在純凈水中Park等認為fsc=1。

圖10 不同等效體積直徑下氣泡末速度的對比Fig.10 Comparison of bubble terminal velocity with different equivalent volume diameters

為獲得更加通用的精確顯式關系式，借鑒Park等的關系式重新構建一新的顯式關系式，引入A2和B2兩個參數修正關系式中的分母，并根據實驗數據擬合出如下關系式：

(11)

采用均方根誤差RMS來衡量各關系式間的預測精度，其定義如下：

(12)

式中：vei為實驗獲得的氣泡末速度；vsi為通過關系式擬合得到的氣泡末速度。經計算，新關系式、Park、Tomiyama和Wallis等的關系式的RMS分別為2.75%、3.63%、10.01%和20.12%，新關系式的RMS明顯更小。因此在純凈水中，新關系式在預測等效體積直徑介于2～8 mm之間的氣泡末速度時精度更高。

4 結論

針對單氣泡上升過程，本文從氣泡運動軌跡、氣泡瞬時速度、氣泡瞬時縱橫比、氣泡末速度等方面對氣泡運動特性進行了研究，研究結論如下。

1) 氣泡橫向振蕩幅度與氣泡橫向速度峰值呈正比，兩者均隨氣泡等效體積直徑的增加而降低。

2) 氣泡瞬時縱橫比變化頻率與氣泡瞬時縱向速度變化頻率一致，是氣泡瞬時橫向速度變化頻率的兩倍，與氣泡等效體積直徑無關。

3) 隨氣泡等效體積直徑增大，氣泡瞬時縱橫比振蕩幅度不斷加大，瞬時橫向速度振蕩幅度不斷減小，縱向速度振蕩幅度減小。

4) 大多數氣泡在橫向偏移為0時橫向速度趨向于最大值，縱橫比趨向于最小值，反之亦然。氣泡下表面渦的不對稱脫離是導致氣泡縱橫比、上升速度、橫向位置周期性變化的原因。

5) 氣泡縱橫比大小與縱向速度大小呈反比關系，即氣泡縱橫比最大時氣泡縱向速度最小，氣泡縱橫比最小時氣泡縱向速度最大。

6) 本文針對氣泡末速度關系式目前絕大多數為隱式關系式的弊端，提出了兼顧精確性和便利性的顯式關系式，新關系式的RMS明顯更小。