宋 濤，蔣開喜，周俊武,3，馮昱清

（1.東北大學 材料與冶金學院，遼寧沈陽110819；2.北京礦冶研究總院，北京市100160；3.礦冶過程自動控制技術(shù)北京市重點實驗室，北京市102628；4.CSIRO Mineral Resource Flagship,Clayton VIC 3168,Australia）

機械攪拌反應器內(nèi)氣體分布的數(shù)值模擬

宋 濤1,2,3，蔣開喜2，周俊武2,3，馮昱清4

（1.東北大學 材料與冶金學院，遼寧沈陽110819；2.北京礦冶研究總院，北京市100160；3.礦冶過程自動控制技術(shù)北京市重點實驗室，北京市102628；4.CSIRO Mineral Resource Flagship,Clayton VIC 3168,Australia）

機械攪拌反應器內(nèi)通常包含氣、液、固等多相體系，利用計算流體力學方法針對其內(nèi)部復雜物料運動的模擬近年來取得了很大進展。歐拉雙流體模型被用來模擬一種裝配Rushton型葉輪機械攪拌反應器內(nèi)包含氣泡的氣液兩相流，并通過不同的相間力模型及自定義湍流子模型描述槽內(nèi)氣泡復雜運動。分析了不同曳力、氣泡誘導湍流、升力、湍流分散力等模型在模擬氣泡分布時的效果，并將模擬結(jié)果與Barigou和Greaves的試驗結(jié)果進行比較。結(jié)果表明，提出的模擬方法在機械攪拌反應器模擬中能夠得到準確的氣體分布結(jié)果。

機械攪拌反應器；氣泡流；氣液兩相流；數(shù)值模擬

機械攪拌反應器廣泛應用于礦物加工、濕法冶金、化工工業(yè)、食品工業(yè)、廢水處理等過程工業(yè)。反應器內(nèi)通常包含氣、液、固等多相體系，各相之間的混合交互也在一定程度上影響著槽內(nèi)流場。近年來使用計算流體力學方法對機械攪拌反應器內(nèi)多相流場的研究發(fā)展很快，該方法在投入較少的情況下能獲得較完備的流場信息，成為進行設(shè)備設(shè)計和優(yōu)化的有力工具。

許多學者使用歐拉—歐拉方法進行機械攪拌反應器內(nèi)氣液兩相流的模擬，并取得了不錯的結(jié)果[1-3]。通常情況下,曳力對氣泡在反應器內(nèi)的運動影響很大[4]。常用的幾種曳力模型基本上來源于單個氣泡在靜止液體中的運動,對于類似于機械攪拌反應器內(nèi)的湍流環(huán)境能否適用還需要探討[5]。Ishii和Zuber[6]提出的曳力系數(shù)模型,主要考慮不同氣泡直徑下氣泡運動形狀和氣泡簇效應對氣泡終端速度的影響,但沒有考慮湍流環(huán)境對氣泡運動的影響。Pinelli等[7]通過對小顆粒在湍流環(huán)境中的運動總結(jié)得到湍流環(huán)境中顆粒終端速度變化規(guī)律,得到相關(guān)曳力模型,并且獨立于葉輪樣式,反應器尺寸和液體黏度。Fajner等[8]則是在Pinelli等的基礎(chǔ)上考慮了相間密度差的影響,并重新修正了式中參數(shù)。盡管Pinelli模型和Fajner模型一般用于顆粒在湍流環(huán)境中運動的計算,但在機械攪拌反應器中的湍流環(huán)境中,尺寸小的氣泡與顆粒運動狀態(tài)類似，可以考慮采用此種模型進行氣泡運動的模擬。另外,在存在大量氣泡的情況下,氣泡對液相流場的影響也不應該忽視。氣泡運動帶來的液相湍流脈動可稱作氣泡誘導湍流,一些氣泡誘導湍流模型在機械攪拌反應器內(nèi)氣液流場模擬中的表現(xiàn)值得關(guān)注。其他氣泡與液相之間的相互作用力,例如升力、湍流分散力對氣體分布的模擬結(jié)果影響也很大，如何在機械攪拌反應器內(nèi)混合環(huán)境中使用也需要進行試驗。

