劉鑫，張煜，張麗，靳海波

（1北京石油化工學(xué)院化學(xué)工程學(xué)院，北京 102617；2中科合成油技術(shù)有限公司研發(fā)中心，北京 101407）

基于氣泡群相間作用力模型的加壓鼓泡塔流體力學(xué)模擬

劉鑫1,2，張煜2，張麗2，靳海波1

（1北京石油化工學(xué)院化學(xué)工程學(xué)院，北京 102617；2中科合成油技術(shù)有限公司研發(fā)中心，北京 101407）

目前，多數(shù)文獻(xiàn)報(bào)道了冷態(tài)加壓湍動鼓泡塔內(nèi)流動特征，并且通過實(shí)驗(yàn)數(shù)據(jù)回歸相關(guān)經(jīng)驗(yàn)關(guān)聯(lián)式。然而，此類關(guān)聯(lián)式適用范圍有限，難以直接外推到工業(yè)鼓泡塔反應(yīng)器條件。因此，在FLUENT平臺上建立了基于氣泡群相間作用力的、動態(tài)二維加壓鼓泡塔計(jì)算流體力學(xué)模型。通過數(shù)值模擬考察了操作壓力為0.5～2.0 MPa，表觀氣速為0.20～0.31 m·s?1，內(nèi)徑0.3 m鼓泡塔內(nèi)流場特性參數(shù)分布，并且與冷態(tài)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。結(jié)果表明，采用修正后的氣泡群曳力模型、徑向力平衡模型以及壁面潤滑力模型描述氣泡群相間作用力，能夠較為準(zhǔn)確地反映平均氣含率和氣含率徑向分布隨操作壓力和表觀氣速變化的規(guī)律。

鼓泡塔；氣液兩相流；計(jì)算流體力學(xué)；氣泡群曳力模型；徑向力平衡模型；壁面潤滑力模型

引 言

煤間接液化技術(shù)作為煤炭清潔高效利用的重要手段之一，它是以煤炭為原料制取目標(biāo)產(chǎn)品為液態(tài)燃料的技術(shù)。隨著以費(fèi)托合成工藝為核心的煤炭間接液化技術(shù)的重要性日益突顯，作為費(fèi)托合成反應(yīng)的載體，費(fèi)托合成反應(yīng)器得到了快速發(fā)展。費(fèi)托合成反應(yīng)器先后經(jīng)歷了固定床、流化床和漿態(tài)床 3個(gè)不同階段，其中采用漿態(tài)床反應(yīng)器生產(chǎn)的柴油，具有高能量密度、低排放、品質(zhì)好等優(yōu)點(diǎn)，因此，漿態(tài)床已成為煤制油進(jìn)行費(fèi)托合成反應(yīng)的主流反應(yīng)器[1]。

工業(yè)鼓泡塔反應(yīng)器通常帶壓操作，并且表觀氣速較高，處于湍流鼓泡狀態(tài)。目前對于加壓鼓泡塔的研究大多集中于冷態(tài)實(shí)驗(yàn)測量和流動規(guī)律的總結(jié)，根據(jù)實(shí)驗(yàn)數(shù)據(jù)回歸經(jīng)驗(yàn)關(guān)聯(lián)式，缺乏能夠詳細(xì)描述加壓鼓泡塔的流體力學(xué)模型。Krishna等[2]根據(jù)加壓氣泡界面的Kelvin-Helmholtz不穩(wěn)定性分析和氣泡曳力模型，將壓力影響作為氣相密度修正項(xiàng)ρG/ρG,1atm引入曳力系數(shù)計(jì)算公式。建立了雙氣泡二維軸對稱CFD模型，用CFX軟件模擬了塔徑0.15 m加壓湍動鼓泡塔充分發(fā)展段的液相軸向擴(kuò)散系數(shù)及氣含率，并與冷態(tài)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。在高氣速下，該模型計(jì)算的氣含率明顯高于實(shí)驗(yàn)測量值。Chen等[3]建立了結(jié)合氣泡群體平衡方程（PBM）的 3D雙流體和混合物模型，用 FLUENT模擬了塔徑為0.162 m加壓鼓泡塔流場?？紤]到Krishna等[2]的模型在高氣速下對氣含率高估，認(rèn)為造成偏差的原因可能是，高氣速下兩相湍流和尾渦的卷吸作用影響顯著。為了改進(jìn)模型的計(jì)算準(zhǔn)確度，Chen等[3]引入(ρG/ρG,1atm)0.25修正氣泡曳力項(xiàng)。該模型采用耦合PBM的3D模擬，為保證計(jì)算準(zhǔn)確度，離散網(wǎng)格數(shù)量隨鼓泡塔尺寸提高而快速增加，需要大量計(jì)算資源。以現(xiàn)有計(jì)算機(jī)運(yùn)算能力，仍然難以采用此模型進(jìn)行大型加壓鼓泡塔反應(yīng)器 3D模擬和放大設(shè)計(jì)。

