槽道湍流內(nèi)氣泡瞬態(tài)受力數(shù)值研究

龐明軍，孫濤，代軍榮（常州大學(xué)機(jī)械工程學(xué)院，江蘇 常州 213164）

?

研究開發(fā)

槽道湍流內(nèi)氣泡瞬態(tài)受力數(shù)值研究

龐明軍，孫濤，代軍榮
（常州大學(xué)機(jī)械工程學(xué)院，江蘇 常州 213164）

摘要：為了深入理解氣泡在湍流場(chǎng)內(nèi)的運(yùn)動(dòng)機(jī)制，使用歐拉-拉格朗日單向耦合數(shù)值方法，詳細(xì)分析了氣泡在低雷諾數(shù)槽道湍流場(chǎng)內(nèi)的瞬態(tài)受力情況；其中液相湍流速度場(chǎng)采用直接數(shù)值模擬方法求解，氣泡的瞬態(tài)受力由牛頓運(yùn)動(dòng)方程計(jì)算；計(jì)算時(shí)，考慮了相間阻力、剪切升力、壓力梯度力、虛擬質(zhì)量力、重力對(duì)氣泡運(yùn)動(dòng)的影響。目前的計(jì)算結(jié)果表明：氣泡所受的瞬態(tài)作用力分量同時(shí)取決于重力作用方向、液相流向以及氣泡所處的法向位置；不同方向上、不同位置處影響氣泡運(yùn)動(dòng)的主要作用力分量是不同的；相比較而言，與重力方向垂直的剪切升力分量在近壁面區(qū)域?yàn)橛绊憵馀葸\(yùn)動(dòng)的主要作用力，壓力梯度力的法向分量在近壁面區(qū)域之外為影響氣泡運(yùn)動(dòng)的主要作用力，相間阻力分量在整個(gè)槽道區(qū)域內(nèi)均為影響氣泡運(yùn)動(dòng)的主要作用力，除了豎直槽道近壁面處之外、虛擬質(zhì)量力也均為影響氣泡運(yùn)動(dòng)的主要作用力。

關(guān)鍵詞：歐拉-拉格朗日模型；單向耦合；直接數(shù)值模擬；氣泡；受力分析

鑒于良好的物理特性，泡狀流被廣泛應(yīng)用于金屬冶煉、礦物浮選、泡狀反應(yīng)堆、微氣泡減阻和流體輸送等工程領(lǐng)域。對(duì)于泡狀流而言，氣泡在液相中的分布會(huì)直接影響液相的流動(dòng)特性、壓降分布以及傳質(zhì)傳熱等特性，為了理解氣泡的分布模式和分布機(jī)理，國(guó)內(nèi)外學(xué)者開展了大量的研究[1-4]。然而由于問題的復(fù)雜性，目前關(guān)于氣泡的分布模式仍未取得一個(gè)普遍認(rèn)可的機(jī)理解釋。特別是從氣泡受力角度出發(fā)，研究氣泡運(yùn)動(dòng)和分布的研究者為數(shù)不多。陳鳳等[5]研究了電場(chǎng)力作用下沸騰氣泡的電場(chǎng)力情況，閆超星等[6]研究了擺動(dòng)情況下氣泡慣性力和浮力的變化情況。據(jù)本文作者所知，僅 GIUSTI等[7]給出了湍流場(chǎng)內(nèi)氣泡受力時(shí)均值的分布情況。事實(shí)上，氣泡的瞬態(tài)運(yùn)動(dòng)與其瞬態(tài)受力有關(guān)，研究氣泡在湍流場(chǎng)的瞬態(tài)受力情況，有助于了解處于湍流場(chǎng)不同位置的氣泡，支配其運(yùn)動(dòng)的主要相間力有哪些，也有助于理解氣泡是如何從湍流外層運(yùn)動(dòng)到湍流內(nèi)層。鑒于實(shí)驗(yàn)測(cè)量的困難性，本文使用數(shù)值模擬方法對(duì)槽道湍流內(nèi)氣泡的瞬態(tài)受力進(jìn)行了深入的分析和研究。

