張 妍楊 帆,2

（1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；

2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

單個(gè)氣泡上浮過程的數(shù)值模擬

張 妍1楊 帆1,2

（1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；

2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

為了了解氣泡上浮過程中形態(tài)的變化，采用FLUENT軟件中“流體體積”（Volume of Fluid，VOF）模型對(duì)單個(gè)氣泡在靜止液體中的上浮情況進(jìn)行數(shù)值模擬，得到了奧特斯數(shù)（Eo）在O(10-1)～O(102)，莫頓數(shù)(Mo)在O(10-9)～O(104)范圍內(nèi)氣泡的形態(tài)。將計(jì)算得到的結(jié)果與氣泡形狀圖譜做對(duì)比，印證了無量綱參數(shù)Eo數(shù)、Mo數(shù)和Re數(shù)的數(shù)值大小與氣泡在上浮過程中形狀的變化和最終速度密切相關(guān)。

氣泡上浮 VOF 數(shù)值模擬 氣泡形狀

引言

氣液兩相流廣泛存在于人類的生產(chǎn)、生活等各個(gè)領(lǐng)域。其中，氣泡的產(chǎn)生及特性對(duì)船舶運(yùn)輸、石化工業(yè)、食品、醫(yī)療以及能源的開發(fā)和利用都有著重要的影響。對(duì)氣泡的形狀和運(yùn)動(dòng)特性的掌握，對(duì)生產(chǎn)過程中的參數(shù)設(shè)定、控制運(yùn)行以及生產(chǎn)效率的提高等方面起著重要的作用。

國(guó)內(nèi)外學(xué)者對(duì)氣泡上浮過程中的形態(tài)變化、速度變化、受力情況、破碎現(xiàn)象以及氣泡融合等現(xiàn)象做了大量研究。1959年，Young等[1]采用線性模型模擬研究了小雷諾數(shù)下氣泡和水滴的移動(dòng)。1976年，Grace等[2]人給出了著名的氣泡形狀圖譜，提出了三個(gè)無量綱數(shù)奧特斯數(shù)（Eo）、莫頓數(shù)（Mo）、雷諾數(shù)（Re），以表征氣泡在黏性溶液中的上升形態(tài)。Bhaga和Weber[3]接著完善了氣泡圖譜。此外，Li等[4]運(yùn)用VOF（Volume of Fluid）法模擬了氣泡在液體中的形成和上升。付攀等[5]對(duì)微小氣泡的運(yùn)動(dòng)進(jìn)行了仿真計(jì)算，并推導(dǎo)了運(yùn)動(dòng)方程。除了模擬方法之外，高速攝影也是主要采用的研究方式。例如，張建生等[6]運(yùn)用高速攝影技術(shù)研究了氣泡在水中的運(yùn)動(dòng)規(guī)律；郭容等[7]采用高速CCD成像技術(shù)實(shí)驗(yàn)，研究了氣泡在不同黏度溶液中的上升速度和形狀變化。

氣液兩相流中對(duì)相界面的追蹤是研究重點(diǎn)。目前，追蹤界面的方法主要有VOF法、level set法和Front Tracking等。其中，VOF方法是Hirt和Nichols[8]于1981年提出的。此種方法是在流場(chǎng)網(wǎng)格中定義某一種流體的相函數(shù)F，若其值為1，則說明此種流體完全占據(jù)網(wǎng)格，流體體積與網(wǎng)格體積之比為1；若其值為0，則表示網(wǎng)格內(nèi)沒有此種流體。相界面處網(wǎng)格內(nèi)有多相流體，因此取值為（0,1）。確定相函數(shù)F值為0到1的網(wǎng)格位置，即可實(shí)現(xiàn)對(duì)相界面的追蹤。VOF法計(jì)算量小、模擬精度高且容易實(shí)現(xiàn)，所以本文采用VOF方法來模擬靜止黏性液體中單個(gè)氣泡上升過程中形態(tài)變化等特性。

1 計(jì)算模型與方法

1.1 控制方程

本研究首先做以下幾點(diǎn)假設(shè)：（1）所模擬多相流等溫絕熱，不與外界發(fā)生熱交換；（2）氣體和液體均是不可壓縮牛頓流體；（3）由界面所分隔的兩相流體為單相的流體系。此時(shí)，不可壓縮黏性流體的控制方程為：

