賈翔飛， 錢嘉澍， 吳幼青， 陳劍佩

（華東理工大學(xué) 1. 化學(xué)工程聯(lián)合國家重點實驗室；2. 煤氣化及能源化工教育部重點實驗室，上海 200237）

我國的自然資源呈現(xiàn)貧油、少氣、富煤的現(xiàn)象，每年需要進口大量的石油。為了保障能源安全，我國對煤液化技術(shù)進行了技術(shù)儲備并且成功將該技術(shù)工業(yè)化。煤液化技術(shù)的關(guān)鍵設(shè)備是液化反應(yīng)器，其本質(zhì)上屬于鼓泡塔反應(yīng)器。該反應(yīng)器具有良好的混合性能并且生產(chǎn)能力大、結(jié)構(gòu)簡單、易于操作[1-2]。

煤液化反應(yīng)器中涉及到氫氣（氣）、供氫溶劑（液）、煤粉顆粒（固）的流動問題，但反應(yīng)器中的煤粉顆粒較小并且濃度較低，煤粉顆粒實際上是伴隨著液相一起運動[3]，因此通常將煤粉顆粒與液相作為混合均勻的漿相處理，將氣-液-固三相流動簡化為氣-漿兩相流動。根據(jù)反應(yīng)器操作的表觀氣速不同，反應(yīng)器內(nèi)的狀態(tài)可以分為均勻鼓泡狀態(tài)、過渡狀態(tài)、湍動鼓泡狀態(tài)。在均勻鼓泡狀態(tài)下，鼓泡塔內(nèi)的氣泡尺寸小而均勻，氣含率與表觀氣速呈現(xiàn)近似線性關(guān)系[4]。工業(yè)生產(chǎn)過程中的鼓泡塔反應(yīng)器一般處于湍動鼓泡狀態(tài)[5]，在該狀態(tài)下，流動結(jié)構(gòu)由于氣泡大量的聚并和破碎行為而變得非常復(fù)雜。強烈的聚并和破碎行為使得鼓泡塔內(nèi)的氣泡尺寸分布范圍較寬，并對鼓泡塔內(nèi)與反應(yīng)器設(shè)計和放大有關(guān)的氣含率、液速、熱量和質(zhì)量傳遞速率等參數(shù)產(chǎn)生影響[1]。因此，反應(yīng)器內(nèi)的氣泡破碎聚并以及基礎(chǔ)流體力學(xué)行為對反應(yīng)器的設(shè)計放大以及安全運行有著重要影響。

在對反應(yīng)器進行設(shè)計和放大時，計算流體動力學(xué)（CFD）逐漸成為一種重要的工具[2,6-8]。在綜合考慮計算量和計算精度的前提下，通常采用Eulerian-Eulerian 方法[9]對鼓泡塔進行模擬。使用該方法研究氣液鼓泡塔的流體力學(xué)行為關(guān)鍵在于合理地描述氣-液相間作用力和氣泡的大小分布。相間作用力包括曳力、升力、壁面潤滑力及湍流擴散力。曳力是由氣液之間的摩擦產(chǎn)生的，曳力的大小決定了氣泡在液相中的上升速度并最終影響鼓泡塔內(nèi)的氣含率的大小，而升力、壁面潤滑力和湍流擴散力會影響氣含率在徑向上的分布。

由于氣泡與氣泡之間存在復(fù)雜的相互作用，鼓泡塔內(nèi)的氣泡存在局部氣泡直徑分布。早期的研究者通常將該氣泡直徑分布簡化為一個恒定的當(dāng)量氣泡直徑[10-12]。在均勻鼓泡狀態(tài)下，氣泡之間的相互作用較弱，氣泡的尺寸分布較窄，氣泡直徑均一，通過調(diào)節(jié)曳力模型及氣泡直徑會獲得與實驗值較吻合的結(jié)果。但在湍動鼓泡狀態(tài)下，氣泡之間強烈的相互作用造成氣泡較寬的尺寸分布，因而上述對氣泡直徑的簡化無法獲得與實驗值吻合的結(jié)果。