1 計(jì)算條件和數(shù)值方法

1.1 計(jì)算條件

因?yàn)橄嚅g作用力（如相間阻力和剪切升力）與氣泡的滑移速度密切相關(guān)，而氣泡的滑移速度又與氣泡所受浮力有關(guān)，因此，為了理解重力在槽道內(nèi)作用方向?qū)馀菟矐B(tài)受力的影響，本文設(shè)計(jì)了兩大類計(jì)算工況：一類是讓氣泡浮力方向與液相主流方向一致；另一類是讓氣泡浮力方向與液相主流方向垂直。同時(shí)為了研究氣泡位置對(duì)氣泡瞬態(tài)受力的影響，每一大類計(jì)算工況又分別包含了3種計(jì)算工況，即讓氣泡分別位于黏性底層（即氣泡位置y+=3.75）、過渡層（即氣泡位置 y+=10）和湍流外層（即氣泡位置 y+=145），以便觀察不同重力方向和不同氣泡位置氣泡瞬態(tài)受力情況。其中y+為量綱為1的壁面位置坐標(biāo)，氣液相的具體計(jì)算參數(shù)見表 1。目前使用的網(wǎng)格已在文獻(xiàn)[8-9]進(jìn)行了詳細(xì)的精度驗(yàn)證，本文不再贅述，表明網(wǎng)格精度能夠滿足計(jì)算要求。

計(jì)算區(qū)域?yàn)橐粋€(gè)二維槽道，幾何形狀和詳細(xì)尺寸見圖1，圖中x、y、z方向分別對(duì)應(yīng)液相的流向、法向和展向。如圖所示，對(duì)于豎直槽道，重力方向與流向平行；當(dāng)計(jì)算水平槽道時(shí)，將重力方向逆時(shí)針旋轉(zhuǎn)90°與流向垂直。計(jì)算時(shí)，對(duì)于液相，在流向和展向設(shè)置周期邊界條件，在壁面法向設(shè)置無(wú)滑移邊界條件；對(duì)于氣泡，同樣在流向和展向設(shè)置周期邊界條件，為了更好地觀察氣泡在不同位置時(shí)瞬態(tài)受力情況，氣泡在法向的位置始終保持不變。

表1 氣液相計(jì)算參數(shù)

圖1 計(jì)算區(qū)域和尺寸

1.2 控制方程

目前采用歐拉-拉格朗日單向耦合模型模擬氣液兩相流動(dòng)，計(jì)算時(shí)不考慮氣泡對(duì)液相流動(dòng)的影響。為了準(zhǔn)確求解液相的瞬態(tài)速度場(chǎng)，液相流場(chǎng)采用直接數(shù)值模擬方法求解，氣泡的瞬態(tài)受力由牛頓運(yùn)動(dòng)方程進(jìn)行計(jì)算。計(jì)算時(shí)，所有的控制方程進(jìn)行了量綱歸一化處理，與長(zhǎng)度相關(guān)的量用半槽道高度（h）量綱為1化，與速度相關(guān)的量用摩擦速度（uτ）量綱為1化，與時(shí)間相關(guān)的量用h/uτ量綱為1化，具體過程詳見文獻(xiàn)[8]。

量綱歸一化液相的連續(xù)性和動(dòng)量方程見式(1)、式(2)。

量綱歸一化的氣泡運(yùn)動(dòng)方程見式(3)。

式(3)中，等號(hào)左邊項(xiàng)為氣泡慣性力項(xiàng)，等號(hào)右邊項(xiàng)分別為壓力梯度力、相間阻力、虛擬質(zhì)量力、剪切升力和浮力項(xiàng)；CD、CV和CLF分別為相間阻力系數(shù)、虛擬質(zhì)量力系數(shù)和升力系數(shù)。對(duì)于球形氣泡，虛擬質(zhì)量力系數(shù)為 0.5[10]，阻力系數(shù)和升力系數(shù)分別由式(4)[11]和式(5)[12]來計(jì)算。

