基于EMMS方法的鼓泡塔反應器CFD及群平衡模擬

王玨，楊寧

（1中國科學院過程工程研究所多相復雜系統(tǒng)國家重點實驗室，北京 100190；2中國科學院大學，北京 100049）

基于EMMS方法的鼓泡塔反應器CFD及群平衡模擬

王玨1,2，楊寧1

（1中國科學院過程工程研究所多相復雜系統(tǒng)國家重點實驗室，北京 100190；2中國科學院大學，北京 100049）

能量最小多尺度 （energy-minimization multi-scale，EMMS） 方法已經被應用于氣液體系中群平衡（population balance model，PBM）模型的改進。EMMS模型可計算氣泡破碎聚并過程的能量，進而獲得聚并速率的修正因子。應用這一模型對高氣速鼓泡塔進行了模擬計算，并進一步對比了均一尺徑模型、CFD-PBM 模型以及CFD-PBM-EMMS模型的模擬結果與實驗數據。結果表明，在高表觀氣速條件下，基于EMMS方法的群平衡模型可以更加準確地預測鼓泡塔中不同高度的氣泡尺徑分布和軸向液速，同時提高了對整體氣含率和局部氣含率的模擬準確性。在表觀氣速為0.16 m·s-1和0.25 m·s-1時，CFD-PBM-EMMS模型對氣泡尺徑分布的預測精度更高，同時整體氣含率模擬的相對誤差下降為5%和15%，局部氣含率模擬平均相對誤差下降為8%和17%。

計算流體力學；群平衡模型；鼓泡塔；氣含率；氣泡尺徑分布

引 言

鼓泡塔反應器具有結構簡單、操作簡便、良好的傳熱傳質效率等優(yōu)點，被廣泛應用于化學工程、生物工程、石油工程等領域[1]。隨著計算機技術的發(fā)展，計算流體力學（computational fluid dynamics，CFD）方法已經成為鼓泡塔反應器模擬和優(yōu)化的重要工具。對于鼓泡塔中的氣液兩相流動常見的模擬方法有歐拉-歐拉方法和歐拉-拉格朗日方法。其中，歐拉-拉格朗日方法需要對每一個氣泡進行跟蹤，對計算資源要求高，適用于氣含率較小的反應器模擬[2]；而歐拉-歐拉方法假定離散相與連續(xù)相之間可相互滲透，不直接追蹤每一個氣泡的運動，從而顯著地降低了計算量，被廣泛地應用于實驗室規(guī)模和工業(yè)規(guī)模的鼓泡塔模擬之中。

歐拉-歐拉雙流體模型中，氣液相間作用力模型對模擬結果至關重要[3]。相間作用力包括曳力、升力、虛擬質量力、湍流擴散力等，其中曳力對氣液系統(tǒng)的影響占主導地位[4-7]。在曳力模型中，氣泡尺寸常被作為一個重要的輸入參數，直接影響模擬結果。模型設置中，通常依據實驗或經驗將氣泡直徑設定為定值。然而實際系統(tǒng)中，由于氣泡之間存在破碎和聚并的現象，導致氣泡直徑在反應器內呈現動態(tài)時空非均勻分布特征。針對此特征，群平衡模型（population balance model，PBM）提供了一個解決方案，對氣泡數密度建立傳輸方程，考慮氣泡聚并和破碎對氣泡數密度變化的影響。因此，群平衡模型耦合計算流體力學方法（CFD-PBM）對鼓泡塔進行模擬已經成為近年來的熱點之一[8]。

群平衡模型中，聚并核函數和破碎核函數直接決定了氣泡粒徑分布規(guī)律的模擬結果[9]。聚并過程通常簡化為二元聚并的情況[10]，并基于物理現象推導和實驗測量提出了不同的聚并核函數[11-12]。Chen等[13]發(fā)現現有的破碎聚并函數并不能使氣泡群的聚并速率和破碎速率在鼓泡塔穩(wěn)定鼓泡區(qū)達到動態(tài)平衡，這與實驗現象不符。Wang等[14]提出，當氣泡之間的距離比氣泡湍動路徑長時，氣泡的實際碰撞頻率要比Prince模型小，導致模型高估了氣泡間的聚并效率。Bhole等[15]認為現有聚并模型在計算碰撞頻率時使用對應尺度渦的速度代替氣泡速度，而氣泡與湍動渦之間是存在速度差的，所以聚并速率應該乘以一個小于1的系數。其他學者也提出，需要使用小于1的修正因子來修正聚并函數以減小聚并速率[16-18]。肖頎[19]和 Yang等[20]率先提出了基于 EMMS方法的群平衡（PBM）模型模擬鼓泡塔的CFD-PBM-EMMS模型，模擬結果表明了該模型可以合理地模擬實驗室規(guī)模中低氣速鼓泡塔中的氣泡粒徑分布。Qin等[21]采用了EMMS-PBM的方法對液-液體系的轉子定子中的液滴尺徑分布進行了修正，提高了原有模型的準確度。本文針對高表觀氣速的氣液鼓泡塔，應用CFD-PBM-EMMS方法對高氣速氣液鼓泡塔進行了模擬，通過與實驗數據對比驗證了該方法的可行性。