由于氣泡之間的相互作用對鼓泡塔內(nèi)的流體力學(xué)參數(shù)有重要影響，因而群體平衡模型（PBM）逐漸被應(yīng)用于描述鼓泡塔內(nèi)的氣泡尺寸的變化及分布，并與CFD耦合以對鼓泡塔進行模擬[13-16]。多年來，人們進行大量的工作來建立合理的聚并及破碎模型[16-19]，并希望將其應(yīng)用于CFD-PBM耦合模型中。

本文將一個新穎的氣泡破碎模型[19]植入PBM并與歐拉雙流體模型進行耦合。首先在ANSYS Fluent 平臺中將CFD模型所需要的相間作用力模型通過用戶自定義函數(shù)（UDF）加入，然后將PBM 模型中聚并、破碎模型源項通過UDF植入并與CFD模型耦合建立CFD-PBM 模型。采用該模型對實驗室尺寸的圓柱鼓泡塔進行二維軸對稱模擬，并與文獻(xiàn)中的實驗數(shù)據(jù)進行比較。

1 模型方法

1.1 雙流體模型的控制方程及湍流模型

文獻(xiàn)[9]中提到，采用重整化群k-ε（RNGk-ε）湍流模型可以更好地描述鼓泡塔中的流動。本文采用RNGk-ε湍流模型計算液相中的湍動能和湍動耗散率。

1.2 相間作用力模型

1.2.1 曳力 曳力的大小會影響氣泡在液相中的上升速度，并進一步影響鼓泡塔內(nèi)氣含率的大小，曳力通常采用式（3）計算[14]：

1.2.2 升力 液相中速度梯度的存在會使氣泡受到一個橫向升力，這會對流體力學(xué)參數(shù)的徑向分布產(chǎn)生影響，升力系數(shù)的計算采用Tomiyama 模型[20]。

1.2.3 壁面潤滑力 氣泡受到的壁面潤滑力采用Tomiyama 模型[20]。

1.2.4 湍動耗散力 湍動耗散力采用Lopez de Bertodano模型[21]。

1.3 群體平衡模型

1.3.1 群體平衡方程（PBE） 在鼓泡塔中，存在著數(shù)量非常龐大的氣泡，這些氣泡會發(fā)生聚并和破碎，最終使鼓泡塔中的氣泡存在一定的尺寸分布。PBE（式（15））以及合理的聚并、破碎源項的引入使得鼓泡塔中的氣泡尺寸分布得以描述。

其中：Si為第i組氣泡中PBM 的源項[22]，在目前的研究工作中，成核及生長速率被忽略，只考慮聚并生成、聚并消失、破碎生成及破碎消失項（分別表示為BC,i，DC,i，BB,i，DB,i）。各源項的計算采用文獻(xiàn)[22]的方法。

1.3.2 聚并核函數(shù) 鼓泡塔中引起氣泡之間相互碰撞的機理通常有3種，分別是液相湍流引起的碰撞、氣泡上升速度差異引起的碰撞以及大氣泡尾渦引起的碰撞。但不是每次氣泡碰撞都會使兩個氣泡發(fā)生聚并，所以通過引入氣泡聚并效率的概念來表示兩個氣泡聚并的可能性。最終，可以通過式（17）來計算直徑為di和dj的氣泡的聚并速率。

鼓泡塔中總的聚并速率可以認(rèn)為是3種機理聚并（湍流渦引起的聚并、氣泡尾渦引起的聚并、氣泡上升速度不同引起的聚并）速率的線性疊加。

（1）湍流渦引起的聚并

①碰撞頻率

②聚并效率

根據(jù)文獻(xiàn)[25]提出的液膜理論，氣泡的相互作用時間和氣泡聚并所需時間的相對大小決定著氣泡聚并的概率。文獻(xiàn)[25-26]給出了直徑為di、dj的氣泡的聚并效率。

（2）由氣泡尾渦引起的聚并

①碰撞頻率

鼓泡塔中氣泡在尺寸較大時會由球形變成球帽形，并在氣泡的尾部形成渦流。此時若其他氣泡進入該大氣泡的尾渦影響區(qū)域，就會被該氣泡的尾渦夾帶而加速上升。該過程會使大氣泡和夾帶氣泡發(fā)生碰撞并聚并，文獻(xiàn)[16]提出尾渦引起的氣泡碰撞頻率為：