本文使用歐拉雙流體模型模擬了一種機械攪拌反應器內(nèi)的氣液流場和氣體分布。對比了3種不同曳力系數(shù)模型的計算效果,找到適合攪拌環(huán)境中氣泡運動的曳力系數(shù)模型;并通過數(shù)值試驗的方法對比了幾種氣泡誘導湍流、升力和湍流分散力對氣體分布模擬結(jié)果的影響。通過模擬結(jié)果與Barigou和Greaves[9]的試驗值的比較,確定了較為準確的機械攪拌反應器內(nèi)氣液兩相流模擬模型。

1 模型描述

1.1 控制方程

對于機械攪拌反應器槽內(nèi)氣液兩相流系統(tǒng)，采用歐拉雙流體模型進行描述。在不考慮氣泡與液相之間的質(zhì)量與熱量傳遞時，氣相和液相運動分別由以下時均質(zhì)量和動量守恒方程控制[10]（α=l時為液相控制方程，α=g時為氣相控制方程）。

式中：rα為各相的體積分數(shù)，ρα為各相密度，t為時間，Uα為各相平均速度矢量，ρα為各相壓力。Fα相間動量傳遞，包括曳力、升力、虛擬質(zhì)量力、壁面潤滑力和湍流分散力等。這里的μeff為等效黏度，可以認為是分子黏度與湍流黏度的和的形式，即μeff=μo+μtα。湍流黏度μtα=C1ραkα2/εα，k，ε分別為湍動能和湍流耗散率，C1為k-ε湍流模型常數(shù)，取0.09。對于氣相采用零方程模型描述湍流，對于液相采用k-ε兩方程模型描述。

1.2 曳力

曳力是由氣泡周圍液體施加在其上，并對加速上升的氣泡產(chǎn)生相反方向阻礙的力。在機械攪拌反應器內(nèi)氣液兩相流模擬中，曳力（FD）可以看作是決定氣體分布的最主要因素。其可以寫作以下形式：

式中：dg為氣泡直徑,CD為曳力系數(shù)。曳力系數(shù)依賴于流體性質(zhì)、氣泡直徑、重力以及氣液相界面間的污染程度[11]。

Ishii和Zuber針對完全雜質(zhì)水環(huán)境下的氣泡、顆粒、液滴的運動規(guī)律提出的曳力系數(shù)模型考慮了氣泡形狀和氣泡間相互作用的因素,使用分段模型的方式描述不同氣泡的終端速度。

式中:UT為靜止液體環(huán)境下氣泡的終端速度。Mendelson[12]提出了純水中1.3 mm至15 mm氣泡終端速度的計算公式：

式中:sphere、ellipse、cap分別代表小直徑、中等直徑和較大直徑不同形狀氣泡的曳力系數(shù)。Rem為修正的氣泡雷諾數(shù)，E(rg)和(1-rd)2為反映多氣泡環(huán)境下氣泡間相互作用對曳力系數(shù)影響函數(shù)。

依據(jù)靜止液體中氣泡的受力平衡，曳力系數(shù)模型可以用以下形式表示：

式中:σ為液體表面張力系數(shù)。

為了考慮湍流環(huán)境對氣泡終端速度的影響，Pinelli等在一種稀釋溶液的湍流環(huán)境中建立了湍流環(huán)境中的顆粒終端速度（US）與靜止液體環(huán)境中顆粒終端速度（UT）的經(jīng)驗關(guān)系式：

式中:λ/dg為Kolmogorov微尺度長度λ=(v3/ε)0.25（v為液相運動黏度）與氣泡直徑的比值。使用US代替公式（10）中的UT，即可計算湍流環(huán)境下的曳力系數(shù)。

Fajner等在Pinelli等研究的基礎(chǔ)上對模型進行了修改，并在顆粒分布模擬中取得了較好的效果：

雖然Pinelli等和Fajner等提出的曳力計算模型原用來模擬攪拌系統(tǒng)中顆粒運動，許多學者使用類似的模型較好地模擬了在湍流環(huán)境下的氣泡運動分布。所以本文也測試了這些模型對機械攪拌式反應器內(nèi)氣液流場的模擬效果。