1 數學模型和數值模擬

本文中的三維CFD模型以歐拉-歐拉模型為基礎，采用重整化群（RNG）k-ε湍流模型和Schiller-Naumann曳力模型。同時，耦合群平衡模型（PBM）計算氣泡尺徑分布，其中描述氣泡聚并和破碎過程分別采用 Luo聚并模型和 Luo &Svendsen破碎模型，并使用離散法（class method，CM）求解群平衡方程。

1.1 控制方程

氣液系統(tǒng)模擬采用歐拉-歐拉模型，歐拉-歐拉模型又稱雙流體模型（two-fluid model，TFM）。模擬中把液體設置為連續(xù)相，氣體設置為離散相，并且不考慮氣體的壓縮和相間質量、熱量的傳遞。控制方程可以表示為以下形式。

質量守恒方程

動量守恒方程

式中，α、ρ、u、t分別表示相含率、密度、速度矢量和時間；P表示壓力；μeff表示流體有效黏度；Mi指相間動量交換項。

式中，各項對應的有效黏度μeff是液相分子黏度μL和湍流黏度μt之和。其中湍流黏度為

式中，Cμ=0.0845，湍動能 k和湍動能耗散率ε可從湍流模型耦合計算得到。

1.2 湍流模型

湍流模型按照文獻[22-23]中的選擇，采用重整化群（RNG）k-ε混合相模型（Mixture），模擬湍動狀態(tài)下的氣液系統(tǒng)[13]。

式中，C1ε=1.42，C2ε=1.68。

在此 Mixture模型中，氣液兩相共用一套方程，方程中采用混合密度ρm和速度um，計算公式如下

湍流動能產生項Gk,m用以下公式計算得到

式中，η≡Sk/ε(S 是張量的積)，η0=4.38，β=0.012。

1.3 相間作用力

根據Sanyal等[4]和Van Baten等[5]的研究結果，本文的模擬計算只考慮曳力，忽略其他相間作用力，曳力系數選擇 Schiller-Naumann[24]模型模擬氣液體系相間作用力[11]。

其中，曳力系數為

式中，CD是氣泡群曳力系數，CD,0為單氣泡曳力系數，CD=CD,0(1-αg)p，本模擬中p=1。

1.4 群平衡模型

群平衡模型對氣泡在時間和空間中進行的變化進行追蹤，可以很好地預測氣泡尺徑的分布。在群平衡模型中氣泡尺徑的變化與聚并速率和破碎速率密切相關。而沸騰、冷凝還有氣泡在邊界上的生成和減小在本文中暫不考慮。

式中，n是氣泡的數密度；c，b函數分別為聚并核函數和破碎核函數；β是母氣泡破碎時子氣泡的分布函數；式（13）右端的第1項和第2項分別表示由聚并引起的氣泡產生和消亡；式（13）右端的第3項和第4項則分別表示由破碎引起的氣泡產生和消亡。

1.4.1 破碎模型 本文選用Luo & Svendsen[25]的破碎模型

其中

1.4.2 聚并模型 Luo[26]的聚并模型為

式中，聚并速率c等于碰撞頻率?c與聚并效率Pc的積。

其中

1.5 聚并速率修正因子

能量最小多尺度（EMMS）模型由Li[27]提出，現已廣泛應用于氣固流化床反應器中。Yang等[28]基于EMMS模型提出了適用于氣-液體系的DBS模型，并且應用于鼓泡塔、攪拌槽等氣液反應器。對于鼓泡塔體系，DBS模型將氣相簡化為“大氣泡相”和“小氣泡相”兩相，通過大小氣泡的直徑(dL和dS)、氣含率(αL和αS)和表觀速度(Ug,L和 Ug,S)來描述氣液系統(tǒng)的狀態(tài)。DBS模型的主要模型方程[28]如下。

