氣液兩相流擺動特性實驗的數(shù)值模擬

馬 霞，李國棟，范翠香，高 揚

西安理工大學(xué)，陜西 西安 710048

氣液兩相流擺動特性實驗的數(shù)值模擬

馬 霞，李國棟，范翠香，高 揚

西安理工大學(xué)，陜西 西安 710048

為了給鼓泡塔反應(yīng)器設(shè)計提供依據(jù)，運用計算流體力學(xué)（CFD）軟件模擬了鼓泡塔氣液兩相流動態(tài)行為。采用雙歐拉法對鼓泡塔矩形反應(yīng)器內(nèi)不同曝氣量下氣液兩相流的擺動特性進(jìn)行了模擬考察，液相采用標(biāo)準(zhǔn)κ-ε紊流模型，氣相采用分散相零方程模型，分析了網(wǎng)格尺寸、時間步長以及相間作用力對模擬結(jié)果的影響，模擬的曝氣量為42.5～237 mL/s。結(jié)果表明，當(dāng)相間作用力僅考慮阻力時，氣液兩相流呈現(xiàn)周期性擺動規(guī)律；隨著氣流量的增加，氣泡羽流的擺動幅度和頻率增大，同時液體的氣含率也在增加；模擬的氣液兩相流擺動頻率數(shù)據(jù)與實驗值吻合較好，兩者的相對誤差為7.2%～12.9%。

鼓泡塔 氣液兩相流 擺動特性 雙歐拉κ-ε模型

鼓泡塔是一種氣液兩相流反應(yīng)器，氣體從塔底部噴入，形成分散的氣泡群，氣泡群在液體中上升，同時帶動液體向上運動，氣泡羽流兩側(cè)會產(chǎn)生向下的液體流動，流動現(xiàn)象復(fù)雜多變，已有大量關(guān)于此類流動問題的文獻(xiàn)發(fā)表[1-5]，其中矩形反應(yīng)器內(nèi)氣泡群流動一直是實驗和數(shù)值模擬研究的熱點，但是由于流動特性與反應(yīng)器結(jié)構(gòu)及操作方式有很大的關(guān)系，此類反應(yīng)器的設(shè)計和放大仍然是一個問題。

與傳統(tǒng)的實驗方法相比，數(shù)值模擬具有耗費小、可重復(fù)等優(yōu)點，能夠預(yù)測實驗難以測量的局部流場細(xì)節(jié)。目前氣液兩相流模擬主要有歐拉-拉格朗日法[6-8]和歐拉-歐拉法[9-14]（雙歐拉法）兩種方法，歐拉-歐拉法的計算不受氣含率的限制，計算資源相對低[14]。利用雙歐拉κ-ε模型進(jìn)行三維模擬可以復(fù)現(xiàn)氣泡羽流的周期性運動特性[1,11-17]，但對于相間作用力仍存在爭議，Sokolichin等[11,12,15]認(rèn)為，阻力對氣液兩相流的擺動特性影響較大，虛擬質(zhì)量力對模擬結(jié)果的影響不明顯；Mudde等[14]報道，為了準(zhǔn)確模擬氣泡羽流的振蕩周期和速度場，還應(yīng)考慮虛擬質(zhì)量力與阻力；并且，不同的學(xué)者采用了不同的阻力模型。本研究采用計算流體力學(xué)（ANSYS14.0中CFX5）軟件，采用歐拉-歐拉法結(jié)合標(biāo)準(zhǔn)κ-ε紊流模型，基于Enrique Juliá等[18]的實驗結(jié)果，對鼓泡塔內(nèi)氣液兩相流的擺動特性進(jìn)行三維數(shù)值模擬，考察不同的網(wǎng)格尺寸、時間步長、相間作用力（阻力、升力和虛擬質(zhì)量力）的影響，選擇出適宜的參數(shù)和模型以能夠較為準(zhǔn)確地模擬鼓泡塔內(nèi)氣液兩相流的擺動特性，以期為提高鼓泡塔反應(yīng)器混合效率提供依據(jù)。

1 反應(yīng)器結(jié)構(gòu)及模型