式中，u為流體的速度；ρ為流體密度；t為時(shí)間；p為壓強(qiáng)；μ為流體的動(dòng)力黏度；Fs為流體所受表面張力。

根據(jù)質(zhì)量守恒，得到相界面的運(yùn)動(dòng)所滿足輸運(yùn)方程為：

1.2 計(jì)算區(qū)域及網(wǎng)格

使用CFD前處理軟件ANSYS ICEM，建立模擬模型為底邊長(zhǎng)0.08m，高0.2m的矩形區(qū)域。經(jīng)拓?fù)浜蠼⑼暾哪Ｐ?，模型兩邊為固壁，上邊界和下邊界為周期性邊界條件。整個(gè)矩形區(qū)域內(nèi)充滿液體，液體內(nèi)部有一直徑為0.02m的圓形氣泡。圖1為模擬計(jì)算模型示意圖。

圖1 計(jì)算模型示意圖

圖2 網(wǎng)格示意圖

計(jì)算區(qū)域采用四邊形單元進(jìn)行離散，允許的最大網(wǎng)格尺寸為0.08。網(wǎng)格節(jié)點(diǎn)數(shù)為25219，網(wǎng)格數(shù)為25010。圖2所示為網(wǎng)格示意圖。

1.3 計(jì)算方法

模擬采用非穩(wěn)態(tài)基于壓力的VOF兩相模型進(jìn)行計(jì)算。重力加速度設(shè)為鉛錘向下的9.8m/s2。采用PISO算法實(shí)現(xiàn)壓力與速度的耦合。非穩(wěn)態(tài)項(xiàng)采用一階隱式時(shí)間推進(jìn)法處理，計(jì)算時(shí)間步長(zhǎng)為0.1×10-3s，時(shí)間步數(shù)為20000。

1.4 氣泡運(yùn)動(dòng)概述

氣泡在上升過程中主要受到重力、曳力、表面張力等各種作用力，其綜合作用決定了氣泡的形狀。各作用力的相對(duì)大小可以用Eo數(shù)、Mo數(shù)和Re數(shù)表征。

Eo數(shù)表征重力與表面張力的比值，定義如下：

其中，g為重力加速度，m/s2；Δρ為液體與氣體的密度差，kg/m3；dc為氣泡當(dāng)量直徑，m；σ為表面張力系數(shù)，N/m。

Mo數(shù)表征連續(xù)相的物理性質(zhì)，突出了粘度的影響。定義如下：

其中，lμ為液體的動(dòng)力粘度，Pa·s；lρ為液體的密度，kg/m3。

Re數(shù)表征慣性力與粘性力的比值，定義如下：

其中，v為氣泡穩(wěn)定后的最終速度，m/s。

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

本文首先選取了五組Eo數(shù)和Mo數(shù)，給定一個(gè)密度差分別計(jì)算了各對(duì)應(yīng)參數(shù)，并根據(jù)計(jì)算所獲得的氣泡最終速度，確定該工況所對(duì)應(yīng)的Re數(shù)。計(jì)算參數(shù)如表1所示。

表1 計(jì)算工況

2.1 A4工況下氣泡上升過程中形態(tài)的變化

圖3（a）-（f）表示的是A4工況（即Eo＝5.21×101，Mo＝1.45×10-4）氣泡上升中幾個(gè)典型的形狀變化圖。在質(zhì)量力作用下，氣泡在上浮初期如圖3（b）所示，由下邊界中心點(diǎn)處向內(nèi)凹陷，凹陷區(qū)域隨著上浮過程不斷增大，氣泡形狀越來越扁。從圖4（e）可看出，氣泡上升到一定程度后，凹陷處開始變平坦，氣泡稍微變厚直至氣泡穩(wěn)定，最終保持球帽形勻速上升。這一過程是氣泡所受各力平衡的過程。氣泡所受的浮力、重力、曳力和表面張力是影響氣泡上升過程中形狀變化和速度的重要因素。所受重力基本不變，其他三個(gè)力隨著氣泡上浮的位置、速度的大小等逐漸變化而達(dá)到平衡。Mo數(shù)較小（黏度較?。?、Eo數(shù)較大（表面張力較?。┒y以維持氣泡的初始圓形形狀，最終變成下邊界水平的球帽形。