質量守恒方程

力平衡方程

對氣液體系的能量關系的描述如下。

總能量

介尺度能耗

微尺度能耗

式中，CD,p為固體顆粒曳力系數；CD,b為氣泡曳力系數；?為渦與氣泡碰撞頻率；αb為氣含率；Pb為氣泡破碎概率；fv為破碎前后子氣泡與母氣泡體積比；cf為氣泡破碎時表面積增加量，等于σ為表面張力。

在鼓泡塔的氣液體系中，存在氣相控制和液相控制兩種機制。當氣相控制為主導時，液相的湍動能趨于最小（Nturb→ min），而液相控制時，氣泡在湍流中表面振蕩能耗最小（Nsurf→ min）。在這兩種控制機制都發(fā)揮作用的復雜流型之下由兩種機制協(xié)調控制，DBS模型穩(wěn)定條件為

DBS模型在當前的氣液體系研究中的新進展見文獻[29-30]。Yang等[31]在DBS的基礎上提出了氣液體系DBS曳力模型，在CFD模擬鼓泡塔的氣含率方面優(yōu)勢明顯，并推廣至攪拌槽中的氣含率預測[32]。最近，肖頎[19]和Yang等[20]將氣液體系內介尺度能量的描述和 PBM 模型相結合，提出了對聚并模型的修正方法，如圖1所示。氣液體系在穩(wěn)定狀態(tài)下，微尺度能量趨于最小，同時不同尺徑之間的氣泡破碎和聚并處于動態(tài)平衡。此時在聚并相前加入抑制聚并的修正因子 C，在使用 PBM 求解的介尺度能耗 Nbreak_PBM和 DBS求解的介尺度能耗Nbreak_DBS相等時，就得到由DBS模型進行修正的PBM修正因子。

破碎核函數：Luo & Svendsen[25]；聚并核函數：Luo 對應修正因子關系式[26]。

圖1 聚并核函數修正因子求解流程[20]Fig. 1 Flow chart for calculation of coalescence corrector[20]

2 建模對象和方法

2.1 結構模型設置和網格無關性驗證

本文用 GAMBIT建立鼓泡塔的三維物理模型并且劃分網格，在進氣口采用面進氣的方式，對分布器上部鼓泡區(qū)域模擬計算。鼓泡塔塔徑390 mm，塔高2000 mm，液面高為1000 mm。為了獲得合適尺寸的網格，本文對粗（coarse）、中（medium）、細（fine）3種不同尺寸的網格尺寸進行了分析。圖2表示了三維網格劃分方式，并給出了中等尺寸網格的俯視圖和側視圖。表1為分布器不同的劃分尺寸。模型分布器與McClure等[33]的保持一致，中間部分為分布器管道壁面，兩邊為進氣區(qū)域。進氣區(qū)域采取面進氣的模型設置，各個區(qū)域網格結構如圖2所示。

圖2 網格結構Fig. 2 Mesh configuration

表1 網格無關性驗證模型網格尺寸及數目Table 1 Mesh size and cells adopted for grid independence test

圖3給出了3種不同尺寸網格模擬得到的不同高度截面的平均氣含率分布。3種不同尺寸的網格均能大致反映出鼓泡塔中氣含率的軸向變化趨勢。粗網格的網格數目為17萬個，因為網格數目偏低，導致計算時產生了一定的誤差。而在網格數目大于36萬個之后，網格數目對模擬結果的影響較小。因此，考慮到計算速度和模擬結果準確性，本文后續(xù)模擬工作均采用中等尺寸的網格。

圖3 不同網格尺寸模擬所得局部氣含率分布Fig.3 Local gas holdup at different column height of different mesh size

2.2 數學模型求解設置：

模擬計算基于Ansys Fluent 15.0軟件，采用非穩(wěn)態(tài)模型方程。鼓泡塔入口處采用gas-velocity邊界條件，出口處使用pressure-outlet邊界條件。壁面使用無滑移邊界。為了保持計算的穩(wěn)定性，計算初始時使用時間步長Δt=0.0001 s，計算1000步，增大時間步長至Δt=0.001 s，計算1000步至1.1 s，最終使用Δt=0.005 s穩(wěn)定計算至120 s。壓力速度耦合方式采用SIMPLE格式離散，動量方程、體積分數以及湍流方程均采用一階迎風格式離散，松弛因子使用0.75的默認設置。