在式(5)中

式中Srb、Reb分別為量綱歸一化的剪切率和氣泡雷諾數(shù)，ζ=SrbReb。

1.3 數(shù)值方法

目前的計(jì)算由自編程序完成，液相控制方程采用有限差分法進(jìn)行離散，空間項(xiàng)采用二階精度的差分格式；在流向和展向進(jìn)行均勻網(wǎng)格劃分，在法向進(jìn)行非均勻網(wǎng)格劃分，越靠壁面網(wǎng)格越密；為了避免壓力場(chǎng)的不穩(wěn)定性，采用交錯(cuò)網(wǎng)格系統(tǒng)，速度存儲(chǔ)于單元格的邊界上，而壓力存儲(chǔ)于網(wǎng)格中心；時(shí)間的推進(jìn)，其壓力項(xiàng)采用隱式格式，速度的時(shí)間推進(jìn)采用二階精度的Adams-Bashforth格式；速度和壓力的耦合采用MAC方法；壓力泊松方程由多重網(wǎng)格進(jìn)行求解。采用二階精度的 Crank-Nicholson格式對(duì)微氣泡的加速度方程進(jìn)行積分，得到微氣泡的速度；利用同樣的方法對(duì)速度方程積分得到氣泡流向和展向上的位移。有關(guān)計(jì)算方法和計(jì)算模型的詳細(xì)介紹和準(zhǔn)確性驗(yàn)證詳見文獻(xiàn)[8]，該作者已采用該模型成功模擬了微氣泡對(duì)液相湍流的調(diào)制情況。

2 結(jié)果分析

下面給出不同工況氣泡瞬態(tài)受力隨計(jì)算時(shí)間的變化情況，圖 2、圖 3中縱坐標(biāo)為氣泡瞬態(tài)受力情況，橫坐標(biāo)為量綱歸一化的計(jì)算時(shí)間[t*=t/(h/uτ)]。圖中負(fù)的瞬態(tài)力表示該力指向坐標(biāo)軸負(fù)方向。

2.1 豎直槽道氣泡瞬態(tài)受力分析

圖2給出豎直槽道內(nèi)氣泡瞬態(tài)受力情況，即氣泡所受浮力方向與液相流向一致情況下，氣泡在流向和法向上的受力情況，有關(guān)相間力的具體計(jì)算公式詳見式(3)。從圖2可以看出，因湍流場(chǎng)的不穩(wěn)定性，處于不同位置的氣泡所受相間作用力隨著計(jì)算時(shí)間的推移在作不穩(wěn)定的波動(dòng)，特別是處于過渡層和外層上的氣泡，其法向分力波動(dòng)幅度較大；而且氣泡所受相間力的大小與其所處位置密切相關(guān)。在流向上，隨著遠(yuǎn)離壁面，氣泡所受的相間阻力（FD）和虛擬質(zhì)量力（FVM）在增大；在法向上，隨著遠(yuǎn)離壁面，氣泡所受的剪切升力（FL）在減小，相間阻力（FD）是先增大后減小，而虛擬質(zhì)量力（FVM）和壓力梯度力（Fp）是在逐漸增大。