Joshi[4]認(rèn)為氣泡所受的橫向作用力決定其沿塔徑向的運(yùn)動模式，是湍動鼓泡塔形成氣含率不均勻分布的重要影響因素；Tabid等[5]通過3D動態(tài)模擬鼓泡塔流場，比較了不同相間作用力的影響，認(rèn)為曳力、升力和湍流擴(kuò)散力對氣含率和軸向液速分布影響顯著。文獻(xiàn)[6-8]通過引入氣泡徑向力平衡機(jī)制，確定氣含率徑向分布，一維及二維數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合較好。綜上所述，二維軸對稱模型耦合適當(dāng)?shù)南嚅g作用力模型，能夠在可接受的計(jì)算代價(jià)下，得到較為準(zhǔn)確的常壓鼓泡塔流場模擬結(jié)果。本文擬建立動態(tài)二維軸對稱加壓湍動鼓泡塔模型，包括徑向力平衡機(jī)制、氣泡群曳力模型和壁面潤滑力，根據(jù)內(nèi)徑0.3 m加壓冷態(tài)鼓泡塔冷態(tài)實(shí)驗(yàn)數(shù)據(jù)，擬合氣泡群曳力系數(shù)參數(shù)，使模型具備預(yù)測加壓湍動鼓泡塔反應(yīng)器內(nèi)多相流場特性，為工業(yè)漿態(tài)床反應(yīng)器放大和內(nèi)構(gòu)件優(yōu)化提供指導(dǎo)。

1 數(shù)學(xué)模型

1.1 守恒方程

本文所提出的雙流體加壓鼓泡塔模型，采用雷諾時(shí)均的控制方程[9-11]。與操作壓力相比，加壓鼓泡塔靜液壓相對較小，可近似認(rèn)為氣體密度基本不隨軸向高度變化，冷態(tài)條件下可忽略兩相間的傳熱和傳質(zhì)作用，將相間作用力以體積源項(xiàng)的方式，添加至動量守恒方程組中，控制方程可表示如下：

連續(xù)性方程為

文獻(xiàn)[9, 11]給出了氣液兩相表達(dá)形式一致的動量守恒方程

文獻(xiàn)[12-17]報(bào)道采用標(biāo)準(zhǔn)k-ε模型能夠描述鼓泡塔內(nèi)液相湍流場，本文也采用該模型計(jì)算液相湍流動能k和湍流動能耗散率ε的分布，以及湍流黏度μt，模型方程為

式中，參數(shù)Cμ=0.09，C1ε=1.44，C2ε=1.92，σk=1.0，σε=1.3。

1.2 相間作用力模型

鼓泡塔雙流體模型中，主要通過體積平均的相間作用力源項(xiàng)描述氣液相間動量傳遞[4-5,10-11]。本文所提出的加壓鼓泡塔流體力學(xué)模型，主要考慮氣泡曳力、升力、湍流擴(kuò)散力和壁面潤滑力。另外，虛擬質(zhì)量力表現(xiàn)為顆粒在加速運(yùn)動時(shí)，使周圍流體也加速運(yùn)動所需要額外的力，Khopkar等[18]認(rèn)為虛擬質(zhì)量力在氣體分布器附近更重要，而在流動主體區(qū)域的作用遠(yuǎn)小于曳力，且本文通過實(shí)際工況模擬發(fā)現(xiàn)，添加虛擬質(zhì)量力模型對模擬結(jié)果幾乎無影響，因此本文忽略了虛擬質(zhì)量力。

1.2.1 曳力 鼓泡塔內(nèi)氣泡曳力是相間作用力研究的重點(diǎn)之一，也是現(xiàn)有各種流體力學(xué)模型主要考慮的相間作用之一，甚至是許多模型中唯一考慮的相間作用力。對于鼓泡塔雙流體模型，動量守恒方程中的曳力源項(xiàng)通常表述為單位體積內(nèi)氣泡群所產(chǎn)生的合曳力，表示為

在加壓湍動鼓泡流狀態(tài)下，氣含率相對較高，氣泡在浮升過程中會受到附近氣泡的影響，需要引入氣泡群曳力系數(shù)進(jìn)行修正，文獻(xiàn)中通常將氣泡群曳力系數(shù)CD與單個(gè)氣泡曳力系數(shù)CD0與氣含率相關(guān)聯(lián)進(jìn)行修正。

文獻(xiàn)中給出的單氣泡曳力系數(shù)模型通常與氣泡Reynolds數(shù)Re相關(guān)。Schiller & Naumann模型常用于氣液兩相流模擬，在鼓泡塔流場數(shù)值模擬過程中穩(wěn)定性較好，模型方程如下

與 Schiller & Naumann模型相比，Morsi & Alexander模型更為完整，其曳力系數(shù)對氣泡Reynolds數(shù)的劃分更詳細(xì)，曳力系數(shù)與Reynolds數(shù)的關(guān)系可表示為

式中，a1、a2、a3為模型參數(shù)，隨Reynolds數(shù)取值區(qū)間的改變而對應(yīng)不同的數(shù)值。但是，該模型的計(jì)算穩(wěn)定性不及其他模型。

Lau等[19]和 Roghair等[20]采用直接數(shù)值模擬（front-tracking model，前緣追蹤模型）單個(gè)氣泡及氣泡群浮升過程，通過模擬結(jié)果擬合得到氣泡群曳力系數(shù)與氣含率和Eo數(shù)的關(guān)聯(lián)式

但是，Roghair等[20]未考慮操作壓力對氣泡群曳力系數(shù)的影響，因此，本文在式（9）的基礎(chǔ)上，根據(jù)冷態(tài)實(shí)驗(yàn)[21]測量的全塔平均氣含率（壓差法測量）和氣泡直徑（電導(dǎo)探針測量），調(diào)整單氣泡曳力模型中的曳力系數(shù)，使平均氣含率的模擬結(jié)果符合實(shí)驗(yàn)測量值。根據(jù)調(diào)定后的曳力系數(shù)、氣含率、氣泡直徑及新添加的量綱1壓力項(xiàng)，擬合模型參數(shù)，得到適用于加壓條件下改進(jìn)的氣泡群曳力系數(shù)計(jì)算公式