圖4為鼓泡塔中不同高度下的Sauter平均尺徑（d32）和7 mm的氣泡尺徑體積分數隨著時間變化的模擬結果。由圖4（a）可以看出，在40 s之后，不同高度截面的 d32基本相同，認為此時鼓泡塔已經進入較為穩(wěn)定的狀態(tài)。由圖中可以觀察到，氣泡群的平均尺徑沿著塔高逐漸增加。由圖4（b）可以看出，在40 s之后，7 mm氣泡的體積分數呈現出周期性的變化趨勢，流場數據均取40～120 s的數據進行時間平均。

本文中群平衡模型的求解方法采用離散法，離散法具有概念簡單、容易理解的特點，能夠給出直觀的顆粒粒級分布。

離散法處理各個氣泡粒級

圖4 氣泡尺寸隨時間的變化趨勢Fig. 4 Instantaneous variation of bubble size

同時氣泡個數和氣含率之間存在

式中，Ni、αi、Vi分別表示第i個氣泡粒級的氣泡數密度函數、氣含率、體積。

TFM 中用于描述氣泡變化的氣含率αi項和PBM 方程中描述不同氣泡濃度的 Ni的關系即可獲得。結合McClure等[33]實驗數據，將氣泡尺徑范圍設定為1，3，5，……，39 mm 20個粒級，入口的初始顆粒粒徑使用Bhole等[15]使用的5 mm。

3 模擬結果與討論

本文使用3種模型進行模擬，如表2所示。case 1采用均一氣泡尺徑的雙流體（TFM）模型；case 2采用群平衡（PBM）模型耦合雙流體模型，群平衡（PBM）核函數采用 Luo聚并函數和 Luo &Svendsen破碎模型；case 3也采用群平衡（PBM）模型耦合雙流體（TFM）模型，不同的是加入了EMMS方法對聚并模型修正。經計算，表觀氣速為0.16 m·s-1時修正因子為 0.1230，表觀氣速為 0.25 m·s-1時計算出的修正因子為0.1051。3種模型中均使用Schiller-Naumann曳力模型，入口處氣泡的尺徑均為5 mm。

表2 不同模型算例的設置Table 2 Different model settings

3.1 氣含率

圖5為不同模型鼓泡塔整體氣含率模擬和實驗結果的對比。由圖可知高氣速下傳統(tǒng)的雙流體模型明顯高估氣含率，0.16 m·s-1和0.25 m·s-1表觀氣速下相對誤差分別為40%和80%。可能原因在于，高氣速條件下，氣泡受到液相湍動的影響較為明顯，氣泡在不同高度段有著不同的尺徑分布，氣泡的實際尺徑不能使用入口氣泡的尺徑近似處理。Schiller-Naumann曳力模型在計算氣含率方面，也受到尺徑的影響。Xu等[34]發(fā)現，在表觀氣速為0.1 m·s-1時，氣泡尺徑設置為2.5、4.5、6.5 mm均一氣泡尺徑時，采用Schiller-Naumann曳力模型的條件下發(fā)現，氣含率會隨著氣泡尺徑的增加而降低。case 1假定了均一的氣泡尺徑，沒有描述氣泡經過聚并破碎后平均尺徑沿塔高增大的現象，導致模擬結果中的氣含率偏高。相比較而言 PBM 模型描述了氣泡上升過程中的破碎與聚并現象，得到了氣泡的尺徑分布。未經修正的群平衡模型略低估了氣含率，在0.16 m·s-1和0.25 m·s-1的表觀氣速下整體氣含率模擬數據和實驗結果之間的相對誤差為20%和15%，而經過修正之后的群平衡模型因提供了更為準確的氣泡尺徑信息，從而預測的氣含率模擬結果和實驗數據之間的相對誤差為5%和15%。

圖5 整體氣含率模擬數據和實驗結果[33]的對比Fig. 5 Comparison between different models and experiments[33]

圖6 不同模型整體氣含率云圖Fig. 6 Contour plots for different models