從圖2中可以清楚得看到，氣泡所處位置不同，不同方向上控制氣泡運(yùn)動(dòng)的相間力也不同。當(dāng)氣泡位于黏性底層時(shí)（y+=3.75），在流向上控制氣泡運(yùn)動(dòng)的力主要為浮力（FG）和相間阻力（FD），在法向上控制氣泡運(yùn)動(dòng)的力主要為剪切升力（FL）和相間阻力（FD）。當(dāng)氣泡位于過渡層時(shí)（y+=15），在流向上控制氣泡運(yùn)動(dòng)的力主要為浮力（FG）和相間阻力（FD），虛擬質(zhì)量力（FVM）也不容忽視；在法向上控制氣泡運(yùn)動(dòng)的力主要為壓力梯度力（Fp）和虛擬質(zhì)量力（FVM），但相間阻力（FD）也較為重要。氣泡位于湍流外層與氣泡處于過渡層相似，流向上控制氣泡運(yùn)動(dòng)的力也主要為浮力（FG）和相間阻力（FD），同樣虛擬質(zhì)量力（FVM）也較為重要；但在法向上，不同于氣泡位于過渡層時(shí)的情況，控制氣泡運(yùn)動(dòng)的力主要為虛擬質(zhì)量力（FVM）和壓力梯度力（Fp）。

圖2 豎直槽道內(nèi)氣泡瞬態(tài)受力情況

2.2 水平槽道氣泡瞬態(tài)受力分析

圖3給出了水平槽道內(nèi)氣泡瞬態(tài)受力情況，即液相主流方向與浮力（重力）方向相垂直情況下，氣泡的瞬態(tài)受力。類似于豎直槽道內(nèi)氣泡的受力情況，水平槽道內(nèi)氣泡的瞬態(tài)受力也隨著計(jì)算時(shí)間和氣泡位置在變化。當(dāng)氣泡位于黏性底層時(shí)，在流向上控制氣泡運(yùn)動(dòng)的力主要有阻力（FD）、剪切升力（FL）和虛擬質(zhì)量力（FVM）；在法向上控制氣泡運(yùn)動(dòng)的力主要為浮力（FG）和相間阻力（FD）。當(dāng)氣泡位于過渡層和湍流外層時(shí)，流向上控制氣泡運(yùn)動(dòng)的力主要為相間阻力（FD）和虛擬質(zhì)量力（FVM）；而法向上控制氣泡運(yùn)動(dòng)的力主要為浮力（FG）、相間阻力（FD）、虛擬質(zhì)量力（FVM）和壓力梯度力（FP），而且阻力的方向基本與浮力的方向相反。

2.3 槽道內(nèi)氣泡受力分析總結(jié)

為了清楚展示氣泡在槽道內(nèi)的受力情況，根據(jù)圖2和圖3給出了微氣泡受力分析示意圖，如圖4所示?？梢钥闯?，控制氣泡運(yùn)動(dòng)主要的力與重力（浮力）方向密切相關(guān)。當(dāng)氣泡位于黏性底層時(shí)，在重力方向上，控制氣泡運(yùn)動(dòng)主要的力為浮力（FG）和相間阻力（FD）；在與重力相垂直的方向上，控制氣泡運(yùn)動(dòng)的力主要為剪切升力（FL）和相間阻力（FD），對(duì)于水平槽道，虛擬質(zhì)量力（FVM）不容忽視。當(dāng)氣泡位于過渡層和湍流外層時(shí)，在重力方向上，控制氣泡運(yùn)動(dòng)的力主要有浮力（FG）、相間阻力（FD）和虛擬質(zhì)量力（FVM），對(duì)于水平槽道，壓力梯度力（Fp）也很重要；在與重力相垂直的方向上，重力的影響比較復(fù)雜，對(duì)于豎直槽道而言，控制氣泡運(yùn)動(dòng)的力主要為壓力梯度力（Fp）和虛擬質(zhì)量力（FVM），對(duì)于水平槽道而言，主要的力為阻力（FD）和虛擬質(zhì)量力（FVM）。

圖3 水平槽道內(nèi)氣泡瞬態(tài)受力情況