其中，CD,0表示單個(gè)氣泡的曳力系數(shù)，本文采用Schiller & Naumann模型計(jì)算CD,0，湍動鼓泡塔中氣泡的相對Reynolds數(shù)Re通常均大于1000，CD,0可取常數(shù)0.44。

1.2.2 升力 鼓泡塔內(nèi)氣泡浮升運(yùn)動過程中，由于液相在其運(yùn)動方向的兩側(cè)流動不對稱，會造成氣泡兩側(cè)壓力不平衡，從而產(chǎn)生垂直于氣泡運(yùn)動方向的升力。大氣泡在鼓泡塔內(nèi)受到液相剪切流動、曳力以及渦旋的作用后，易發(fā)生變形，并且在氣泡的后方產(chǎn)生偏斜尾渦，推動氣泡橫向運(yùn)動。Drew等[22]給出了單位體積內(nèi)，氣泡群所受合升力的計(jì)算公式

在鼓泡塔充分發(fā)展段，如果氣泡升力系數(shù) CL所取的符號不同，氣泡升力會指向塔中心或者塔壁，使氣含率徑向分布對應(yīng)形成中心高的拋物線分布，或者形成壁面峰的鞍形分布，此外，升力系數(shù) CL取值大小會顯著影響氣含率徑向分布的梯度。

1.2.3 湍流擴(kuò)散力 對于湍動鼓泡塔流場，湍流擴(kuò)散力體現(xiàn)了氣液兩相隨機(jī)運(yùn)動所產(chǎn)生的擴(kuò)散作用，在該力的作用下氣含率徑向分布趨于均勻。Lopez de Bertodano[23]將氣泡的隨機(jī)運(yùn)動類比于分子熱運(yùn)動，推導(dǎo)出單位體積氣泡受到液相渦旋所受的擴(kuò)散力表達(dá)式

考慮到湍流擴(kuò)散力主要在氣含率較高的區(qū)域起作用，因此，引入極限函數(shù) fTD,limiting對湍流擴(kuò)散力作用區(qū)域進(jìn)行調(diào)整，湍流擴(kuò)散力可表示為

其中，εG,1和εG,2為湍流擴(kuò)散力起作用的局部氣含率下限及上限值，本文所述模型中將其分別取值為0.3和0.7。

文獻(xiàn)[7-8]中認(rèn)為鼓泡塔內(nèi)徑向力作用力對氣含率的分布起著主導(dǎo)作用，當(dāng)升力和湍流擴(kuò)散力達(dá)到動態(tài)平衡，氣含率的分布也隨之趨于穩(wěn)定。張煜等[7]通過常壓鼓泡塔冷態(tài)實(shí)驗(yàn)數(shù)據(jù)和一維模型計(jì)算結(jié)果，確定了氣泡升力系數(shù) CL和湍流擴(kuò)散力系數(shù)CTD符合如下函數(shù)關(guān)系

式（15）中的系數(shù)0.2是一個(gè)固定的值，不隨實(shí)驗(yàn)條件而改變。本模型假設(shè)該徑向力平衡機(jī)制同樣適用于加壓湍動鼓泡塔流場模擬，湍流擴(kuò)散力系數(shù) CTD取 1.0，液含率與升力系數(shù)對應(yīng)情況如表 1所示。

表1 不同液含率下的升力系數(shù)Table 1 Lateral lift coefficient under different liquid holdup

1.2.4 壁面潤滑力 當(dāng)浮升氣泡接近壁面時(shí)，氣泡表面液體速度分布受到壁效應(yīng)的影響，表現(xiàn)為氣泡與壁面之間的液相排出速率下降，而在另一側(cè)的排液速率增加，從而產(chǎn)生將氣泡推離壁面的力，稱之為壁面潤滑力[24]。本文模型采用 Tomiyama[25]的壁面潤滑力模型，單位體積的氣泡所受的壁面潤滑力可表述為

式中，壁面潤滑力系數(shù)CWL可表示為

式中，yw為壁面距離，Cw可表示為氣泡Eo數(shù)的分段函數(shù)

在進(jìn)行加壓鼓泡塔軸對稱二維動態(tài)模擬時(shí)發(fā)現(xiàn)，當(dāng)網(wǎng)格密度較高時(shí)，在充分發(fā)展段經(jīng)常出現(xiàn)壁面氣含率高中心低的情況，這與實(shí)驗(yàn)測量結(jié)果不符。當(dāng)模型中考慮壁面潤滑力后，能夠有效避免出現(xiàn)氣含率分布的非物理解。

2 模型設(shè)置及求解方法