圖6為不同模型的時均氣含率分布云圖。所有case均能反映出鼓泡塔中氣含率分布的結構特征。進氣口處受分布器受非進氣區(qū)域影響，塔中心氣含率偏低。在鼓泡塔中上高度處，case 1的氣含率分布較為均勻，這與使用均一氣泡尺徑相關。在Usup=0.25 m·s-1的情況下，較高的氣含率導致模擬液面高于鼓泡塔頂部。同時在case 2、case 3中，可以明顯地看到鼓泡塔中的氣含率在軸向和徑向不同位置的不均勻分布；其中case 2由于氣含率偏低，導致鼓泡液位高度與未鼓泡前的靜液位高度相比差距并不明顯。

圖7 局部氣含率模擬數據和實驗結果[33]的對比Fig.7 Comparison of local gas holdup between different models and experiments[33]

圖7是不同高度處的局部氣含率沿徑向的分布。McClure等[33]選取了鼓泡塔中段和接近相界面處的兩個截面記錄了實驗數據。case 1、case 2和case 3的模擬數據同實驗結果相比，氣含率的徑向分布結果均呈現出中間高，邊壁低的分布規(guī)律。分布器非均勻分布結構使得軸心處的氣含率略有偏低，3個case也都捕捉到了這一現象。與整體氣含率相同的是，局部氣含率的結果也表現出了case 1明顯高估了氣含率，case 2、case 3在誤差允許的范圍內都可以較為可靠地反映鼓泡塔z=550 mm處的局部氣含率徑向分布。case 2的整體氣含率低于case 3。

3.2 氣泡大小分布

3.2.1 修正因子對氣泡尺徑分布的影響 圖8給出了 case 2、case 3 在 0.16 m·s-1和 0.25 m·s-1表觀氣速下鼓泡塔z=550 mm和z=1050 mm處氣泡尺徑分布的模擬值和實驗值的對比，其中case 2中的氣泡峰值位于39 mm的粒級。小于39 mm的其他粒級的氣泡，分布較為平緩。分析原因：① case 2 中的氣泡粒級受到1～39 mm粒級區(qū)間設置的限制和聚并核函數對聚并速率的高估，使粒徑的峰值為 39 mm，而且從37 mm到39 mm的粒級間陡然上升，由于原始聚并模型高估了聚并速率[11]，模擬得到的粒徑分布顯示出大氣泡較多的趨勢；② 從曲線整體分析，1～37 mm的氣泡尺徑分布較為平緩，與實驗結果相比，說明了模擬中氣泡的聚并被高估了。case 3嘗試使用修正因子抑制聚并函數，修正氣泡聚并速率。修正后的模型中，不同尺徑氣泡的體積分布范圍集中于5～15 mm，9 mm的氣泡具有最高的體積分數。兩種不同的模型對比，修正后的群平衡模型預測的粒徑分布更加符合實際。實驗中體積分數最高的氣泡粒級為 5 mm，與模型預測結果仍有一定差距。McClure等[33]分析實驗結果認為，在H/d=2、H/d=3截面上，氣泡尺徑分布在實驗允許的誤差范圍內是相同的，而帶有修正因子的模擬結果也反映了這一點。因此，在H/d≥2截面上，氣泡的聚并和破碎達到了動態(tài)平衡。

圖8 氣泡尺徑分布模擬數據和實驗結果[33]的對比Fig.8 Comparison of bubble size distribution between different models and experiments[33]

3.2.2 Sauter平均直徑 除模擬出高氣速下氣泡尺徑在不同位置處的動態(tài)變化外，群平衡模型可進一步獲得氣泡的Sauter平均尺徑（d32）信息，此平均直徑將用于雙流體模型中的耦合計算。圖9給出了case 1中的雙流體模型采用均一氣泡的設定，氣泡直徑保持了5 mm的初始值不變。未經修正的PBM和 CFD耦合模型計算出的氣泡在不同徑向位置上的d32在8～18 mm范圍變化，而修正后的PBM模型由于聚并受到了抑制，平均尺徑減小，H/d=2到H/d=3徑向位置上平均粒徑在7～8 mm范圍變化。

3.3 軸向液速

圖10 中可以看出，不同氣泡尺徑的條件下，3個不同case模擬的軸向液速相差不大，說明氣泡尺徑對軸向液速的模擬的影響不大。分布器結構導致鼓泡塔軸心處液速較為平緩，與塔中心處的氣含率較低（圖7）情況相同。在z=1050 mm處，模擬的軸向液速都較實際值偏低，與徐琰等[35]對于離子液體的氣液兩相流的模擬類似，其模擬結果也出現接近塔頂處液速減小的趨勢。z=550 mm處，不同模型之間差別很小，對于液速的預測大體能符合實驗中獲得的真實流速。