基于Enrique Juliá的實驗結(jié)果[18]，建立圖1所示反應(yīng)器計算模型，其中鼓泡塔寬度W（沿X方向）為264 mm、厚度B（沿Z方向）為31 mm、高度H（沿Y方向）為792 mm，坐標(biāo)原點位于模型底部

式中，σ為湍流Prandtl數(shù)。

(4) 相間動量傳遞模型：

兩相之間的作用力可以表示為：

式中，CD為阻力系數(shù)，阻力系數(shù)分別采用CD為0.66[12]和Grace阻力模型：

式中，Vb為氣泡的最終速度，db為氣泡直徑，M為莫頓數(shù)，E0為奧托斯數(shù)。μref為0.000 9 kg/(m·s)。

式中，CL為升力系數(shù)；CV為虛擬質(zhì)量力系數(shù)。

(5) 初始條件和邊界條件

本研究采用計算流體力學(xué)（ANSYS14.0中CFX5）軟件模擬求解。

邊界條件設(shè)置：液相（水）的進(jìn)口速度為0 m/s，空氣的進(jìn)口速度是氣體流量的函數(shù)；出口邊界為自由液面邊界條件，氣相：允許逸出，液相：認(rèn)為是自由滑移的壁面，不允許流出；壁面：液相采用無滑移邊界條件，氣相采用自由滑移邊界條件。

初始計算條件：t為0 s時刻計算區(qū)域內(nèi)充滿水，即設(shè)水的體積分?jǐn)?shù)為1，空氣的體積分?jǐn)?shù)為0；初場內(nèi)水靜止不動，速度為0 m/s。

進(jìn)口氣流量Q0為42.5，82.5，125，237 mL/s，根據(jù)文獻(xiàn)[18]，對應(yīng)的當(dāng)量氣泡直徑分別為3.8，4.3，5，6.1 mm。水動力方程組采用有限體積法進(jìn)行求解。離散格式采用混合格式。

3 結(jié)果與討論

3.1 網(wǎng)格尺寸與時間步長的選擇

隨著網(wǎng)格尺寸變小，軸向速度的時均模擬值和實驗數(shù)據(jù)出現(xiàn)了較大離散[1,3,12]。為了研究網(wǎng)格尺寸對求解結(jié)果的影響，在如下計算工況計算了網(wǎng)格的影響（見表1）：時間步長t為0.025 s，進(jìn)口氣流量Q0為237 mL/s（相應(yīng)地氣體表觀速度為2.9 cm/s），計算當(dāng)量氣泡直徑D為6.1 mm。從表1可以看出，粗、中、細(xì)三種網(wǎng)格下氣泡羽流擺動頻率f的模擬結(jié)果與實驗結(jié)果的相對誤差分別為7.2%、8.8%、12.8%，在點A1（0, 0.2ymax, 0）位置處液相垂向速度的時均模擬值與實驗值的相對誤差分別為2.4%、15%、12.7%。與粗網(wǎng)格相比，中、細(xì)網(wǎng)格下氣泡羽流的擺動頻率、時間平均的液相軸向速度的模擬值表現(xiàn)出相對誤差增加的趨勢，但從總體上三者差異不大?？紤]到粗網(wǎng)格下模擬結(jié)果與實驗結(jié)果吻合得最好，除非特別指出，以下研究的網(wǎng)格均為粗網(wǎng)格。

表1 網(wǎng)格尺寸對計算氣體氣含率、擺動頻率和液相垂向速度的影響Table 1 Effect of mesh size on the calculated plume gas hold-up, oscillating frequency and vertical liquid velocity

時間步長對模擬結(jié)果的影響見表2，時間步長為0.25，0.025，0.01 s。時間步長與CFL數(shù)（柯朗數(shù)）有關(guān)，即：其中v是氣相垂向速度，Δy是單元尺寸；Δt是時間步長。考慮到單元尺寸為粗網(wǎng)格尺寸，且氣泡最終的上升速度經(jīng)驗值為0.3 m/s，則依據(jù)CFL準(zhǔn)則可得出：Δt≤0.036 s。因此，根據(jù)表2的計算結(jié)果，計算步長最終選為0.025 s，一方面Δt為0.025 s證明了CFL準(zhǔn)則，另一方面，更小的時間步長對模擬結(jié)果的影響不明顯。