本文參照文獻(xiàn)[21]所述內(nèi)徑0.3 m冷模加壓鼓泡塔裝置生成二維計(jì)算網(wǎng)格，計(jì)算區(qū)域?qū)挾葹?.15 m、高度為6.6 m。根據(jù)表觀氣速和操作壓力設(shè)置初始靜液位高度H0，防止膨脹床層高度超過計(jì)算域高度，造成液體溢出，靜液位設(shè)定高度為2.5～3.0 m。設(shè)定H0以下的液含率為1.0，H0以上的液含率為 0，氣液相均處于靜止?fàn)顟B(tài)。邊界條件為：鼓泡塔底部的孔板分布器（r/R＜0.8）簡化為氣相速度進(jìn)口邊界條件；塔頂氣相出口設(shè)置為壓力出口；中心軸采用軸對稱邊界條件；壁面的氣液相均設(shè)定為無滑移條件，湍流方程設(shè)定為標(biāo)準(zhǔn)壁面函數(shù)，如圖 1所示。本文采用文獻(xiàn)[21]的實(shí)驗(yàn)數(shù)據(jù)作為檢驗(yàn)?zāi)Ｐ偷囊罁?jù)，其中局部氣含率、氣泡上升速度及氣泡Suater直徑，是采用雙電導(dǎo)探針，根據(jù)氣液兩相電導(dǎo)率差異較大的原理來測量的，對氣泡上升速度的測量是通過考察電導(dǎo)針處，一段時(shí)間時(shí)間內(nèi)統(tǒng)計(jì)所有氣泡上升速度的數(shù)量平均獲得的。

通過用戶自定義方程UDF，將前文所述的各氣泡界面力模型嵌入FLUENT計(jì)算程序，全塔設(shè)定單一氣泡直徑等于秦玉建[21]的冷態(tài)實(shí)驗(yàn)測量值。模型采用 FLUENT 15.0進(jìn)行動態(tài)模擬，采用Phase-Coupled SIMPLE算法求解壓力-速度耦合，離散空間梯度使用Least Squares Cell based方法，導(dǎo)數(shù)采用一階迎風(fēng)格式，時(shí)間導(dǎo)數(shù)采用一階隱式格式，為防止發(fā)散，初始時(shí)間步長為0.001 s，待流場充分發(fā)展后，可以逐漸增加步長至0.004 s。模擬結(jié)果表明，當(dāng)計(jì)算物理時(shí)間達(dá)到20 s后，加壓鼓泡塔充分發(fā)展段的徑向氣含率分布基本不變，可近似認(rèn)為流場充分發(fā)展，以表觀氣速0.31 m·s?1，操作壓力為0.5 MPa，距氣體分布器2.4 m處氣含率徑向分布為例，如圖2所示。本文取模擬時(shí)間50 s后的流場進(jìn)行時(shí)均，圖3給出了不同采樣時(shí)間對應(yīng)的時(shí)均氣含率分布，可以看出當(dāng)流場基本穩(wěn)定后，采樣時(shí)間大于25 s后的氣含率分布基本一致。因此，文中所述模擬均提取流場發(fā)展50 s后的數(shù)據(jù)，再時(shí)均50 s作為數(shù)值模擬結(jié)果。

圖4給出了網(wǎng)格數(shù)量對充分發(fā)展段氣含率徑向分布模擬結(jié)果的影響，考察了網(wǎng)格數(shù)量為3960、5940、17160和47520個(gè)4種不同的網(wǎng)格，并且均在近壁區(qū)域加密了4層邊界層網(wǎng)格。從圖中可以看出，加密網(wǎng)格對中心氣含率分布影響較小，主要使近壁面氣含率分布變得更為平滑。綜合考慮計(jì)算代價(jià)和數(shù)值模擬精度，采用網(wǎng)格數(shù)為5940個(gè)。

圖1 加壓鼓泡塔結(jié)構(gòu)及邊界條件Fig.1 Schematic diagram of pressurized bubble column and boundary conditions

圖2 模擬的氣含率徑向分布隨時(shí)間的變化Fig. 2 Simulated radial gas holdup profiles varied with time

圖3 模擬時(shí)均氣含率分布隨采樣時(shí)間的變化Fig. 3 Simulated radial gas holdup profiles varied with sampling time

圖4 網(wǎng)格密度對氣含率徑向分布的影響Fig. 4 Influence of mesh density on radial gas holdup distribution (VG=0.27 m·s?1，P=0.5 MPa，H=2.6 m)

3 結(jié)果與討論

3.1 流場分布

鼓泡塔流場的最主要特征是氣含率在徑向上的不均勻分布，即中心區(qū)域氣含率高，邊壁處氣含率較低。這種氣含率的不均勻分布形成密度差，從而驅(qū)動液相流動，形成液相塔內(nèi)的大尺寸循環(huán)流動，呈現(xiàn)出中心區(qū)液相向上運(yùn)動，近壁區(qū)液體向下回流的規(guī)律。圖5是CFD模擬得到的，表觀氣速為0.27 m·s?1，操作壓力為1.5 MPa，內(nèi)徑0.3 m加壓鼓泡塔內(nèi)時(shí)均氣含率二維分布。模擬結(jié)果表明，在鼓泡塔底部，回流液體在壁面附近富集，使氣泡向塔中心集中。經(jīng)過約1倍塔徑發(fā)展后，氣含率徑向分布呈現(xiàn)出中心高、邊壁低的趨勢，并且不隨高度變化，形成充分充分發(fā)展段，直至接近塔上部液位處。在氣體分布器上方2.6 m和3.0 m處，冷態(tài)測量得到的氣含率徑向分布基本相同，可認(rèn)為處于充分發(fā)展段，模擬的氣含率分布與實(shí)驗(yàn)結(jié)果符合較好，這表明，本文所提出的基于氣泡群界面力的加壓鼓泡塔二維模型，能夠確定氣含率分布。