4 結 論

準確地模擬氣泡尺徑分布，描述鼓泡塔中不同空間位置的氣泡大小，對高氣速下氣液兩相流體系有著更為重要的意義。本文采用CFD-PBM模型對高氣速下的氣液鼓泡塔體系進行了模擬，應用EMMS/DBS模型計算出不同氣速下的聚并速率修正因子，研究了不同群平衡模型下的氣液兩相流的流場信息，獲得以下結論。

（1）高氣速條件下的鼓泡塔中發(fā)生強烈的湍動，氣泡的破碎和聚并導致入口處氣泡尺徑和塔中穩(wěn)定鼓泡段的氣泡尺徑有著很大的區(qū)別，此時采用均一氣泡尺徑模型不能準確模擬塔內的所有流場信息。基于EMMS方法修正群平衡模型，考慮了氣泡的聚并和破碎，獲得了具有穩(wěn)定尺徑分布的氣泡群信息，捕捉到氣泡群尺徑由入口處的均一粒徑(5 mm)上升至平均粒徑7 mm的動態(tài)變化。

圖9 Sauter平均直徑在不同表觀氣速下的徑向分布Fig.9 Radial profiles of Sauter mean diameter at different superficial velocities

（2）傳統(tǒng)的 PBM聚并模型需要經過修正才能準確地描述反應器中的實際氣泡分布。一般情況下傳統(tǒng)的聚并核函數模型均高估了聚并速率。使用EMMS-PBM 耦合方法計算出的修正因子有效地抑制了聚并，得到單峰分布的氣泡尺徑曲線。與未修正前的曲線相比，模擬預測更加接近實驗結果。

（3）CFD-PBM-EMMS群平衡模型可以模擬氣泡粒徑從入口到出口處不同高度上的平均直徑的變化，并且在穩(wěn)定鼓泡區(qū)給出了較為準確的氣泡平均尺徑。在表觀氣速0.16 m·s-1和0.25 m·s-1時，使用CFD-PBM-EMMS模型模擬結果中，整體氣含率的相對誤差為 5%和 15%，局部氣含率的平均相對誤差為8%和17%；未經修正的CFD-PBM模型，整體氣含率的相對誤差為20%和15%，局部氣含率的相對誤差為23%和20%。CFD-PBM-EMMS方法優(yōu)化了對高氣速下氣液流場的模擬。

（4）模擬結果表明，EMMS-PBM 耦合方法作為一種改進群平衡方程的思路是可靠的，此方法仍有進一步的改進空間。

符 號 說 明

b——破碎核函數，s-1

C——聚并修正因子

c——聚并核函數，m3·s-1

cf——破碎前后氣泡表面積增加率

D——氣泡直徑，m

dL, dS——分別為大、小氣泡直徑，m

d32——Sauter平均直徑，m

FD——氣泡群曳力，N·m-3

k ——湍動能，m2·s-2

M ——相間動量交換項，kg·m-2·s-2

N——數密度函數

Nbreak——破碎聚并能耗，m2·s-3

圖10 軸向液速模擬結果[33]與實驗數據的對比Fig.10 Comparison of upward liquid velocity between different models and experiment[33]

Nsurf——氣泡表面能，m2·s-3

NT——懸浮輸送能，W·kg-1

Nturb——液相黏性耗散能，m2·s-3

P ——壓強，Pa

Re ——Reynolds數

Ug,L, Ug,S——分別為大、小氣泡表觀氣速，m·s-1

Usup——表觀氣速，m·s-1

V ——氣泡體積，m3

We ——Weber數

α,αg——分別為相含率、氣含率

αL,αS——分別為大、小氣泡氣含率

ε ——湍動能耗散率，m2·s-3

λ ——湍流渦尺度，m

μ ——動力學黏度，Pa·s

ξ ——渦尺寸和氣泡直徑比值

ξij——兩聚并氣泡直徑之比

ρ ——密度，kg·m-3

σ——表面張力，N·m-1

ω——渦與氣泡碰撞頻率，s-1

ωc——兩氣泡碰撞頻率，m3·s-1

[1] KANTARCI N, BORAK F, ULGEN K O. Bubble column reactors[J].Process Biochem., 2005, 40(7)： 2263-2283.

[2] SIMCIK M, MOTA A, RUZICKA M C, et al. CFD simulation and experimental measurement of gas holdup and liquid interstitial velocity in internal loop airlift reactor[J]. Chem. Eng. Sci., 2011,66(14)： 3268-3279.