由于小氣泡的尾渦相對較小，因此這里只考慮大氣泡尾渦引起的聚并，所以在碰撞頻率中加入因子Θ。

②聚并效率

尾渦引起的聚并效率采用文獻(xiàn)[27]提出的關(guān)聯(lián)式。

（3）由于氣泡上升速度不同引起的聚并

①碰撞頻率

由氣泡上升速度不同引起的聚并模型與液相湍流引起的聚并模型相似[23]，其氣泡碰撞頻率表達(dá)式為：

②聚并效率

由于氣泡上升速度不同造成的聚并，文獻(xiàn)[14, 24]將聚并效率設(shè)置為0.5，即

1.3.3 破碎核函數(shù) 對于氣泡的破碎過程，本文只考慮湍流渦引起的破碎，并且破碎過程為二元破碎，即每個母氣泡破碎成兩個子氣泡，故氣泡破碎速率可采用式（26）進行計算[16,18]。

（1）氣泡破碎碰撞頻率

目前，一些研究者[17,28-30]認(rèn)為大尺寸的湍流渦會把部分能量傳遞給氣泡從而使氣泡發(fā)生破碎。因此將λ≤d的渦和λ>d的渦與氣泡的碰撞頻率分別進行計算[19]。

氣泡破碎時必須同時滿足能量約束和壓力約束[13,15-16,31]。當(dāng)λ>d時[28]，湍流渦中的能量只能部分作用于氣泡并導(dǎo)致氣泡的破碎，因此在量綱為一湍流渦能量中引入修正因子Ψ。

對氣泡頸處的受力進行分析最終得到壓力破碎因子為：

能量約束條件的推導(dǎo)可見文獻(xiàn)[18]，能量約束因子為：

1.4 模型方程求解

基于文獻(xiàn)[33]設(shè)置相關(guān)實驗并進行二維軸對稱模擬,采用空氣-水作為實驗介質(zhì)，鼓泡塔的相關(guān)幾何模型見圖1（a）。鼓泡塔網(wǎng)格采用ANSYSICEM CFD劃分并對進口處以及壁面處進行加密處理。最終經(jīng)過網(wǎng)格無關(guān)化檢驗后確定四邊形網(wǎng)格數(shù)量為9 581，計算過程中的時間步長取0.002 s。氣體均勻地從底部r/R<0.9的中心區(qū)域網(wǎng)格進入（圖1（b））。實驗的氣速范圍為0.02～0.10 m/s，覆蓋了均勻鼓泡流和湍動鼓泡流。所有計算均在ANSYSFluent 平臺上進行，壓力-速度方程使用SIMPLE算法處理，動量方程采用QUICK 格式，體積分?jǐn)?shù)方程等均采用一階迎風(fēng)格式離散，松弛因子均為默認(rèn)值。模擬中所用到的相間作用力模型均采用UDF加入Fluent 求解平臺。

圖1 鼓泡塔反應(yīng)器幾何模型（a）及網(wǎng)格截面（b）Fig.1 Geometry (a)and mesh (b)of bubble column reactor

鼓泡塔反應(yīng)器的入口及出口邊界條件為速度進口和壓力出口，壁面邊界條件為無滑移邊界。PBE方程的求解采用離散法，同時破碎和聚并模型均采用UDF實現(xiàn)。本文氣泡尺寸為0 .5≤db≤40.3 mm，分為20組，各組氣泡大小的計算如式（36）所示。

2 結(jié)果與討論

2.1 概述

采用CFD-PBM 耦合模型對鼓泡塔反應(yīng)器內(nèi)的氣含率及軸向液速分布情況進行模擬，氣泡聚并及破碎模型通過UDF接口植入ANSYSFluent 中。其中鼓泡塔反應(yīng)器入口表觀氣速分別為0.02 m/s和0.10 m/s，并與相應(yīng)表觀氣速下的實驗數(shù)據(jù)[33]進行對比。