作為一個(gè)事實(shí)，重力）（浮力）的存在會(huì)使氣泡在重力方向上比其周圍的液體運(yùn)動(dòng)得快，即存在一個(gè)負(fù)的滑移速度，導(dǎo)致在重力方向上，氣泡所受相間阻力的方向始終與浮力方向相反，且二者的數(shù)值也很大；另外，滑移速度的存在，再加上局部液相的速度梯度，必然會(huì)在與重力相垂直的方向上產(chǎn)生一個(gè)剪切升力，因液相的局部速度梯度在壁面附近最大，所以氣泡在黏性底層時(shí)其所受的剪切升力較大，對(duì)其運(yùn)動(dòng)起主導(dǎo)作用。隨著遠(yuǎn)離壁面，液相的速度梯度急劇下降，剪切升力變得不再重要。虛擬質(zhì)量力是因氣泡運(yùn)動(dòng)拖動(dòng)局部液體一起運(yùn)動(dòng)而產(chǎn)生的一種力，所以對(duì)氣泡而言是一種附加阻力，因此方向也始終與浮力方向相反。

圖4 氣泡主要受力示意圖

對(duì)于充分發(fā)展的槽道湍流，當(dāng)流動(dòng)達(dá)到統(tǒng)計(jì)穩(wěn)態(tài)時(shí)，液相在流向上的速度梯度基本為零，對(duì)應(yīng)的短距離內(nèi)的壓力梯度也為零，因此在流向上氣泡的壓力梯度力可忽略不計(jì)；而液相在法向上的速度梯度較大，所以在法向上，氣泡所受壓力梯度力較大，且隨著流場(chǎng)的變動(dòng)在不斷調(diào)整和變化。

從目前的分析可以看出，不論是水平槽道還是豎直槽道，氣泡從槽道中心區(qū)域向壁面附近區(qū)域的運(yùn)動(dòng)與壓力梯度力和虛擬質(zhì)量力的法向分量有直接的關(guān)系。

3 結(jié) 論

利用歐拉-拉格朗日單向耦合數(shù)值方法詳細(xì)調(diào)查了豎直和水平槽道湍流場(chǎng)內(nèi)、不同法向位置處氣泡的瞬態(tài)受力情況，可以總結(jié)以下幾點(diǎn)結(jié)論。

（1）氣泡瞬態(tài)受力隨湍流場(chǎng)不斷變化，其大小和方向同時(shí)取決于重力方向、液相流向以及氣泡在槽道內(nèi)壁面法向上的位置；不同方向上、不同位置處，影響氣泡運(yùn)動(dòng)的主要作用力分量是不同的。

（2）與重力相垂直的剪切升力分量在近壁面區(qū)域處為控制氣泡運(yùn)動(dòng)的主要作用力，在近壁面之外區(qū)域、壓力梯度力法向分量為控制氣泡運(yùn)動(dòng)的主要作用力，不同方向上的相間阻力和虛擬質(zhì)量力分量均為控制氣泡運(yùn)動(dòng)的主要作用力。

（3）相同工況下，法向上的虛擬質(zhì)量力分量與法向上的壓力梯度力分量有關(guān)；重力方向上的相間阻力分量取決于氣泡所受浮力，非重力方向上的相間阻力在遠(yuǎn)離壁面區(qū)域處較小且隨湍流場(chǎng)波動(dòng)。

（4）氣泡從槽道中心區(qū)域向壁面附近區(qū)域的運(yùn)動(dòng)與壓力梯度力和虛擬質(zhì)量力的法向分量有直接的關(guān)系。

參考文獻(xiàn)

[1] SHEN X Z，MISHIMA K，NAKAMURA H.Two-phase phase distribution in a vertical large diameter pipe[J].Int.J.Heat Mass Transfer，2005，48（1）：211-225.

[2] OKAWA T，KATAOKA I，MORI M.Numerical simulation of lateral phase distribution in turbulent upward bubbly two-phase flow[J].Nucl.Eng.Des.，2002，213（2/3）：183-197.

[3] 劉國(guó)強(qiáng)，孫立成，田道貴，等.豎直圓管內(nèi)空泡份額徑向分布特性形成機(jī)制[J].核動(dòng)力工程，2014，35（1）：42-45.