2.2 A4工況下氣泡上升過程中的流線變化

圖3 A4工況氣泡上浮形狀變化圖

圖4中（a）-（f）為A4工況下氣泡上浮流線變化過程。為便于觀察流線與氣泡的相對(duì)位置，圖中只標(biāo)示出左半部分，而右半部分與之對(duì)稱。由圖可直觀觀察到氣泡從靜止開始逐漸變形過程中氣泡的回流現(xiàn)象。

如圖4所示，（a）、（b）圖為氣泡運(yùn)動(dòng)初期，周圍很快形成渦旋。上邊界和下邊界形成一個(gè)壓力差而推動(dòng)氣泡上浮。隨著氣泡的上升渦旋不斷增大，氣泡下部向內(nèi)凹陷程度越大。到達(dá)一定程度后，回流渦旋不再增大，氣泡下邊界逐漸平坦，最終氣泡形狀基本不變，以勻速垂直上浮。

圖4 A4工況氣泡上浮流線圖

2.3 計(jì)算結(jié)果在圖譜上的分布

圖5為氣泡形態(tài)圖譜簡(jiǎn)圖以及計(jì)算工況在圖譜上的分布。圖譜被劃分為三個(gè)區(qū)域，即球形區(qū)、橢球形區(qū)和球帽形區(qū)。所計(jì)算工況A1、A2分布在橢球形區(qū)，可以看出，Mo數(shù)和Eo數(shù)較小、Re數(shù)較大時(shí)（A2），較難維持球形而變得更扁。這是由于液體黏度小，阻止氣泡的變形能力較差，從而使得氣泡變得更扁。相反，液體黏度較大時(shí)，氣泡較易維持球形而更圓。工況A3分布在球形區(qū)，上升過程中始終維持球形而達(dá)到穩(wěn)定的速度。工況A4、A5分布在球帽形區(qū)域。與橢球形區(qū)域相似，Mo數(shù)和Eo數(shù)較小、Re數(shù)較大時(shí)（工況A4），更難維持球形而成為球帽形，氣泡最終下邊界變得平坦。而A5工況Mo數(shù)和Eo數(shù)較大，下邊界部分向內(nèi)凹陷，上部分基本維持球形而變成橢球帽形。由此可見，三個(gè)無量綱參數(shù)Eo數(shù)、Mo數(shù)和Re數(shù)對(duì)氣泡上浮形狀有著密切的關(guān)系。

圖5 氣泡形態(tài)圖譜簡(jiǎn)圖以及計(jì)算結(jié)果在圖譜上的分布圖

2.4 球帽形區(qū)內(nèi)氣泡的最終速度與形狀的變化關(guān)系

為了研究球帽形區(qū)域內(nèi)氣泡的最終形狀以及速度與Mo數(shù)的關(guān)系，固定Eo數(shù)不變，對(duì)通過改變液體黏度所得的五個(gè)Mo數(shù)工況點(diǎn)進(jìn)行計(jì)算，計(jì)算工況如表2所示。

表2 計(jì)算工況

圖6則為表2中所示計(jì)算工況下氣泡的最終速度和形狀隨Mo數(shù)變化圖。上述計(jì)算條件均在球帽形區(qū)域，只改變了液體黏度的大小，Mo數(shù)中其他參數(shù)保持不變?？梢钥闯觯S著Mo數(shù)變小，氣泡下表面向里凹陷的程度越大，最終難以維持基本的球形而變成球帽形。Mo數(shù)較大時(shí)，則凹陷程度越小，更接近于球形。結(jié)果表明，液體的黏度越大，Mo數(shù)越大，氣泡的變形程度越小，最終速度越小。高黏度液體更有利于阻止氣泡的變形和減慢運(yùn)動(dòng)速度。同時(shí)，據(jù)圖6所得，在橢球帽區(qū)域中，氣泡上浮過程中的最終速度與Mo數(shù)的對(duì)數(shù)基本成線性反比，即與黏度四次方的對(duì)數(shù)呈現(xiàn)反比的關(guān)系。最終速度及形狀隨Mo數(shù)變化圖