圖5 氣含率二維分布模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的比較Fig.5 Comparison of simulated 2-D gas holdup distribution and experimental data(VG=0.27 m·s?1，P=1.5 MPa)

3.2 曳力模型比較分析

圖6（a）為操作壓力為0.5 MPa時(shí)，不同的曳力模型計(jì)算的充分發(fā)展段氣含率徑向分布與實(shí)驗(yàn)結(jié)果的比較，與Schiller & Naumann模型相比，改進(jìn)的氣泡群曳力模型和Morsi & Alexander曳力模型模擬結(jié)果更接近實(shí)驗(yàn)值。Schiller & Naumann模型所獲得模擬結(jié)果更接近實(shí)驗(yàn)值；而當(dāng)操作壓力增至1.5 MPa時(shí)，Schiller & Naumann模型和Morsi & Alexander模型的氣含率分布與低壓時(shí)的模擬結(jié)果基本相同，不能反映壓力增加后局部氣含率增大的趨勢。然而，改進(jìn)的氣泡群曳力模型能夠描述壓力對局部氣含率的影響，并且模擬結(jié)果與實(shí)測值符合較好，如圖6（b）所示。

圖7為不同操作壓力下，上述3種曳力模型計(jì)算的軸向氣速徑向分布與實(shí)驗(yàn)數(shù)據(jù)的比較。冷態(tài)實(shí)驗(yàn)數(shù)據(jù)顯示，當(dāng)表觀氣速為0.27 m·s?1時(shí)，隨著操作壓力的升高，局部氣相速度略有下降。從圖7中可以看出，3種曳力模型對塔中心區(qū)域（r/R＜0.3）的氣相速度較為接近，越靠近塔壁模擬計(jì)算值與實(shí)驗(yàn)值差別越大。當(dāng)壓力為0.5 MPa時(shí)，改進(jìn)的氣泡群曳力模型的計(jì)算結(jié)果比另兩種模型更接近實(shí)驗(yàn)值。而當(dāng)壓力增加到1.5 MPa時(shí)，Schiller & Naumann模型和Morsi & Alexander模型計(jì)算得到氣相速度與低壓下的計(jì)算值相比，并沒有如實(shí)驗(yàn)數(shù)據(jù)表現(xiàn)出下降趨勢，反而略有增加，然而，改進(jìn)的氣泡群壓力模型能夠定性反映軸向氣速隨操作壓力增大而減小的趨勢。由于本文給出的氣泡群壓力系數(shù)表達(dá)式是基于平均氣含率和氣泡直徑的冷態(tài)數(shù)據(jù)回歸得到的，尚且沒有考慮對氣相上升速度的修正，這可能是導(dǎo)致中心區(qū)域以外的軸向氣速模擬結(jié)果偏差較大的一個(gè)原因。在冷態(tài)實(shí)驗(yàn)中，氣相是以氣泡離散相存在的，測量的氣相速度即是氣泡的上升速度，而模擬出來的氣相上升速度是將氣相作為連續(xù)相考慮的，即作為連續(xù)介質(zhì)的氣相與液相充分混合，因此所計(jì)算出來的氣相上升速度相較于實(shí)驗(yàn)測量氣泡上升速度值偏小，這可能也是氣相上升速度存在偏差的另一個(gè)原因。

圖6 不同曳力模型對氣含率徑向分布的影響Fig.6 Effect of various drag model on radial gas holdup distribution(VG=0.27 m·s?1, H=2.6 m)

圖7 不同曳力模型對氣相上升速度徑向分布的影響Fig.7 Effect of various drag model on gas rising velocity radial distribution(VG=0.27 m·s?1,H=2.6 m)

綜上，改進(jìn)的氣泡群曳力模型能夠用于不同操作壓力下的鼓泡塔流場，盡管對氣相上升速度分布模擬精度不夠理想，但是定性能夠體現(xiàn)壓力對軸向氣相速度的影響，并且高壓條件下模擬得到的氣含率徑向分布與實(shí)驗(yàn)值符合較好。因此，該模型與Schiller & Naumann和Morsi & Alexander模型相比模擬準(zhǔn)確度有所提高，基本能夠描述湍動加壓鼓泡塔內(nèi)，操作壓力對充分發(fā)展段氣含率和軸向氣速分布的影響趨勢。

3.3 表觀氣速對流場的影響

表觀氣速是影響加壓鼓泡塔內(nèi)氣液兩相流動的關(guān)鍵參數(shù)，冷態(tài)實(shí)驗(yàn)結(jié)果表明，全塔平均氣含率隨表觀氣速提高而增大。表2列出了操作壓力為0.5 MPa，表觀氣速為0.20～0.31 m·s?1時(shí)，全塔平均氣含率模擬值與實(shí)驗(yàn)值的比較結(jié)果，可以看出模型對平均氣含率的預(yù)報(bào)較為準(zhǔn)確，最大誤差絕對值不超過 2%，能夠定量描述加壓鼓泡塔平均氣含率隨表觀氣速變化的規(guī)律。