表2 時間步長對計算氣體氣含率、擺動頻率和液相垂向速度參數(shù)的影響Table 2 Effect of time-step on the calculated plume gas hold-up , oscillating frequency and vertical liquid velocity

3.2 阻力、虛擬力和升力的作用

為了研究相間作用的影響，模擬了不同作用力的影響。工況A僅考慮阻力的影響，阻力系數(shù)采用常阻力系數(shù)Cd為0.66和變阻力系數(shù)Grace模型，計算結(jié)果見表3。氣泡尺寸分布的實驗情況表明[18]，氣泡尺寸均勻，且呈橢球狀，尺寸范圍較窄，大多數(shù)氣泡直徑為4.3～5.1 mm，氣泡沒有發(fā)生聚并和破碎。Grace阻力模型考慮到了氣泡的變形，阻力系數(shù)隨著氣泡的形狀而變化?？紤]到實驗中氣泡直徑范圍窄且變形小，可近似認(rèn)為橢球體，阻力系數(shù)按常值考慮更為合理，氣泡羽流的擺動頻率計算結(jié)果也證實了這一點。因此在以下研究中阻力系數(shù)按Cd為0.66計算。

表3 不同的阻力模型對計算氣體空隙率和擺動頻率參數(shù)的影響Table 3 Effect of different drag force model on the calculated gas hold-up and plume oscillating frequency parameters

工況B考慮阻力和虛擬質(zhì)量力的影響，當(dāng)虛擬質(zhì)量力系數(shù)為0.5時，B工況下氣泡羽流擺動頻率的計算值為0.133 Hz。在本研究的進(jìn)口氣流量范圍內(nèi)，相間作用力僅考慮阻力項與考慮阻力和虛擬質(zhì)量力相比，兩者對氣泡羽流擺動頻率的影響差別很小，原因在于由虛擬質(zhì)量力產(chǎn)生的流體加速運動是很小的，可忽略。

工況C考慮阻力和升力的影響，升力系數(shù)取為0.5，對于較高的進(jìn)口氣流量時，工況C會使氣泡羽流擺動頻率的模擬值（f=0.162 Hz）過高估計。事實上，在這種情況下，因為進(jìn)口速度高，使得液相產(chǎn)生較高的剪切率，升力的作用尤其顯得重要。數(shù)值模擬結(jié)果與實驗結(jié)果不一致，主要是升力系數(shù)取值可能不合理，有關(guān)這方面的研究不充分導(dǎo)致。

基于以上分析，本文研究中相間作用力僅考慮阻力，不考慮升力、虛擬質(zhì)量力的影響。

3.3 瞬態(tài)模擬結(jié)果

進(jìn)口氣流量Q0為237 mL/s時xy（z = 0 m）平面上氣相體積分?jǐn)?shù)與液相速度矢量分布瞬時圖見圖2。圖2(a)～(d)顯示了從90～97.3 s時氣泡羽流擺動的一個周期，約7.3 s，呈現(xiàn)了低頻振蕩特性，這是由垂向流動的不穩(wěn)定的旋渦結(jié)構(gòu)誘發(fā)的。氣泡群從底部曝氣區(qū)釋放出來，沿著區(qū)域中心線緩慢上升，周邊的液體不斷地被卷吸到羽流中，在卷吸作用下氣泡向中心聚集，形成帶有蘑菇云的氣泡柱，在氣泡柱兩側(cè)液體向下流動，形成左、右兩旋渦環(huán)流，受氣泡尺寸、氣泡上升速度、周圍壁面等影響，卷吸作用不對稱，迫使羽流兩側(cè)壓差不相等，使得羽流偏離原平衡位置，出現(xiàn)左右旋渦上下交錯、垂向流動，旋渦結(jié)構(gòu)不穩(wěn)定，氣泡柱呈蛇形擺動，從而使氣泡羽流發(fā)生了低頻振蕩。

圖2 氣相體積分?jǐn)?shù)云圖與液相速度矢量分布Fig.2 Distribution of gas hold-up represented by contours and the vector distribution of vertical liquid velocity represented by arrows