圖6 氣泡最終速度與形狀隨Mo數(shù)變化圖

3 結(jié)論

（1）經(jīng)過對(duì)氣泡上浮的二維數(shù)值模擬，驗(yàn)證了Eo數(shù)O(-1)～O(2)，Mo數(shù)從O(-9)～O(4)的流動(dòng)范圍中，氣泡最終形狀包括球形、橢球形、球帽形、橢球帽形等形態(tài)。這些形狀與Bhagar和Weber所得的氣泡圖譜相對(duì)應(yīng)。

（2）三個(gè)無量綱參數(shù)Eo數(shù)、Mo數(shù)和Re數(shù)對(duì)氣泡上浮形狀有著密切關(guān)系。

（3）在橢球形區(qū)域內(nèi)，Mo數(shù)和Eo數(shù)較小時(shí)，較難維持球形而更變得更扁。

（4）在橢球帽區(qū)域，隨著Mo數(shù)變小，氣泡下表面向里凹陷的程度越大，最終難以維持基本的球形而變成球帽形。

（5）高黏度液體更有利于阻止氣泡的變形和減慢氣泡上浮運(yùn)動(dòng)速度。在橢球帽區(qū)域中，氣泡上浮的最終速度與Mo數(shù)的對(duì)數(shù)基本成線性反比，即與黏度四次方的對(duì)數(shù)呈反比關(guān)系。

[1]Young N O，Goldstein J S，Block M J.The Motion of Bubbles in A Vertical Temperature Gradient[J].Journal of Fluid Mechanics，1959，（3）：350-356.

[2]Grace J R，Wairegi T，Nguyen T H.Shapes and Velocities of Single Drops and Bubbles Moving Freely Through Immiscible Liquids[J].Transactions of the Institution of Chemical Engineers，1976，（3）：167-173.

[3]Bhaga D，Weber M E.Bubbles in Viscous Liquids：Shapes，Wakes and Velocities[J].Journal of Fluid Mechanics，1981，（105）：61-85.

[4]Li Y，Yang G Q，Zhang J P，et al. Numerical Studies of Bubble Formation Dynamics in Gas-liquid-solid Fluidization at High Pressures[J].Powder Technology，2001，（2）：246-260.

[5]付攀，王路.水中微氣泡運(yùn)動(dòng)特性的理論研究與仿真[J].艦船電子工程，2009，（4）：155-158.

[6]張建生，呂青，孫傳東，等.高速攝影技術(shù)對(duì)水中氣泡運(yùn)動(dòng)規(guī)律的研究[J].光子學(xué)報(bào)，2000，（10）：952-955.

[7]郭容，蔡子琦，高正明.黏性流體中單氣泡的運(yùn)動(dòng)特性[J].高?；瘜W(xué)工程學(xué)報(bào)，2009，（6）：916-921.

[8]Hirt C W，Nichols B D.Volume of Fluid（VOF）Method for the Dynamics of Free Boundaries[J].Journal of Computational Physics，1981，（1）：201-225.

Numerical Simulation of Single Bubble Rising Process

ZHANG Yan1,YANG Fan1,2
(1.School of energy and power engineering, University of Shanghai for Science and Technology, Shanghai 200093; 2.Key Laboratory of multiphase flow and heat transfer of power engineering, University of Shanghai for Science and Technology,Shanghai 200093)

In order to changes in the process of understanding the bubbles form, using FLUENT software "fluid volume (volume of fluid, VOF model of a single bubble in a quiescent liquid float situation of numerical simulation by Otis number (EO) in O (- 1) -O (2), Morton number (MO) in O (- 9) to o (4) within the scope of the bubble morphology. The calculated results are compared with the bubble shape map, which verifies that the numerical values of the dimensionless parameters Eo number, Mo number and Re number are closely related to the change of the shape and the final velocity of the bubbles.

bubble floating, VOF, numerical simulation, bubble shape

上海市科學(xué)技術(shù)委員會(huì)科研計(jì)劃項(xiàng)目（13DZ2260900）資助。