2.2 破碎模型驗證

其中：L和T分別為長度尺度和時間尺度[34]，可分別通過式（39）、（40）進行計算。

圖2 量綱為一破碎速率預(yù)測值與實驗值對比Fig.2 Predicted and experimental values of dimensionless breakup rate

圖2中實線代表模型預(yù)測的量綱為一破碎速率，數(shù)據(jù)點為文獻(xiàn)[35-37]中的實驗測量值。從圖2中可以看出，在量綱為一氣泡直徑較大時模型預(yù)測值略小于實驗值，在量綱為一氣泡直徑較小時模型預(yù)測值與實驗值基本吻合。

圖3比較了子氣泡分布的模型預(yù)測值與實驗測定值，其中實驗數(shù)據(jù)來自文獻(xiàn)[32]。實線為模型預(yù)測的子氣泡分布概率密度曲線，灰色直方圖為概率密度曲線在不同區(qū)間上的積分值，紅色直方圖為相應(yīng)區(qū)間上的實驗測量值。從圖中可以看出由于模型采用了壓力約束條件，在氣泡破碎體積分?jǐn)?shù)fv趨于0時子氣泡分布的概率密度β（fv,d）曲線也趨于0。同時，從圖3中也可以看出，模型預(yù)測的結(jié)果為氣泡更傾向于不均勻破碎，這也與實驗測定的結(jié)果相吻合[32]。

圖3 子氣泡分布預(yù)測值與實驗值Fig.3 Predicted and experimental valuesof daughter bubble size distribution

2.3 升力的影響

升力是由于液相速度梯度差使得氣泡在徑向方向上受到的力，其決定了氣泡在徑向的運動方向。升力的加入會使大氣泡向反應(yīng)器中心移動，因此會對氣含率和軸向液速在徑向上的分布產(chǎn)生影響。圖4和圖5分別示出了添加升力對氣含率和軸向液速在徑向分布的影響，圖4、5中實線分別為未添加升力時的氣含率、軸向液速徑向分布。從圖4可以看出，當(dāng)不加入升力時，氣含率在徑向上的分布更加平坦，在鼓泡塔中心處模擬值小于實驗測量值，在靠近壁面處模擬值遠(yuǎn)大于實驗測量值。由圖5可以看出，升力的加入使得徑向速度梯度變大，模擬值與實驗值更加吻合?？傮w來說，升力對氣泡在鼓泡塔內(nèi)流體力學(xué)參數(shù)的徑向分布十分重要，升力的加入使徑向氣含率的分布與實驗值更加吻合。

2.4 氣含率

圖6示出了不同表現(xiàn)氣速下鼓泡塔反應(yīng)器內(nèi)氣含率及軸向液速分布云圖。由圖6中氣含率分布云圖和圖4中氣含率徑向分布曲線可以看出，鼓泡塔內(nèi)氣含率沿徑向方向逐漸降低。由于鼓泡塔內(nèi)存在液相的循環(huán)，液體在中心處上升而在壁面處下降，這種液相的循環(huán)使得氣含率在鼓泡塔的中心較高，而在壁面處較低。從氣含率云圖與徑向分布圖可以看出，在低表觀氣速（0.02 m/s）下，氣含率的徑向分布比較均勻，但在高表觀氣速（0.10 m/s）下，鼓泡塔中心與壁面處氣含率的徑向分布差異較大。

由圖4可以看出，在高表觀氣速下，鼓泡塔壁面處氣含率的模擬值大于實驗值；在低表觀氣速（0.02 m/s）下，氣含率模擬值在徑向略小于實驗值。從整體上看，在低表觀氣速下，目前模型的計算值低估了氣含率。

圖4 不同表觀氣速下氣含率徑向分布Fig.4 Radial profiles of gas holdup at different superficial gas velocities

圖5 不同表觀氣速下軸向液速徑向分布Fig.5 Radial profiles of axial liquid velocity at different superficial gasvelocities

2.5 軸向氣、液速

