王玨，楊寧

（1中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2中國科學(xué)院大學(xué)，北京 100049）

基于EMMS方法的鼓泡塔反應(yīng)器CFD及群平衡模擬

王玨1,2，楊寧1

（1中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2中國科學(xué)院大學(xué)，北京 100049）

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

計(jì)算流體力學(xué)；群平衡模型；鼓泡塔；氣含率；氣泡尺徑分布

引 言

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

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

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

1 數(shù)學(xué)模型和數(shù)值模擬

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

1.1 控制方程

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

質(zhì)量守恒方程

動量守恒方程

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

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

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

1.2 湍流模型

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

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

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

湍流動能產(chǎn)生項(xiàng)Gk,m用以下公式計(jì)算得到

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

1.3 相間作用力

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

其中，曳力系數(shù)為

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

1.4 群平衡模型

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

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

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

其中

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

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

其中

1.5 聚并速率修正因子

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

質(zhì)量守恒方程

力平衡方程

對氣液體系的能量關(guān)系的描述如下。

總能量

介尺度能耗

微尺度能耗

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

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

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

破碎核函數(shù)：Luo & Svendsen[25]；聚并核函數(shù)：Luo 對應(yīng)修正因子關(guān)系式[26]。

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

2 建模對象和方法

2.1 結(jié)構(gòu)模型設(shè)置和網(wǎng)格無關(guān)性驗(yàn)證

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

圖2 網(wǎng)格結(jié)構(gòu)Fig. 2 Mesh configuration

表1 網(wǎng)格無關(guān)性驗(yàn)證模型網(wǎng)格尺寸及數(shù)目Table 1 Mesh size and cells adopted for grid independence test

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

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

2.2 數(shù)學(xué)模型求解設(shè)置：

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

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

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

離散法處理各個(gè)氣泡粒級

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

同時(shí)氣泡個(gè)數(shù)和氣含率之間存在

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

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

3 模擬結(jié)果與討論

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

表2 不同模型算例的設(shè)置Table 2 Different model settings

3.1 氣含率

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

圖5 整體氣含率模擬數(shù)據(jù)和實(shí)驗(yàn)結(jié)果[33]的對比Fig. 5 Comparison between different models and experiments[33]

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

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

圖7 局部氣含率模擬數(shù)據(jù)和實(shí)驗(yàn)結(jié)果[33]的對比Fig.7 Comparison of local gas holdup between different models and experiments[33]

圖7是不同高度處的局部氣含率沿徑向的分布。McClure等[33]選取了鼓泡塔中段和接近相界面處的兩個(gè)截面記錄了實(shí)驗(yàn)數(shù)據(jù)。case 1、case 2和case 3的模擬數(shù)據(jù)同實(shí)驗(yàn)結(jié)果相比，氣含率的徑向分布結(jié)果均呈現(xiàn)出中間高，邊壁低的分布規(guī)律。分布器非均勻分布結(jié)構(gòu)使得軸心處的氣含率略有偏低，3個(gè)case也都捕捉到了這一現(xiàn)象。與整體氣含率相同的是，局部氣含率的結(jié)果也表現(xiàn)出了case 1明顯高估了氣含率，case 2、case 3在誤差允許的范圍內(nèi)都可以較為可靠地反映鼓泡塔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處氣泡尺徑分布的模擬值和實(shí)驗(yàn)值的對比，其中case 2中的氣泡峰值位于39 mm的粒級。小于39 mm的其他粒級的氣泡，分布較為平緩。分析原因：① case 2 中的氣泡粒級受到1～39 mm粒級區(qū)間設(shè)置的限制和聚并核函數(shù)對聚并速率的高估，使粒徑的峰值為 39 mm，而且從37 mm到39 mm的粒級間陡然上升，由于原始聚并模型高估了聚并速率[11]，模擬得到的粒徑分布顯示出大氣泡較多的趨勢；② 從曲線整體分析，1～37 mm的氣泡尺徑分布較為平緩，與實(shí)驗(yàn)結(jié)果相比，說明了模擬中氣泡的聚并被高估了。case 3嘗試使用修正因子抑制聚并函數(shù)，修正氣泡聚并速率。修正后的模型中，不同尺徑氣泡的體積分布范圍集中于5～15 mm，9 mm的氣泡具有最高的體積分?jǐn)?shù)。兩種不同的模型對比，修正后的群平衡模型預(yù)測的粒徑分布更加符合實(shí)際。實(shí)驗(yàn)中體積分?jǐn)?shù)最高的氣泡粒級為 5 mm，與模型預(yù)測結(jié)果仍有一定差距。McClure等[33]分析實(shí)驗(yàn)結(jié)果認(rèn)為，在H/d=2、H/d=3截面上，氣泡尺徑分布在實(shí)驗(yàn)允許的誤差范圍內(nèi)是相同的，而帶有修正因子的模擬結(jié)果也反映了這一點(diǎn)。因此，在H/d≥2截面上，氣泡的聚并和破碎達(dá)到了動態(tài)平衡。