[4] 龐明軍，徐一丹，魏進(jìn)家.管道泡狀流相分布模式和分布機(jī)理研究進(jìn)展[J].化工進(jìn)展，2014，33（11）：2829-2841.

[5] 陳鳳，宋耀祖，陳民.電場(chǎng)作用下的氣泡受力分析[J].工程熱物理學(xué)報(bào)，2005，26（6）：146-148.

[6] 閆超星，閻昌琪，孫立成.周期性力場(chǎng)作用下氣泡徑向附加力分析[J].原子能科學(xué)技術(shù)，2013，47（12）：2220-2224.

[7] GIUSTI A，LUCCI F，SOLDATI A.Influence of the lift force in direct numerical simulation of upward/downward turbulent channel flow laden with surfactant contaminated microbubbles[J].Chem.Eng.Sci.，2005，60（22）：6176-6187.

[8] PANG M J，WEI J J，YU B.Numerical study of bubbly upflows in a vertical channel using the Euler–Lagrange two-way model[J].Chem.Eng.Sci.，2010，65（23）：6215-6228.

[9] YU B，KAWAGUCHI Y.Direct numerical simulation of visco-elastic drag-reducing flow: a faithful finite difference method [J].Journal of Non-Newtonian Fluid Mechanics， 2004，116（2/3）：431-466.

[10] TABIB M V， ROY S A，JOSHI J B.CFD simulation of bubble column——an analysis of interphase forces and turbulence models [J].Chemical Engineering Journal，2008，139（8）：589-614.

[11] LAIN S， BRODER D， SOMMERFELD M， et al.Modelling hydrodynamics and turbulence in a bubble column using the Euler-Lagrange procedure [J].Int.J.Multiphase Flow，2002，28（8）：1381-1407.

[12] LEGENDER D， MAGNAUDET J.The lift force on a spherical bubble in a viscous linear shear flow [J].J.Fluid Mech.， 1998，369：81-126.

第一作者及聯(lián)系人：龐明軍（1976—），男，博士，副教授，研究方向?yàn)闅庖簝上嗔髋c湍流減阻。E-mail pangmj@cczu.edu.cn。

中圖分類號(hào)：TK 124

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1000-6613（2016）07-1980-06

DOI：10.16085/j.issn.1000-6613.2016.07.007

收稿日期：2015-11-18；修改稿日期：2015-12-14。

基金項(xiàng)目：國(guó)家自然科學(xué)基金（51376026）及江蘇省大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃重點(diǎn)項(xiàng)目（201510292016Z）。

Numerical investigation on instantaneous forces exerting on bubbles in channel turbulent flows

PANG Mingjun，SUN Tao，DAI Junrong
（School of Mechanical Engineering，Changzhou University，Changzhou 213164，Jiangsu，China）

Abstract：In order to deeply understand the motion mechanism of bubbles in turbulent flows，instantaneous forces exerting on bubbles in a channel turbulence flow with a low Reynolds number were analyzed in details with the Euler-Lagrange one-way method.The liquid-phase turbulence field was computed with direct numerical simulations（DNS） in Euler reference frame，and instantaneous forces on bubbles were solved with Newtonian motion equations including drag，shear lift，pressure gradient，virtual mass，gravity（i.e.，buoyant） and inertia forces in Lagrange reference frame.The present investigations showed that forces acting on bubbles depended simultaneously on the gravity direction，the flow direction of liquid and the bubble place in the wall-normal direction.The instantaneous force components controlling bubble motions changed with bubbles location and motion direction.Comparatively speaking，the component of the shear lift force perpendicular to gravity had an important influence on the bubble motion in the region near walls.The wall-normal component of the pressure gradient force was dominant for bubble motions in the region away from walls.The drag force components was very important for bubble motions in the whole channel region，and the virtual mass force played a dominant role in bubble motions in all regions except one near the walls for the vertical channel.

Key words：Euler-Lagrange model；one-way；direct numerical simulations；bubbles；force analysis