圖8為上述表觀氣速下，模擬得到的局部氣含率與實(shí)驗(yàn)結(jié)果的比較。冷態(tài)實(shí)驗(yàn)數(shù)據(jù)顯示，隨著表觀氣速的增加，局部氣含率提高，氣含率沿徑向分布的陡峭程度基本不變。從圖8中可以看出，模擬得到的氣含率分布雖然略高于實(shí)驗(yàn)測量結(jié)果，尤其是表觀氣速為0.20 m·s?1時(shí)偏差相對較大，但是模擬的氣含率分布趨勢與實(shí)驗(yàn)結(jié)果較為接近。模擬和實(shí)驗(yàn)結(jié)果定量上的差異，可主要?dú)w結(jié)于兩個(gè)方面：一方面，改進(jìn)的氣泡群曳力系數(shù)關(guān)聯(lián)式[式（10）]中的參數(shù)是基于壓差法測量得到的平均氣含率的數(shù)據(jù)回歸而來；另一方面，由于氣泡在浮升過程中受到液相湍流渦旋的擾動，會發(fā)生一定隨機(jī)擺動，造成氣泡與電導(dǎo)探針接觸點(diǎn)的氣泡表面法向與探針角度過大，氣泡易從探針尖端脫離，使得探針測量的局部氣含率偏小。秦玉建[21]的冷態(tài)實(shí)驗(yàn)數(shù)據(jù)顯示，通過電導(dǎo)針測量的局部氣含率積分得到的截面平均氣含率，通常略低于差壓法的測量值，尤其當(dāng)表觀氣速為0.20 m·s?1時(shí)，兩者相差較為明顯。

表2 不同表觀氣速下平均氣含率模擬值與實(shí)驗(yàn)值比較Table 2 Comparison between simulation results and experiment values for different superficial gas velocity (P=0.5 MPa)

圖9為不同表觀氣速下，軸向氣相速度模擬與實(shí)驗(yàn)結(jié)果的比較，從圖中可以看出，隨著表觀氣速由0.20 m·s?1增加到0.31 m·s?1，實(shí)驗(yàn)測量的氣相速度略有提高，而徑向分布趨勢基本相同；而模擬得到的氣相速度隨著表觀氣速增加而增大，且分布變得更為陡峭，與實(shí)測結(jié)果相差較多。

圖9 表觀氣速對軸向氣速分布的影響Fig.9 Effect of superficial gas velocity on axial gas rising velocity along radius (P=0.5 MPa, H=2.6 m)

3.4 操作壓力對流場的影響

操作壓力對流場的影響顯著，實(shí)驗(yàn)表明隨著操作壓力的增加，氣泡尺寸略有減小，氣泡上升速率降低，使得氣泡在流場內(nèi)停留時(shí)間增加，使氣含率明顯提高。表3給出了在表觀氣速為0.27 m·s?1，壓力為0.5 MPa～2.0 MPa下，平均氣含率模擬結(jié)果與實(shí)驗(yàn)值的比較，模擬結(jié)果略低于實(shí)測值，最大相對誤差的絕對值為3.3%。

表3 不同操作壓力下平均氣含率模擬值與實(shí)驗(yàn)值比較Table 3 Comparison between simulation results and experimental values under various operating pressure(VG=0.27 m·s?1)

從圖10中可以看出，壓力為0.5～1.5 MPa時(shí)，模擬的氣含率分布略高于實(shí)驗(yàn)值，而當(dāng)壓力達(dá)到2.0 MPa時(shí)，塔中心處實(shí)驗(yàn)結(jié)果高于模擬值。數(shù)值模擬與實(shí)測的氣含率分布曲線變化趨勢基本一致，隨著操作壓力的增加兩者的氣含率均逐漸增加。但是，在部分操作壓力下，模擬的局部氣含率與實(shí)驗(yàn)結(jié)果相比還存在一定偏差。

圖 11表明模擬的氣相速度隨著壓力升高而下降，與實(shí)驗(yàn)結(jié)果所呈現(xiàn)的趨勢基本一致。然而，模擬的軸向氣相速度分布比實(shí)測分布更為陡峭，并且模擬結(jié)果在總體上也略低于實(shí)驗(yàn)測量結(jié)果。由于模型假設(shè)全塔氣泡尺寸不變，而實(shí)際上湍動鼓泡塔內(nèi)氣泡直徑存在一定空間分布，而氣泡運(yùn)動速度通常與其尺寸相關(guān)，使得模型雖然能夠定性描述軸向氣速分布隨操作壓力變化的基本趨勢，但是在定量上還存在明顯偏差。在后續(xù)工作中將引入數(shù)量衡算方程（population balance equation），模擬得到加壓鼓泡塔內(nèi)氣泡尺度分布，并以此校正氣相速度分布，提高模型的預(yù)測精度。

圖10 壓力對氣含率徑向分布的影響Fig.10 Effect of operating pressure on radial gas holdup distribution（VG=0.27 m·s?1, H=2.6 m）

圖11 壓力對軸向氣速分布的影響Fig.11 Effect of operating pressure on gas rising velocity radial distribution（VG=0.27 m·s?1）

4 結(jié) 論

（1）本文采用雙流體模型對內(nèi)徑0.3 m加壓湍動鼓泡塔內(nèi)的氣液兩相流動，進(jìn)行了動態(tài)二維軸對稱數(shù)值模擬。針對加壓條件下，采用改進(jìn)的氣泡群曳力系數(shù)關(guān)聯(lián)式，同時(shí)考慮氣泡群徑向力平衡機(jī)制，以及壁面潤滑力建立的加壓鼓泡塔流體力學(xué)模型，數(shù)值計(jì)算的穩(wěn)定性和效率較高，模擬得到的氣含率分布與實(shí)驗(yàn)結(jié)果符合較好，能夠描述加壓鼓泡塔內(nèi)氣含率和軸向氣速分布隨操作條件變化的趨勢。