3.4 擺動頻率

圖3給出了Q0為237 mL/s時液相垂向速度在位置點A（0.5xmax, 0.2ymax, 0）處隨時間變化的模擬結(jié)果與實驗結(jié)果曲線。當(dāng)Q0為237 mL/s時，液相垂向速度的振幅模擬值在±0.6 m/s內(nèi)，數(shù)值模擬結(jié)果與實驗值基本一致。采用快速傅立葉變換對圖3(a)中數(shù)據(jù)進(jìn)行處理，相應(yīng)的頻譜特性曲線見圖3(b)：對于模擬數(shù)據(jù)，在f為0.131 Hz處達(dá)到峰值，即該羽流的擺動頻率為0.131 Hz，擺動周期為1/f為7.63 s；對于實驗數(shù)據(jù)，在f為0.125 Hz處達(dá)到峰值。

圖3 點A（0.5xmax,0.2ymax,0）相應(yīng)的液相垂向速度隨時間的變化和頻譜特性Fig.3 Time series of the horizontal liquid velocity and corresponding frequency series of power spectral density function for point A(0.5xmax,0.2ymax,0)

其他工況下氣泡羽流的擺動頻率也可以采用上面的方法獲得，Q0為125，82.5，42.5 mL/s時位置點A（0.5xmax, 0.2ymax, 0）處液相的垂向速度隨時間變化的模擬結(jié)果和頻譜特性分別見如圖4所示。Q0為125，82.5，42.5 mL/s時液相垂向速度的振幅模擬值分別在±0.4，±0.3，±0.2 m/s內(nèi)，隨著曝氣量的減少，氣泡羽流的擺動幅度在減小。

圖4 不同Q0下點A相應(yīng)的液相垂向速度隨時間的變化和頻譜特性Fig.4 Time series of the horizontal liquid velocity and corresponding frequency series of power spectral density function for point A under different Q0values

Enrique Juliá[18]通過對點A（0.5xmax, 0.2ymax, 0）、B（0.2xmax, 0.2ymax, 0）、C（-0.2xmax, 0.2ymax, 0）、D（-0.5xmax, 0.2ymax, 0）四個位置處的氣泡羽流擺動頻率進(jìn)行平均化處理，給出了平均的擺動頻率數(shù)據(jù)。為了和實驗數(shù)據(jù)更好的對比，四種工況下擺動頻率模擬數(shù)據(jù)也采用對點A（0.5xmax, 0.2ymax, 0）、B（0.2xmax, 0.2ymax, 0）、C（-0.2xmax, 0.2ymax, 0）、D（-0.5xmax, 0.2ymax, 0）四個位置處的擺動頻率值進(jìn)行平均化處理獲得，對比結(jié)果見表4，表中結(jié)果顯示：對于同一工況，四個位置處氣泡羽流的擺動頻率基本是一樣的，氣泡羽流擺動頻率的計算值與實驗值相對誤差為7.2%～12.9%。

表4 氣泡羽流的擺動頻率與實驗數(shù)據(jù)的比較Table 4 Comparison of experimental and calculated plume oscillating frequency

不同表觀氣速Ug（Ug為0.52～2.90 cm/s，對應(yīng)的氣流量為42.5～237 mL/s）下氣泡羽流的擺動頻率見圖5，隨著氣速增加，氣泡羽流的擺動頻率增加。對不同表觀氣速（0.52～2.90 cm/s）下總氣含率在100～200 s內(nèi)進(jìn)行時均處理，可得到時均總氣含率變化規(guī)律，見圖6，隨著氣速增加，液體的氣含率在增加。

圖5 不同表觀氣速時的擺動頻率Fig.5 Oscillating frequency under different superficial gas velocities

圖6 不同表觀氣速時的總氣含率Fig.6 Global gas hold-up under different superficial gas velocities

隨著氣流速度的增加（曝氣量的增大），流場中心區(qū)域的速度變大，流場的紊動由中心區(qū)域為主擴(kuò)散到了中心區(qū)域和氣泡柱的兩側(cè)區(qū)域，紊動強度增大，同時氣泡尺寸總體增加，也加劇了流場的紊動，由此導(dǎo)致了氣泡羽流的擺動幅度加大，擺動頻率加快；同時，環(huán)流渦旋結(jié)構(gòu)的發(fā)展更充分，渦兩側(cè)的氣相速度變大，更多的氣泡被帶入循環(huán)中，液體的氣含率增加。