圖8 氣泡尺徑分布模擬數(shù)據(jù)和實(shí)驗(yàn)結(jié)果[33]的對比Fig.8 Comparison of bubble size distribution between different models and experiments[33]

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

3.3 軸向液速

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

4 結(jié) 論

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

（1）高氣速條件下的鼓泡塔中發(fā)生強(qiáng)烈的湍動，氣泡的破碎和聚并導(dǎo)致入口處氣泡尺徑和塔中穩(wěn)定鼓泡段的氣泡尺徑有著很大的區(qū)別，此時(shí)采用均一氣泡尺徑模型不能準(zhǔn)確模擬塔內(nèi)的所有流場信息?；贓MMS方法修正群平衡模型，考慮了氣泡的聚并和破碎，獲得了具有穩(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聚并模型需要經(jīng)過修正才能準(zhǔn)確地描述反應(yīng)器中的實(shí)際氣泡分布。一般情況下傳統(tǒng)的聚并核函數(shù)模型均高估了聚并速率。使用EMMS-PBM 耦合方法計(jì)算出的修正因子有效地抑制了聚并，得到單峰分布的氣泡尺徑曲線。與未修正前的曲線相比，模擬預(yù)測更加接近實(shí)驗(yàn)結(jié)果。

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

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

符 號 說 明

b——破碎核函數(shù)，s-1

C——聚并修正因子

c——聚并核函數(shù)，m3·s-1

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

D——?dú)馀葜睆剑琺

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

d32——Sauter平均直徑，m

FD——?dú)馀萑阂妨?，N·m-3

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

M ——相間動量交換項(xiàng)，kg·m-2·s-2

N——數(shù)密度函數(shù)

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

圖10 軸向液速模擬結(jié)果[33]與實(shí)驗(yàn)數(shù)據(jù)的對比Fig.10 Comparison of upward liquid velocity between different models and experiment[33]

Nsurf——?dú)馀荼砻婺?，m2·s-3

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

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

P ——壓強(qiáng)，Pa

Re ——Reynolds數(shù)

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

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

V ——?dú)馀蒹w積，m3

We ——Weber數(shù)

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

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

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

λ ——湍流渦尺度，m

μ ——動力學(xué)黏度，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]. 化學(xué)進(jìn)展, 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] 肖頎. 攪拌釜及鼓泡塔內(nèi)氣液兩相流的多尺度模型[D]. 西安： 西安交通大學(xué), 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]. 北京： 中國科學(xué)院大學(xué), 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模型的攪拌釜內(nèi)氣液兩相流數(shù)值模擬[J].化工學(xué)報(bào), 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]. 化工學(xué)報(bào), 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收到修改稿。

聯(lián)系人：楊寧。

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

國家自然科學(xué)基金項(xiàng)目（91434121, 91634203）。