1.3 氣泡誘導湍流

Kataoka和Serizawa[13]認為在包含氣泡的氣液兩相流模擬中，液相流體會受到氣泡運動的影響，如果僅僅使用通過單相液體運動特性得到的湍流模型，不能夠準確預測流場內(nèi)液相湍流信息。需要考慮在單相流湍流模型的基礎(chǔ)上加入氣泡誘導湍流。宋濤等[14]曾在氣升式環(huán)流反應器氣液兩相流模擬中使用向k和ε的守恒方程中添加源項的方法，來描述氣泡誘導湍流所引起的湍動能的對流和擴散現(xiàn)象，并取得了較好的模擬效果。本文也嘗試將此模型用于機械攪拌反應器內(nèi)氣液兩相流的模擬，分別在k-ε湍流方程中添加以下源項：

式中:Ck和Cε為可調(diào)整的系數(shù)，這里都取1。

Sato和Sekoguchi[15]提出的模型，通過將湍流應力張量中的等效黏度分為剪切誘導項和氣泡誘導項的方法進行模擬：

式中:μ1,SI=μtl為剪切誘導黏度，μ1，BI為氣泡誘導黏度。本文也將考察此模型在機械攪拌反應器內(nèi)氣液兩相流模擬中的效果。

1.4 升力

升力描述的是氣泡與液相剪切場之間的相互作用，是氣泡在上升過程中水平方向運動的決定因素之一。在氣液兩相流模擬中，升力模型的選取往往會對計算結(jié)果造成較大影響。單位體積升力（FL）的計算公式為：

式中:CL為升力系數(shù)。這里考慮考察兩種升力模型。式(19)為使用固定的升力系數(shù)；式(20)則按照Tomiyama[16]提出的基于單氣泡的升力系數(shù)模型進行計算。

1.5 湍流分散力

湍流分散力是由湍流渦動和相間曳力共同作用的結(jié)果。有些氣泡在上升過程中會被湍流渦捕獲,并由相間曳力從高濃度區(qū)域輸送到低濃度區(qū)域。

Lopez模型是最早提出的湍流分散力模型之一，其表達式如下[17]：

湍流分散力系數(shù)CTD取值在0.1到0.5之間的模型應用于小直徑氣泡運動模擬時效果較好。

Burns等[18]給出了在歐拉兩相流體系下的湍流分散力（FTD）的一般形式（Favre Averaged Drag Model）：

式中:湍流分散力系數(shù)CTD取值為1；CD為相間曳力系數(shù)；νt1為液相運動黏度；σt1為湍流施密特數(shù)，此處取值0.9。

1.6 模型實現(xiàn)

本文反應器幾何模型使用與Barigou和Greaves相同的機械攪拌反應器，兩相介質(zhì)分別為空氣和水。該反應器為平底反應器，高度和直徑都為1 m。在桶壁四周均勻分布4個擋板，擋板寬度為0.1 m。反應器攪拌系統(tǒng)為直徑0.33 m的標準Rushton葉輪，葉輪距底面0.25 m。充氣裝置位于葉輪正下方槽底中心處，為直徑19 mm的圓形充氣孔。操作條件分別為充氣量Qg=0.001 64 m3/s，葉輪轉(zhuǎn)速180 rpm和充氣量Qg=0.006 87 m3/s，葉輪轉(zhuǎn)速285 rpm。為了節(jié)省計算時間，這里假設(shè)槽內(nèi)流場對稱，建立了完整攪拌槽的1/2區(qū)域，兩個擋板和三片葉輪以反應器的中心截面為對稱。氣泡直徑使用試驗中觀察的平均值2 mm。模擬試驗采用的網(wǎng)格數(shù)為223 393，并利用網(wǎng)格數(shù)為343 533的模型做網(wǎng)格無關(guān)性試驗。具體網(wǎng)格結(jié)構(gòu)如圖1所示。使用多重參考坐標系法模擬葉輪的旋轉(zhuǎn)，即葉輪區(qū)域在旋轉(zhuǎn)坐標系下計算，外層帶有擋板的區(qū)域使用靜止坐標系?？刂品匠淌褂肁NSYS CFX 14.5進行求解，采用自定義子模型的方式對相間力和氣泡誘導湍流進行描述。中心對稱面設(shè)置為旋轉(zhuǎn)循環(huán)界面。模型出口使用一種對于氣相為壓力出口，對于液相為壁面的邊界條件。