4 結(jié) 論

a）本研究采用雙歐拉法結(jié)合標(biāo)準(zhǔn)的κ-ε方程模型，對鼓泡塔內(nèi)氣液兩相流動的動態(tài)行為進(jìn)行了三維數(shù)值模擬，能預(yù)測氣液兩相流的擺動特性，獲得了氣/液相速度、氣含率、擺動頻率等數(shù)據(jù)，模擬結(jié)果與實驗結(jié)果基本一致。

b）相間作用力模型影響氣泡流擺動頻率的計算結(jié)果表明：僅考慮阻力的影響，模擬的氣泡羽流擺動頻率與實驗結(jié)果吻合較好，兩者的相對誤差為7.2%～12.9%；虛擬質(zhì)量力的影響相對較小，可忽略；升力系數(shù)為0.5時，與實驗值相比，氣泡羽流擺動頻率的模擬值偏大，升力系數(shù)如何合理取值有待進(jìn)一步研究。

c）相間作用力僅考慮阻力的影響下模擬結(jié)果表明：隨著曝氣量的增加（42.5～237 mL/s），氣泡羽流的擺動頻率加快，由0.061 Hz增至0.134 Hz，擺動幅度加大，由0.2 m/s增至0.6m/s；同時液體的氣含率也相應(yīng)地增加了。

[1] Becker S, Sokolichin A, Eigenberger G. Gas-liquid flow in bubble columns and loop reactors. Part II. Comparison of detailed experiments and flow simulations [J]. Chem Eng Sci, 1994, 49(24):5747-5762.

[2] Buwa V V, Ranade V V. Characterization of dynamics of gas-liquid flows in rectangular bubble columns [J]. AIChE Journal, 2004, 50(10): 2394-2407

[3] Cheng W, Yuichi M, Fujio Y. Estimation of the liquid velocity field in two-phase flows using inverse analysis and particle tracking velocimetry [J]. Flow Measurement and Instrumentation, 2005, 16(5):303-308.

[4] 宋 策, 程 文, 胡保衛(wèi), 等. 氣泡羽流空隙率的計算及其不穩(wěn)定規(guī)律的研究 [J]. 水利學(xué)報, 2011, 42(4):419-424. Song Ce, Cheng Wen, Hu Baowei, et al. Research on the calculation of void fraction of bubble plume and its instability pattern [J]. Journal of Hydraulic Engineering, 2011, 42(4):419-424.

[5] 肖柏青, 張法星, 劉春艷, 等. 曝氣池內(nèi)氣泡羽流附壁效應(yīng)的試驗研究 [J]. 水力發(fā)電學(xué)報, 2012, 31(4):104-107,112. Xiao Boqing, Zhang Faxing, Liu Chunyan, et al. Canada effect of bubble plume in aeration tanks [J]. Journal of Hydroelectric Engineering, 2012, 31(4):104-107,112.

[6] Delnoij E, Lammers F A, Kuipers J A M, et al. Dynamic simulation of dispersed gas-liquid two-phase flow using a discrete bubble model [J]. Chem Eng Sci, 1997, 52(9):1429-1458.

[7] Darmana D, Deen N G, Kuipers J A M. Detailed modeling of hydrodynamics, mass transfer and chemical reactions in a bubble column using a discrete bubble model [J]. Chem Eng Sci, 2005, 60(12):3383-3404.

[8] Buwa V V, Deo D S, Ranade V V, et al. Eulerian-Lagrangian simulations of unsteady gas-liquid flows in bubble columns [J]. Int J Multiphase Flow, 2006, 32(7):864-885.

[9] Borchers O, Busch C, Sokolichin A, et al. Applicability of the standard k-epsilon turbulence model to the dynamic simulation of bubble columns. Part II. Comparison of detailed experiments and flow simulations [J]. Chem Eng Sci, 1999, 54(24):5927-5935.