（2）與其他采用單氣泡曳力系數(shù)的模型相比，改進(jìn)的氣泡群曳力模型能夠較為準(zhǔn)確地模擬加壓鼓泡塔總體氣含率和氣含率分布，正確描述軸向氣速分布隨表觀氣速和壓力的變化趨勢。對處于加壓湍動鼓泡狀態(tài)的工業(yè)漿態(tài)床反應(yīng)器內(nèi)流場特性參數(shù)分布的預(yù)測和變化規(guī)律的認(rèn)識具有重要意義。

符 號 說 明

CD——曳力系數(shù)

CD,0——單氣泡曳力系數(shù)

CWL——壁面潤滑力系數(shù)

CTD——湍流擴(kuò)散力系數(shù)

D——塔徑，m

dB——?dú)馀葜睆?，m

Fex——相間作用力，N·m?3

FL——?dú)馀萆?，N·m?3

FTD——?dú)馀萃牧鲾U(kuò)散力，N·m?3

FWL——?dú)馀荼诿鏉櫥Γ琋·m?3

g——重力加速度，m·s?2

H——軸向高度，m

k, kL——湍流動能, m2·s?2

Nc——計(jì)算網(wǎng)格數(shù)量

nW——指向壁面外的單位法向量

P——壓力，Pa

P0——標(biāo)況下大氣壓，Pa

R——塔半徑，m

r——徑向位置，m

u——?dú)庀嗷蛞合嗨俣?，m·s?1

uG——?dú)庀嗨俣?，m·s?1

uL——液相速度，m·s?1

VG——表觀氣速，m·s?1

yw——壁面距離，m

ε——湍流耗散率，m2·s?3

εG——局部氣含率

εL——局部液含率

ε?G——平均氣含率

ε?L——平均液含率

μt——湍流黏度，Pa·s

ρG——?dú)怏w密度，kg·m?3

ρL——液體密度，kg·m?3

σ——液體表面張力，N·m?3

下角標(biāo)

G——?dú)庀?/p>

L——液相

i——液相(1)或氣相(2)

[1]許世峰, 王斯民, 李彩霞, 等. FT合成漿態(tài)床反應(yīng)器的研究進(jìn)展[J]. 化工進(jìn)展, 2013, 32(S1): 4-8.XU S F, WAMG S M, LI C X, et al. Recent advances of the slurry-bed reactor FT synthesis[J]. Chemical Industry and Engineering Progress, 2013, 32(S1): 4-8.

[2]KRISHNA R,VAN BATEN J M. Eulerian simulations of bubble columns operating at elevated pressures in the churn turbulent flow regime[J]. Chemical Engineering Science, 2001, 56(21): 6249-6258.

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

[4]JOSHI J B. Computational flow modeling and design of bubble column reactors[J]. Chemical Engineering Science, 2001, 56: 5893-5933.

[5]TABID 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, 2007, 139(3): 589-614.

[6]LUCAS D, KREPPER E, PRASSER H M. Prediction of radial gas profiles in vertical pipe flow on the basis of buble size distribution[J]. International Journal of Thermal Sciences, 2001, 40(3): 217-226.

[7]張煜, 李紅波, 李兆奇, 等. 湍動漿態(tài)床流體力學(xué)研究(Ⅳ): 帶垂直列管束的漿態(tài)床流體力學(xué)模型與模擬[J]. 化工學(xué)報(bào), 2011, 62(12): 3373-3380. ZHANG Y, LI H B, LI Z Q, et al. Hydrodynamics of turbulent slurry bubble column (Ⅳ): Modeling and simulation of bubble column with internal pipe bundles[J]. CIESC Journal, 2011, 62(12): 3373-3380.

[8]李兆奇, 王麗軍, 管小平, 等. 基于徑向力平衡的鼓泡塔二維流體力學(xué)模型[J]. 化工學(xué)報(bào), 2014, 65(11):4221-4230. LI Z Q, WANG L J, GUAN X P, et al. 2-D hydrodynamic model of bubble column based on lateral-force balance[J]. CIESC Journal, 2014, 65(11): 4221-4230.

[9]DREW D A. Mathematical modeling of two-phase fow[J]. Annual Review of Fluid Mechanics, 1983, 15(1): 261-291.

[10]RAFIQUE M, CHEN P, DUDUKOVI? M P. Computational modeling of gas-liquid flow in bubble column[J]. Reviews in Chemical Engineering, 2004, 20(3): 225-375.

[11]DREW D A, PASSMAN S L. Applied Mathematical Sciences[M]. New York: Springer, 1999.

[12]PFLEGER D, BECKER S. Modeling and simulation of the dynamic flow behaviour in a bubble column[J]. Chemical Engineering Science, 2001, 56(4): 1737-1747.

[13]MUDDE R F, SIMONIN O. Two and three dimensional simulation of a bubble plume using a two fluid model[J]. Chemical Engineering Science, 1999, 54(21): 5061-5069.

[14]BUWA V V, RANADE V V. Dynamics of gas iquid flows in rectangular bubble columns: experiments and single/multi groups imulations[J]. Chemical Engineering Science, 2002, 57(22): 4715-4736.

[15]PFLEGER D, GOMES S, GILBERT N, et al. Hydrodynamic simulation of laboratory scale bubble columns, fundamental studies of Eulerian-Eulerian modeling approach[J]. Chemical Engineering Science, 1999, 54(21): 5091-5099.