圖1 幾何模型及網(wǎng)格劃分

2 結(jié)果與討論

模擬計算分別在不同的操作條件下進行,以涵蓋不同充氣量及葉輪轉(zhuǎn)速條件下裝備Rushton葉輪的機械攪拌反應器內(nèi)的氣液流場。計算結(jié)果與Barigou和Greaves試驗中得到的氣體體積分數(shù)值進行比較,從而考察不同子模型的計算效果。

曳力系數(shù)模型決定氣泡在液體中所受阻力的大小,所以也決定著氣泡在液相中的終端速度。合適的曳力系數(shù)模型能夠改進模擬中氣泡在反應器內(nèi)的分布。在進行模擬實驗之前，首先利用不同的曳力系數(shù)模型計算重力環(huán)境下單個氣泡在液體中的氣泡終端速度,以觀察不同曳力模型所計算的氣泡速度的區(qū)別。圖2是在分別使用Ishii&Zuber模型、Fajner模型和Pinelli模型計算的不同直徑氣泡的終端速度。

圖2 水中不同曳力系數(shù)模型下氣泡終端速度與氣泡直徑之間的關(guān)系

圖2中可以看出2 mm直徑氣泡在使用Ishii& Zuber模型計算曳力系數(shù)時所獲得的終端速度最大, Fajner模型得到的2 mm氣泡終端速度要小很多,而Pinelli模型計算得到的2 mm氣泡終端速度最小。顯然依據(jù)Fajner模型和Pinelli模型,湍流環(huán)境下單氣泡終端速度要小于靜止流體中單氣泡終端速度。

表1是在使用三種不同曳力系數(shù)模型時槽內(nèi)充氣量模擬值與試驗值的比較。從對比結(jié)果可以看出,使用Ishii&Zuber模型時得到的槽內(nèi)氣體體積分數(shù)最小,同時在兩種操作條件下都與試驗值有較大差別，嚴重低估了槽內(nèi)氣體體積分數(shù)。使用Fajner模型時，氣體體積分數(shù)預測值有了較大改善，接近試驗值，但在兩種操作條件下都較試驗值小一些。使用Pinelli模型時，在較高充氣量條件下，氣體體積分數(shù)模型預測值要比試驗值高約9.3%，而在較低充氣量條件下預測值在幾種模型中最接近試驗值，但依然小于試驗值。

表1 不同曳力模型預測值與Barigou和Greaves試驗值[9]的比較

為了進一步改進模型對氣體分布的模擬效果，這里還嘗試使用了不同氣泡誘導湍流模型，升力模型和湍流分散力模型。表2為使用不同模型時氣體體積分數(shù)預測值與Barigou和Greaves試驗值的比較，所有的模型均使用Fajner模型計算曳力系數(shù)。

表2 不同模擬模型預測值與Barigou和Greaves試驗值[9]的比較

表2模型1使用固定升力系數(shù)0.5和系數(shù)為0.1的Lopez湍流分散力模型。Zadghaffari和Moghaddas[19]曾使用此模型計算Rushton葉輪攪拌槽中的氣體分布并取得較好結(jié)果。對比結(jié)果顯示,在兩種充氣量條件下氣體體積分數(shù)預測值都較試驗值稍小。模型2在模型1的基礎(chǔ)上增加宋濤等使用的氣泡誘導湍流模型，得到的氣體體積分數(shù)預測值較模型1減小，此模型在機械攪拌反應器氣液兩相流模擬中效果并不理想。模型3考慮使用Sato和Sekoguchi提出的氣泡誘導湍流模型，計算結(jié)果較模型1有所改進，預測值更加接近試驗值。模型4在模型3的基礎(chǔ)上增加Tomiyama升力模型和Burns等提出的湍流分散力模型，計算結(jié)果顯示使用Burns湍流分散力模型大大降低了槽內(nèi)氣體體積分數(shù)，效果并不理想。模型5使用Sato和Sekoguchi提出的氣泡誘導湍流模型，Tomiyama升力模型，Lopez湍流分散力模型（系數(shù)調(diào)整為0.25），由計算結(jié)果看出預測槽內(nèi)氣體體積分數(shù)值在兩種操作條件下都與試驗值很接近，尤其是在較高充氣量條件下。