由圖6中軸向液速分布云圖及圖5中軸向液速徑向分布圖可以看出，鼓泡塔內(nèi)的軸向液速沿徑向方向逐漸降低。由圖6可以看出，在徑向方向上，中心處的軸向液速最大，沿徑向方向，軸向液速逐漸降低并變?yōu)樨?fù)值。此時的負(fù)值表示液體的流動方向變?yōu)榕c原來流動方向相反。造成這種現(xiàn)象的原因是鼓泡塔內(nèi)的氣相在浮力的驅(qū)動下沿鼓泡塔上升，進而造成鼓泡塔內(nèi)的液相產(chǎn)生這種軸向循環(huán)。

圖5中將現(xiàn)有模型的模擬值與實驗值進行了比較。在高表觀氣速下，即處于湍動鼓泡狀態(tài)時[9]，模擬值與實驗值吻合良好，但在中心處略微低估了軸向液速。在低表觀氣速下，即均勻鼓泡狀態(tài)時，模擬值與實驗值大致吻合，但在鼓泡塔中心處，模擬值要高于實驗值。

圖6 不同表觀氣速下鼓泡塔反應(yīng)器內(nèi)氣含率及軸向液速分布云圖Fig.6 Gas holdup and axial liquid velocity contours at different superficial gasvelocitiesin bubblecolumn reactor

圖7 不同表觀氣速下氣泡軸向速度徑向分布Fig.7 Radial profiles of bubble axial velocity at different superficial gasvelocities

圖7所示為鼓泡塔內(nèi)不同表觀氣速下氣泡軸向速度的徑向分布，可以看出，隨表觀氣速增加，氣泡的軸向速度增加，并且氣泡軸向速度徑向分布的不均勻性也在增加。結(jié)合圖5和圖7也可以看出，氣液之間的相對速度變大。這是因為隨著表觀氣速的增加，反應(yīng)器狀態(tài)由均勻鼓泡狀態(tài)變?yōu)橥膭庸呐轄顟B(tài)[4]，氣泡大小分布變得更寬，反應(yīng)器內(nèi)大氣泡數(shù)量增加[16]，大氣泡尾渦產(chǎn)生的加速效應(yīng)使得氣液之間的相對速度變大。

2.6 湍動耗散率

氣泡的聚并與破碎速率均與湍動耗散率有關(guān)，并且氣泡的聚并與破碎速率會對局部氣泡直徑產(chǎn)生影響，進而影響氣液相間作用力的大小，因此合理地描述鼓泡塔內(nèi)的液相湍動耗散率對于流場的準(zhǔn)確模擬非常重要。

圖8所示為不同表觀氣速下液相的湍動耗散率的徑向分布?？梢钥闯?，湍動耗散率與表觀氣速呈正比關(guān)系。同時，由于壁面的存在以及流體黏性的作用使得壁面處液相的速度梯度較大，因此壁面處的湍動耗散率較大。

圖8 不同表觀氣速下湍動耗散率徑向分布Fig.8 Radial profiles of turbulent dissipation rate at different superficial gasvelocities

2.7 氣泡尺寸分布

圖9所示為鼓泡塔H/D=3截面處不同表觀氣速下氣泡尺寸分布，該圖橫坐標(biāo)為氣泡直徑，縱坐標(biāo)為基于氣泡體積處理得到的氣泡尺寸分布概率密度函數(shù)。從圖中可以看出，隨著表觀氣速增加，H/D=3截面處平均氣泡尺寸分布曲線變得更加低矮，同時大氣泡的體積分?jǐn)?shù)增加。

圖9 H/D=3處不同表觀氣速下氣泡尺寸分布Fig.9 Bubble size distribution at H/D=3 at different superficial gasvelocities