[10] Bech K. Dynamic simulation of a 2D bubble column [J]. Chem Eng Sci, 2005, 60(19):5294-5304.

[11] Sokolichin A, Eigenberger G. Applicability of the standard κ-ε turbulence model to the dynamic simulation of bubble columns. Part I. detailed numerical simulations [J]. Chem Eng Sci, 1999, 54(13-14):2273-2284

[12] Pfleger D, Gomes S, Gilbert N, et al. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian-Eulerian modeling approach [J]. Chem Eng Sci, 1999, 54(21):5091-5099.

[13] Krishna R, Van Baten J M, Urseanu M I. Three-phase Eulerian simulations of bubble column reactors operating in the churn-turbulent regime: a scale up strategy [J]. Chem Eng Sci, 2000, 55(16):3275-3286.

[14] Mudde R F, Simonin O. Two- and three-dimensional simulations of a bubble plume using a two-fluid model [J]. Chem Eng Sci, 1999, 54(21):5061-5069.

[15] Oey R S, Mudde R F, Van H E A. Sensitivity study on interfacial closure laws in two-fluid bubbly flow simulations [J]. AIChE Journal, 2003, 49(7): 1621-1636.

[16] Buwa V V, Ranade V V. Dynamics of gas-liquid flow in a rectangular bubble column: experiments and single/multi-group CFD simulations [J]. Chem Eng Sci, 2002, 57(s 22-23):4715-4736.

[17] 馬 霞, 李國棟, 高 揚, 等. 不同曝氣量和縱橫比下針孔噴射氣泡羽流擺動特性的數(shù)值模擬研究 [J]. 水動力學(xué)研究與進(jìn)展, 2016, 31(4):433-440. Ma Xia, Li Guodong, Gao Yang, et al. Numerical simulation of the oscillating characterization of needle sparging bubble plume under different aeration rates and aspect ratios[J]. Chinese Journal of Hydrodynamics, 2016, 31(4):433-440.

[18] Enrique Juliá J, Hemandez L, Chiva S, et al. Hydrodynamic characterization of a needle sparger rectangular bubble column: Homogeneous flow, static bubble plume and oscillating bubble plume [J]. Chem Eng Sci, 2007, 62(22):6361-6377.

[19] David C W. Turbulence Modeling for CFD [M]. DCW Industries, 2006:1-522.

[20] Jakobsen H A, Sannaes B H, Grevskott S, et al. Modeling of vertical bubble-driven flow [J]. Industrakl and Engineering Chemistry Research, 1997, 36(10):4052-4074.

Numerical Simulation of the Oscillating Gas-Liquid Flow Experiment

Ma Xia, Li Guodong, Fan Cuixiang, Gao Yang
Xi’an University of Technology, Xi’an 710048, China

The numerical simulation of computational fluid dynamics (CFD) software was used to simulate the dynamic behavior of gas liquid two-phase flow in the bubble column reactor to provide the basis for the design of bubbling tower. Simulation study by Eulerian-Eulerian approach has been made for the oscillation characteristics of gas-liquid two-phase flow in different aeration units in a rectangular reactor. The influence of grid size, time step and interphase force on the simulation was analyzed using the standard k-ε turbulence model for liquid phase and the dispersed phase zero equation model for gas phase. The gas flow rates ranging from 42.5 to 237 mL /s were used for simulation. The results revealed that the gas-liquid two-phase flow exhibited the periodic oscillation law if resistance was considered as the only interphase force. With the increase of air flow, the amplitude and frequency of the bubble plume increased, and the gas hold up increased. The simulated frequency of gas-liquid two-phase bubble flow oscillation was in good agreement with the experimental results, and the relative error was 7.2%-12.9%.

bubble column; gas-liquid flow; the oscillating characteristic; Eulerian-Eulerian κ-ε model

O359

A

1001—7631 ( 2017 ) 01—0073—09

10.11730/j.issn.1001-7631.2017.01.0073.09

2016-12-14;

2017-01-28。

馬 霞（1976—），女，高級工程師。E-mail:maxia@xaut.edu.cn。

國家自然科學(xué)基金項目（31500340）；環(huán)境工程國家重點學(xué)科培育學(xué)科項目（106-5X1204）。