[3] HIBIKI T, ISHII M. Lift force in bubbly flow systems[J]. Chem. Eng.Sci., 2007, 62(22)： 6457-6474.

[4] SANYAL J, VASQUEZ S, ROY S, et al. Numerical simulation of gas-liquid dynamics in cylindrical bubble column reactors[J]. Chem.Eng. Sci., 1999, 54(21)： 5071-5083.

[5] VAN BATEN J M, KRISHNA R. CFD simulations of a bubble column operating in the homogeneous and heterogeneous flow regimes[J]. Chem. Eng. Technol., 2002, 25(11)： 1081-1086.

[6] KRISHNAR, 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.

[7] SCARGIALI F, D’ORAZIO A, GRISAFI F, et al. Modelling and simulation of gas-liquid hydrodynamics in mechanically stirred tanks[J]. Chem. Eng. Res. Des., 2007, 85(A5)： 637-646.

[8] WANG T F, WANG J F, JIN Y. A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow[J]. Chem. Eng. Sci.,2003, 58(20)： 4629-4637.

[9] 覃成鵬, 楊寧. 多相分散體系中氣泡/液滴聚并和破碎的群平衡模擬[J]. 化學進展, 2016, (8)： 1207-1223.QIN C P, YANG N. Population balance modeling of breakage and coalescence of dispersed bubbles of droplets in multiphase systems[J].Progress in Chemistry, 2016, (8)： 1207-1233.

[10] RAMKRISHNA D. Population Balances： Theory and Applications to Particulate Systems in Engineering[M]. San Diego： Academic Press,2000.

[11] CHEN P, SANYAL J, DUDUKOVI? M P. Numerical simulation of bubble columns flows： effect of different breakup and coalescence closures[J]. Chem. Eng. Sci., 2005, 60(4)： 1085-1101.

[12] RIBEIRO C P, MEWES D. On the effect of liquid temperature upon bubble coalescence[J]. Chem. Eng. Sci., 2006, 61(17)： 5704-5716.

[13] CHEN P, DUDUKOVIC M P, SANYAL J. Three-dimensional simulation of bubble column flows with bubble coalescence and breakup[J]. AIChE J., 2005, 51(3)： 696-712.

[14] WANG T F, WANG J F, JIN Y. Theoretical prediction of flow regime transition in bubble columns by the population balance model[J].Chem. Eng. Sci., 2005, 60(22)： 6199-6209.

[15] BHOLE M R, JOSHI J B, RAMKRISHNA D. CFD simulation of bubble columns incorporating population balance modeling[J]. Chem.Eng. Sci., 2008, 63(8)： 2267-2282.

[16] YAO W, MOREL C. Volumetric interfacial area prediction in upward bubbly two-phase flow[J]. Int. J. Heat Mass Trans., 2004, 47(2)：307-328.

[17] MITRE J F, TAKAHASHI R S M, RIBEIRO C P, et al. Analysis of breakage and coalescence models for bubble columns[J]. Chem. Eng.Sci., 2010, 65(23)： 6089-6100.

[18] MUKIN R V. Modeling of bubble coalescence and break-up in turbulent bubbly flow[J]. Int. J. Multiphas Flow, 2014, 62： 52-66.

[19] 肖頎. 攪拌釜及鼓泡塔內氣液兩相流的多尺度模型[D]. 西安： 西安交通大學, 2015.XIAO Q. Multi-scale modeling of the gas-liquid flow in mixing tanks and bubble columns[D]. Xi’an： Xi’an Jiaotong University, 2015.

[20] YANG N, XIAO Q. A mesoscale approach for population balance modeling of bubble size distribution in bubble column reactors[J].Chem. Eng. Sci., 2017, DOI： 10.1016/j.ces.2017. 01.026.

[21] QIN C, CHEN C, XIAO Q, et al. CFD-PBM simulation of droplets size distribution in rotor-stator mixing devices[J]. Chem. Eng. Sci.,2016, 155： 16-26.

[22] 許婷婷. DBS曳力模型與湍流模型對氣液CFD模擬的影響[D]. 北京： 中國科學院大學, 2014.XU T T. Influence of DBS drag model and turbulence models on CFD simulation of gas-liquid flow[D]. Beijing： University of Chinese Academy of Sciences, 2014.

[23] XU T, JIANG X, YANG N, et al. CFD simulation of internal-loop airlift reactor using EMMS drag model[J]. Particuology, 2015, 19：124-132.