使用模型5時,充氣量Qg=0.006 87 m3/s,葉輪轉(zhuǎn)速285 rpm條件下,模型預測的液相流速矢量分布及氣體體積分數(shù)分別如圖3和圖4所示。由圖3的模擬結(jié)果可以看出，液相流場的高速區(qū)分別位于充氣口上方以及葉輪四周向壁面輻射的區(qū)域,并形成以葉輪平面為分界的上下兩個循環(huán)流場。充氣量Qg=0.001 64 m3/s,葉輪轉(zhuǎn)速180 rpm時,槽內(nèi)液相流場分布與圖3類似。槽內(nèi)中心截面處氣體體積分數(shù)分布如圖4所示。由模擬結(jié)果看出,充氣口進入攪拌槽的氣體被葉輪分散到槽內(nèi)四周,在槽內(nèi)流場形成的兩個循環(huán)區(qū)域中心處,氣體體積分數(shù)值較高,終端速度較小的氣泡會聚集在這些區(qū)域。

圖3 槽內(nèi)中心截面處液相速度矢量分布(充氣量Qg=0.006 87 m3/s，葉輪轉(zhuǎn)速285 rpm)

圖4 槽內(nèi)中心截面處氣體體積分數(shù)分布(充氣量Qg=0.006 87 m3/s，葉輪轉(zhuǎn)速285 rpm)

通過比較不同曳力模型的模擬效果，Ishii& Zuber模型在使用實驗觀察平均氣泡直徑時并不能準確模擬機械攪拌環(huán)境下氣泡在槽內(nèi)的分布，而通過固體小顆粒試驗得到的曳力系數(shù)模型則能較為準確地預測氣體體積分數(shù)值。而合適的氣泡誘導湍流模型，升力模型以及湍流分散力模型能夠很好地補充湍流環(huán)境，多氣泡間相互作用對氣泡運動帶來的影響，修正氣泡分布的預測值，取得較好的模擬效果。

3 結(jié)論

本文使用歐拉雙流體模型模擬了一種實驗室機械攪拌反應器內(nèi)的氣液流場和氣體分布。通過比較不同的曳力系數(shù)模型，不同的氣泡誘導湍流、升力和湍流分散力模型的模擬效果，考察了這些模型在氣體分布模擬計算中的作用。通過這些模型模擬值與Barigou和Greaves試驗值的對比，提出一種能夠準確預測反應器內(nèi)氣體分布的模擬模型。本文模擬模型可以為機械攪拌反應器內(nèi)氣液兩相流模擬模型的建立提供參考依據(jù)，并在大型工業(yè)攪拌槽、浮選機等攪拌設(shè)備的放大過程中，提供相關(guān)的理論支持。

[1] Wang Weijing,Mao Zaisha.Numerical Simulation of Gas-Liquid Flow in a Stirred Tank with a Rushton Impeller[J].Chinese J. Chem.Eng.,2002,10(4):385-395.

[2] Lane G L,Schwarz M P,Evans G M.Numerical modelling of gasliquid flow in stirred tanks[J].Chemical Engineering Science, 2005,60(8-9):2203-2214.

[3] Scargiali F,D'Orazio A,Grisafi F,Brucato A.Modelling and simulation of gas-liquid hydrodynamics in mechanically stirred tanks[J].Chemical Engineering Research and Design,2007,85(5): 637-646.

[4] Sokolichin A,Eigenberger G,Lapin A.Simulation of buoyancy driven bubbly flow established simplifications and open questions [J].AIChE Journal,2004,50(1):24-45.

[5] Bakker A,Van den Akker H E A.A computational model for the gas-liquid flow in stirred reactors[J].Chemical Engineering Research and Design,1994,72(A4):594-606.

[6] Ishii M,Zuber N.Drag coefficient and relative velocity in bubbly, droplet or particulate flows[J].AIChE Journal,1979,25(5):843-855.