圖10所示為H/D=3處，r/R=0和r/R=0.5處不同表觀氣速下的氣泡尺寸分布。從圖10（a）可以明顯看出，隨著表觀氣速增加，鼓泡塔由均勻鼓泡狀態(tài)轉(zhuǎn)變?yōu)橥膭庸呐轄顟B(tài)，鼓泡塔中心位置的氣泡尺寸分布變得更寬，大氣泡的體積分?jǐn)?shù)增加。這主要是由于湍動鼓泡狀態(tài)下，較大的表觀氣速造成鼓泡塔內(nèi)的氣泡數(shù)目增加。單位體積內(nèi)氣泡數(shù)量的增加使得氣泡之間的相互作用更加頻繁，在液相湍流的誘導(dǎo)下氣泡聚并產(chǎn)生更多的大氣泡，而大氣泡數(shù)量的增加也使得更多的氣泡因為氣泡尾渦效應(yīng)和氣泡上升速度差效應(yīng)而發(fā)生聚并。另一方面，鼓泡塔內(nèi)聚并產(chǎn)生的大氣泡擁有更大的碰撞截面積，因此會有更多的湍流渦體與其發(fā)生相互作用。與此同時，高表觀氣速下湍動耗散率增加使得湍流渦體的湍動能增加，最終表現(xiàn)為氣泡的破碎速率增加。這些因素的綜合作用使得氣泡尺寸的分布范圍變寬。同時，在升力的作用下，大氣泡趨向于向鼓泡塔中心運動，故可以看到高表觀氣速下，r/R=0處氣泡尺寸分布較寬且大氣泡數(shù)量較多。從圖10（b）可以看出，沿徑向方向遠(yuǎn)離鼓泡塔中心處，高表觀氣速與低表觀氣速的氣泡尺寸分布差異在減小。但高表觀氣速下氣泡尺寸分布范圍依然大于低表觀氣速下的范圍。

圖10 不同表觀氣速下r/R=0（a）和r/R=0.5（b）的氣泡尺寸分布Fig.10 Bubblesizedistribution at r/R=0(a),and r/R=0.5(b)at different superficial gasvelocities

圖11所示為表觀氣速為0.10 m/s時不同徑向位置處的氣泡尺寸分布，從圖中可以非常清晰地看出，沿徑向靠近鼓泡塔中心，大氣泡的體積分?jǐn)?shù)在增加，小氣泡的體積分?jǐn)?shù)在減少。

圖11 徑向位置處氣泡尺寸分布Fig.11 Bubblesizedistribution at different radial positions

3 結(jié) 論

將一個新型的破碎模型通過UDF植入ANSYS Fluent 軟件中，并將聚并和破碎模型與雙流體模型進行雙向耦合建立一個通用的CFD-PBM 耦合模型，在考慮了不同機理對氣泡破碎和聚并的影響以及氣液相間作用力后，對處于均勻鼓泡狀態(tài)與湍動鼓泡狀態(tài)下的鼓泡塔流體力學(xué)參數(shù)進行研究。通過數(shù)值模擬得出以下結(jié)論：

（1）在對鼓泡塔進行模擬時，升力是否加入對于氣含率和軸向液速等流體力學(xué)參數(shù)的軸向分布有重要影響，升力的加入使模擬數(shù)據(jù)在徑向上的分布與實驗值吻合得更好。

（2）在對破碎模型的氣泡破碎速率和子氣泡分布預(yù)測值進行驗證之后，將模型在不同表觀氣速下的氣含率、軸向液速模擬值與實驗值進行了對比。結(jié)果顯示CFD-PBM 耦合模型的計算值與實驗值吻合較好，這也證明了該模型可以用于湍動鼓泡狀態(tài)下鼓泡塔的流體力學(xué)模擬。

（3）鼓泡塔內(nèi)液相湍動耗散率隨表觀氣速的增加而增加。同時，壁面的存在以及流體黏性的作用使得壁面處液相的速度梯度較大，因此壁面處的湍動耗散率較大。

（4）隨表觀氣速增加，鼓泡塔內(nèi)的氣泡數(shù)目增加，單位體積內(nèi)氣泡數(shù)量的增加使得氣泡之間的相互作用更加頻繁，氣泡聚并速率增加。另一方面，鼓泡塔內(nèi)聚并產(chǎn)生的大氣泡擁有更大的碰撞截面積，因此會有更多湍流渦體與其發(fā)生相互作用。與此同時，高表觀氣速下湍動耗散率增加使得湍流渦體的湍動能增加，最終表現(xiàn)為氣泡的破碎速率增加。這些因素的綜合作用使得氣泡尺寸的分布范圍變寬。