[24] SCHILLER L, NAUMANN Z. A drag coefficient correlation[J]. Z.Ver. Deutsch. Ing., 1935, 77(1)： 318-320.

[25] LUO H, SVENDSEN H F. Theoretical model for drop and bubble breakup in turbulent dispersions[J]. AIChE J., 1996, 66(5)：766-776.

[26] LUO H. Coalescence, breakup and liquid circulation in bubble column reactors[D]. Trondheim, Norway： the Norwegian Institute of Technology, 1993.

[27] LI J. Particle-fluid Two-Phase Flow： the Energy-Minimization Multi-Scale Method[M]. Beijing： Metallurgical Industry Press, 1994.

[28] YANG N, CHEN J H, ZHAO H, et al. Explorations on the multi-scale flow structure and stability condition in bubble columns[J]. Chem.Eng. Sci., 2007, 62(24)： 6978-6991.

[29] YANG N. Chapter five-mesoscale transport phenomena and mechanisms in gas-liquid reaction systems[J]. Advances in Chemical Engineering, 2015, 46： 245-280.

[30] XIAO Q, YANG N, LI J. Stability-constrained multi-fluid CFD models for gas-liquid flow in bubble columns[J]. Chem. Eng. Sci.,2013, 100： 279-292.

[31] YANG N, WU Z, CHEN J, et al. Multi-scale analysis of gas-liquid interaction and CFD simulation of gas-liquid flow in bubble columns[J]. Chem. Eng. Sci., 2011, 66(14)： 3212-3222.

[32] 肖頎, 楊寧. 基于EMMS模型的攪拌釜內氣液兩相流數值模擬[J].化工學報, 2016, 67(7)： 2732-2739.XIAO Q, YANG N. Numerical simulation of gas-liquid flow in stirred tanks based on EMMS model[J]. CIESC Journal, 2016, 67(7)：2732-2739.

[33] MCCLURE D D, NORRIS H, KAVANAGH J M, et al. Validation of a computationally efficient computational fluid dynamics (CFD)model for industrial bubble column bioreactors[J]. Ind. Eng. Chem.Res., 2014, 53(37)： 14526-14543.

[34] XU L, YUAN B, NI H, et al. Numerical simulation of bubble column flows in churn-turbulent regime： comparison of bubble size models[J].Ind. Eng. Chem. Res., 2013, 52(20)： 6794-6802.

[35] 徐琰, 董海峰, 田肖, 等. 鼓泡塔中離子液體-空氣兩相流的CFD-PBM耦合模擬[J]. 化工學報, 2011, 62(10)： 2699-2706.XU Y, DONG H F, TIAN X, et al. CFD-PBM coupled simulation of ionic liquid-air two-phase flow in bubble column[J]. CIESC Journal,2011, 62(10)： 2699-2706.

CFD-PBM simulation with EMMS correctors for bubble column reactors

WANG Jue1,2, YANG Ning1
(1State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences,Beijing 100190, China;2University of Chinese Academy of Sciences, Beijing 100049, China)

The energy-minimization multi-scale (EMMS) model has been introduced to improve the population balance modeling (PBM) of gas-liquid flows. The energy for bubble breakup and coalescence can be obtained from the EMMS model and then used to derive a correction factor for the coalescence rate. This new model is applied in this study to simulate the bubble columns of high flow rates. Simulations using the three different models, namely, the constant-bubble-size model, the CFD-PBM model and the CFD-PBM-EMMS model, are compared with experimental data. The simulation of CFD-PBM-EMMS gives better prediction for bubble size distribution and liquid axial velocity at different heights as well as the overall and local gas holdup. The relative error of global gas holdup reduces to 5% or 15%, and the mean relative error of local gas holdup reduces to 8% or 17% for 0.16 m·s-1or 0.25 m·s-1of superficial gas velocity.

CFD; population balance modeling; bubble column; gas holdup; bubble size distribution

date:2017-01-08.

Prof. YANG Ning, nyang@ipe.ac.cn

supported by the National Natural Science Foundation of China (91434121,91634203).

TQ 021.1

A

0438—1157（2017）07—2667—11

10.11949/j.issn.0438-1157.20170025

2017-01-08收到初稿，2017-03-13收到修改稿。

聯系人：楊寧。

王玨（1991—），男，碩士研究生。

國家自然科學基金項目（91434121, 91634203）。