[16]李光, 楊曉鋼, 戴干策. 鼓泡塔反應(yīng)器氣液兩相流 CFD數(shù)值模擬[J]. 化工學(xué)報(bào), 2008, 59(8): 1954-1965. LI G, YANG X G, DAI G C. CFD simulation of gas-liquid flow in bubble column[J]. Journal of Chemical Industry and Engineering (China), 2008, 59(8): 1954-1965.

[17]李光. 鼓泡塔反應(yīng)器氣液兩相流數(shù)值模擬模型及應(yīng)用[D]. 上海:華東理工大學(xué), 2010. LI G. Two-phase flow dynamical simulations and modelling of bubble column reactors[D]. Shanghai: East China University of Science and Technology, 2010.

[18]KHOPKAR A R, RAMMOHAN A R, RANADE V V, et al. Gas-liquid flow generated by a Rushton turbine in stirred vessel: CARPT/CT measurements and CFD simulation[J]. Chemical Engineering Science, 2005, 60: 2215-2229.

[19]LAU Y M, ROGHAIR I, DEEN N G, et al. Numerical investigation of the drag closure for bubbles in bubble swarms[J]. Chemical Engineering Science, 2011, 66(14): 3309-3316.

[20]ROGHAIR I VAN SINT ANNALAND M, KUIPERS H J A M. Drag force and clustering in bubble swarms [J]. AIChE Journal, 2013, 59(5): 1791-1800.

[21]秦玉建. 加壓氣液鼓泡塔流體參數(shù)的測量與CFD數(shù)值模擬[D]. 北京: 北京化工大學(xué), 2012. QIN Y J. The measurement of hydrodynamic parameters and CFD simulation of gas-liquid flow behaviors in a pressurized bubble column[D]. Beijing: Beijing University of Chemical Technology, 2012.

[22]DREW D A, LAHEY R T. In Particulate Two-phase Flow[M]. Boston: Butterworth Geinemann, 1993: 509-566.

[23]LOPEZ DE BERTODANO M. Turbulent bubbly two-phase flow in a triangular duct [D]. New York: Rensselaer Rolytechnic Institute, 1992.

[24]ANTAL S P, LAHEY JR R T, FLAHERTY J E. Analysis of phase distribution in fully developed laminar bubbly two-phase flow[J]. Int. J. Multiphase Flow, 1991, 17(5): 635-652.

[25]TOMIYAMA A. Struggle with computational bubble dynamics[C]// ICMF'98, 3rd International Conference on Multiphase Flow. Lyon, France, 1998: 1-18.

Hydrodynamics simulation of pressurized bubble column based on bubble swarm interphase force models

LIU Xin1,2, ZHANG Yu2, ZHANG Li2, JIN Haibo1
(1School of Chemical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China;2Research and Development Center, Synfuels China, Beijing 101407, China)

The commercial bubble column reactors are normally operated under the conditions of pressurized state and churn turbulent bubbly flow which the superficial gas velocity always exceeded 0.1 m·s?1. However, most of the published literatures are focus on the cold model experimental investigation of the flow field characteristics in pressurized churn-turbulent bubble columns, and the empirical or semi-empirical correlations which were on the basis of the cold model experimental data to quantitatively estimate the characteristic parameters. Therefore, it is lack of the systematic study on the flow characteristics in the pressurized turbulent bubble columns by the hydrodynamic models with interactions between bubbles and liquid phase. In this work, a steady 2-D axisymmetric hydrodynamic two-fluid model of pressurized bubble column was developed on the platform of FLUENT 15.0 with the UDF sub-models which described the bubble swarm interface forces. The numerical investigations upon the flow field characteristics of pressurized churn turbulent bubble column with the operation pressure varying from 0.5 MPa to 2.0 MPa, and superficial gas velocity from 0.2 m·s?1to 0.31 m·s?1were performed. Moreover, the comparisons between the simulation results and cold model experimental data were implemented. The comparisons indicated that, the model with the bubble swarm interface forces including improved bubble swarm drag force model, lateral force balance model, and wall lubrication force model, is able to reflect the influences of operation pressure and superficial gas velocity upon flow pattern, that is, the elevatingoperation pressure leads to increase of time averaged gas holdup and decrease of time averaged axial gas phase velocity, whereas for the ascending superficial gas velocity, the gas holdup and gas phase velocity are both raising. Consequently this mathematical model could correctly predict the basic gas-liquid flow pattern in the turbulent pressurized bubble columns, and simulation results are in good agreement with the experimental data.

bubble column; gas-liquid flow; computational fluid dynamics; bubble swarm drag force model; lateral- force balance model; wall lubrication force model

Prof. JIN Haibo, jinhaibo@bipt.edu.cn

TQ 021.1

：A

：0438—1157（2017）01—0087—10

10.11949/j.issn.0438-1157.20160762

2016-06-01收到初稿，2016-10-20收到修改稿。

聯(lián)系人：靳海波。

：劉鑫（1989—），男，碩士研究生。

國家自然科學(xué)基金重大研究計(jì)劃項(xiàng)目（91634101）；北京市屬高等學(xué)校高層次人才引進(jìn)與培養(yǎng)計(jì)劃項(xiàng)目（CIT&TCD20130325）。

Received date: 2016-06-01.

Foundation item: supported by the National Natural Science Foundation of China (91634101) and the Importation and Development of High-Caliber Talents Project of Beijing Municipal Institutions(CIT&TCD20130325).