[7] Pinelli D,Nocentini M,Magelli F.Solids distribution in stirred slurry reactors:influence of some mixer configurations and limits to the applicability of a simple model for predictions[J].Chemical Engineering Communications,2001,188(1):91-107.

[8] Fajner D,Pinelli D,Ghadge R S,Montante G,Paglianti A,Magelli F.Solids distribution and rising velocity of buoyant solid particlesin a vessel stirred with multiple impellers[J].Chemical Engineering Science,2008,63(24):5876-5882.

[9] Barigou M,Greaves M.Bubble-size distributions in a mechanically agitated gas-liquid contactor[J].Chemical Engineering Science, 1992,47(8):2009-2025.

[10] Luo H-P,Al-Dahhan M H.Verification and validation of CFD simulation for local flow dynamics in a draft tube airlift bioreactor [J].Chemical Engineering Science,2011,66(5):907-923.

[11] Tomiyama A,Kataoka I,Zun I,et al.Drag coefficients of single bubbles under normal and micro gravity conditions[J].JSME International Journal Series B,1998,41(2):472-479.

[12] Mendelson H D.The prediction of bubble terminal velocities from wave theory[J].AIChE Journal,1967,13(2):250-253.

[13] Kataoka I,Serizawa A.Basic Equations of turbulence in gasliquid two-phase flow [J].International Journal of Multiphase Flow,1989,15(5):843-855.

[14] 宋濤，蔣開喜，周俊武，等.氣升式環(huán)流反應器內(nèi)氣泡分布的數(shù)值模擬[J].礦冶，2015，24(1)：61-66.

[15] Sato Y,Sekoguchi K.Liquid velocity distribution in two-phase bubble flow[J].International Journal of Multiphase Flow,1975,2 (1):79-95.

[16] Tomiyama A.Struggle with computational bubble dynamics[J]. Multiphase Science and Technology,1998,10(4):369-405.

[17] Lopez de Bertodano M A.Two fluid model for two-phase turbulent jets[J].Nuclear Engineering and Design,1998,179(1): 65-74.

[18] Burns A D,Frank T,Hamill I,et al.The favre averaged drag model for turbulence dispersion in Eulerian multi-phase flows[C] //Proceedings of 5th International Conference on Multiphase Flow.Yokohama:ICMF'04,2004:392.

[19] Zadghaffari R,Moghaddas J S.Evaluation of drag force effect on hold-up in a gas-liquid stirred tank reactor[J].Journal of Chemical Engineering of Japan,2010,43(10):833-840.

Numerical Simulation OF Gas DIstribution In a Stirred Tank Reactor

SONG Tao1,2,3,JIANG Kaixi2,ZHOU Junwu2,3,FENG Yuqing4
(1.School of Materials and Metallurgy,Northeastern University,Shenyang 110819,China; 2.Beijing General Research Institute of Mining and Metallurgy,Beijing 100160,China; 3.Beijing Key Laboratory of Automation of Mining and Metallurgy Process,Beijing 102628,China; 4.CSIRO Mineral Resource Flagship,Clayton VIC 3168,Australia)

Computational fluid dynamics(CFD)contributes a lot in calculations of complicated multi-phase flows in mechanically stirred reactors.To predict the performance of a tank reactor stirred by Rushton turbine,the Eulerian-Eulerian two-fluid method was used for modeling gas-liquid flow in tank,and extra user defined subroutines were incorporated to consider the complex physics. Different modeling methods of drag force,bubble induced turbulence,lift and turbulent dispersion forces were evaluated in the numerical simulations of the bubble distribution in the reactor.The simulation results have been compared with the experimental measurements of Barigou and Greaves.Some reasonable agreements have been obtained both in flow structures and gas distribution.

stirred tank reactor;bubbly flow;gas-liquid flow;numerical simulation

TF351.5+2

A

1004-4345(2015)03-0022-05

2015-04-01

十二五“863”項目《大型微生物冶金生物反應器》(2012AA061503)。

宋濤（1983—），男，東北大學博士研究生，主要從事礦冶設(shè)備、流程建模仿真研究。