三維單旋導(dǎo)向擋板鼓泡流化床內(nèi)氣固流動的CPFD 模擬

曲 悅, 潘 騰, 楊越林, 成有為, 王麗軍, 李 希

(1. 浙江大學(xué) 化學(xué)工程與生物工程學(xué)院, 浙江 杭州310027;2. 浙江大學(xué)衢州研究院, 浙江 衢州 324000)

1 前 言

氣固鼓泡流化床是化學(xué)工業(yè)中常見的流化床型式，例如催化裂化和甲醇制烯烴裝置中的汽提器和部分立式輸送管路均在鼓泡流態(tài)化下操作[1-2]。氣固鼓泡流化床中存在乳相和氣泡相，乳相中的氣體由顆粒之間的空隙通過顆粒床層，氣泡相中的氣體則以氣泡的形式通過顆粒床層，因此，氣泡的生成、發(fā)展和運動是影響鼓泡流化床反應(yīng)器特性的關(guān)鍵因素[3]。大尺寸氣泡的存在可能導(dǎo)致部分氣相無法有效地與固相接觸、加劇顆粒的返混和夾帶等問題，從而使反應(yīng)器內(nèi)的化學(xué)反應(yīng)難以達到預(yù)期的轉(zhuǎn)化率和選擇性[4]。為了提高鼓泡流化床的流態(tài)化性能，可以通過設(shè)置內(nèi)構(gòu)件來抑制氣泡生長、破碎大尺寸氣泡、強化氣固混合等[5]。內(nèi)構(gòu)件技術(shù)作為一個簡單、經(jīng)濟的改善流態(tài)化性能的手段，目前已成功應(yīng)用于丁烯氧化脫氫、苯氧化法制順酐等工業(yè)過程中[3]。

為了研究內(nèi)構(gòu)件對鼓泡流化床內(nèi)氣固流動的影響，許多學(xué)者已成功借助計算流體力學(xué)(CFD)方法獲得了實驗難以測量的微觀流動狀況，例如氣泡的尺寸、形狀和頻率及氣固混合特性等[6-11]。石孝剛等[6]使用雙流體模型 (TFM) 研究帶有傾斜擋板的二維鼓泡流化床內(nèi)的氣固流動特性，結(jié)果表明，擋板影響區(qū)域內(nèi)氣泡的平均尺寸降低、氣泡數(shù)量增加；且氣速越高，擋板的作用越強。Gao 等[7]使用TFM 模型計算對比了帶有V 型擋板和不帶內(nèi)構(gòu)件的FCC 汽提器中的氣固混合特性，結(jié)果表明，V 型擋板汽提器可以有效提高氣固混合效率、降低催化劑返混，從而提高汽提器的汽提效率。Zhu 等[10]利用計算顆粒流體力學(xué)(CPFD)方法對化學(xué)鏈燃燒過程中的鼓泡流化床進行模擬，探究了擋板的數(shù)量、開孔率等參數(shù)對鼓泡流化床中氣泡尺寸、壓力梯度、顆粒分布等方面的影響，結(jié)果表明，擋板顯著改善了流化床中的氣固接觸。Yang 等[11]使用CPFD 方法對二維單旋導(dǎo)向擋板鼓泡流化床進行模擬，探究單旋導(dǎo)向擋板對床內(nèi)固含率、壓差波動、氣泡行為的影響，提出了二維單旋擋板在較高表觀氣速下破碎氣泡的機制。

CPFD 是一種結(jié)合歐拉雙流體模型與拉格朗日離散模型各自優(yōu)點的計算流體力學(xué)方法，通過引入“計算顆?！边@一概念，將多個真實顆粒打包進行計算，在保證計算精度的同時極大地降低了計算資源和成本。目前已有許多學(xué)者采用CPFD 方法研究無內(nèi)構(gòu)件鼓泡流化床內(nèi)的氣固流動特性[12-15]，但針對三維單旋擋板鼓泡流化床的研究相對較少。本研究使用CPFD 方法耦合Igci 曳力模型，對三維單旋導(dǎo)向擋板鼓泡流化床內(nèi)的氣固流動進行模擬，首先將軸向和徑向時均固含率分布模擬結(jié)果與實驗數(shù)據(jù)進行對比來驗證計算方法的可靠性，隨后考察有、無內(nèi)構(gòu)件時鼓泡流化床內(nèi)的氣泡平均當(dāng)量直徑分布和氣體停留時間分布。

2 計算方法

2.1 C PFD 數(shù)學(xué)模型

使用CPFD Barracuda 軟件對單旋導(dǎo)向擋板鼓泡流化床進行模擬。在CPFD 中，氣相被視為連續(xù)相，通過Naiver-Stokes 方程來描述，顆粒被視為離散相，通過拉格朗日框架下的MP-PIC 方法描述。引入“計算顆?！边@一概念，將多個具有同樣物化性質(zhì)和運動狀態(tài)的真實顆粒打包，可大幅度降低計算成本。關(guān)于CPFD 方法的詳細描述及相關(guān)表達式見文獻[16]。

2.2 氣固曳力模型

氣固曳力是氣固兩相流中最重要的作用力之一，因此，選擇合適的曳力模型是準(zhǔn)確模擬和描述鼓泡流化床內(nèi)兩相流動的關(guān)鍵。通常有兩種類型曳力模型被應(yīng)用于模擬氣固兩相流：均勻模型 (homogeneous drag model) 和非均勻模型 (heterogeneous drag model)。前者以Wen&Yu 模型[17]和Gidaspow 模型[18]為代表，在應(yīng)用中往往會高估曳力從而導(dǎo)致床層過度膨脹。為了克服均勻曳力模型的缺點，多種非均勻曳力模型被提出，例如經(jīng)典的過濾曳力模型和EMMS 曳力模型。過濾曳力模型[19]是通過細網(wǎng)格模擬捕捉小尺度的非均勻結(jié)構(gòu)，并在時空尺度上平均化，從而建立曳力系數(shù)與氣固滑移速度、局部固含率等參數(shù)的本構(gòu)關(guān)聯(lián)式。EMMS 曳力模型[20]采用多尺度分析方法并運用能量最小原理，建立氣固曳力與鼓泡流化床整體、局部參數(shù)之間的關(guān)系。

本研究采用Igci 曳力模型即Igci 等[21]基于顆粒動理學(xué)和雙流體模型的細網(wǎng)格模擬推導(dǎo)出的過濾曳力模型對單旋導(dǎo)向擋板鼓泡流化床進行模擬計算，并通過實驗數(shù)據(jù)對模擬結(jié)果進行驗證。以下是Igci 曳力模型的數(shù)學(xué)表達式。

在CPFD 中，單顆粒的曳力通過式(1)和(2)描述：

式中：F為單顆粒曳力，N；mp為顆粒質(zhì)量，kg；gu和pu分別為氣體和顆粒的速度，m·s-1；D代表曳力方程(drag function)；ρg和ρp分別為氣體和顆粒的密度，kg·m-3；dp為顆粒直徑，m；Cd為曳力系數(shù)，在Wen&Yu 模型中Cd通過式(3)描述：

式中：φg為氣相體積分?jǐn)?shù)，Re為雷諾數(shù)。在CPFD Barracuda 軟件中，Re通過式(4)描述：

Igci 曳力模型在Wen&Yu 模型基礎(chǔ)上通過非均勻系數(shù)HIgci進行修正，如式(5)和(6)所示：

式中：DIgci和DWen&Yu分別表示經(jīng)Igci 和Wen&Yu 曳力模型計算得到的式(2)中的曳力方程；、H2D和ξ均為模型參數(shù)，如式(7)～(10)所示：

式中：Δfilter為過濾尺寸，m；ut為顆粒終端速度，m·s-1；μg為氣體黏度，Pa·s；φp為固相體積分?jǐn)?shù)，亦稱為固含率。

2.3 模擬設(shè)置

模擬對象為文獻中的一套帶有兩層單旋導(dǎo)向擋板的三維冷模實驗裝置[9,22]，其結(jié)構(gòu)如圖1 所示。該裝置中流態(tài)化段主體部分為φ140 mm×1 000 mm 的有機玻璃圓筒，底部采用金屬燒結(jié)板作為氣體分布器，內(nèi)置兩層單旋導(dǎo)向擋板，其中葉片寬20 mm、葉片傾角45°、葉片間距14 mm、擋板整體厚度14 mm，圖1 中箭頭表示傾斜方向。

圖1 單旋導(dǎo)向擋板結(jié)構(gòu)圖[9]Fig.1 Configuration of the louver baffle[9]

實驗所采用的顆粒為玻璃微珠，顆粒密度為2 450 kg·m-3，堆積密度為1 395 kg·m-3，平均粒徑為53 μm。該玻璃微珠屬于A 類顆粒，其粒徑分布如圖2 所示。在模擬計算中采用此真實粒徑分布。模擬過程中所設(shè)置的參數(shù)見表1，流態(tài)化介質(zhì)選擇空氣，計算中空氣由底部均勻進入鼓泡流化床 (見圖3)。流化床底部為速度邊界條件，頂部為壓力邊界條件。在CPFD 中使用擋板組件 (Baffles)對兩層單旋導(dǎo)向擋板內(nèi)構(gòu)件進行建模，如圖4 所示。在三維模擬計算中，時間步長初始為5×10-4s，隨后由程序自動調(diào)整以保證收斂條件判斷數(shù)介于0.8 與1.5 之間?？傆嬎銜r間為70 s，取最后20 s 進行時均處理和分析，從而得到計算結(jié)果。

圖2 玻璃微珠累積粒度分布[9]Fig.2 Cumulative particle size distribution of glass beads[9]

表1 模擬參數(shù)Table 1 Simulation parameters

圖3 單旋導(dǎo)向擋板三維鼓泡流化床模擬示意圖Fig.3 Schematic diagram of the simulated 3D bubbling fluidized bed with louver baffles

圖4 CPFD 擋板設(shè)置示意圖Fig.4 Schematic diagram of baffle arrangement in CPFD

3 結(jié)果與討論

3.1 網(wǎng)格數(shù)和計算顆粒數(shù)無關(guān)性分析

為了考察網(wǎng)格分辨率和計算顆粒數(shù)對模擬結(jié)果的影響，采用三種不同精度的網(wǎng)格和計算顆粒數(shù)對兩層單旋導(dǎo)向擋板鼓泡流化床進行模擬，表2 列出了所使用的三種精度的網(wǎng)格數(shù)和計算顆粒數(shù)配置。圖5 為軸向截面時均固含率模擬結(jié)果與實驗數(shù)據(jù)對比圖，圖中H為分布器上方垂直高度，m。由圖5 可見上述三種精度的參數(shù)配置下的模擬結(jié)果區(qū)別不大。綜合考慮計算成本和精度，采用其中網(wǎng)格數(shù)為159 840、計算顆粒數(shù)為2 131 200 進行后續(xù)計算。

表2 網(wǎng)格和計算顆粒數(shù)Table 2 Parameters of grids and particle clouds

圖5 不同分辨率下軸向時均固含率模擬值與實驗值對比Fig.5 Comparison of simulated value and experimental data of time-averaged axial solid volume fraction at different resolutions

3.2 軸向和徑向時均固含率分布

采用Igci 曳力模型分別對有、無兩層單旋導(dǎo)向擋板鼓泡流化床進行模擬，圖6 為軸向分布模擬結(jié)果與實驗數(shù)據(jù)對比，圖7 為H分別為0.15、0.20 和0.25 m 處水平截面上徑向分布模擬結(jié)果與實驗數(shù)據(jù)對比，圖中r/R為歸一化半徑。由圖可見，由Igci 曳力模型得到的徑向分布模擬結(jié)果與實驗值相對接近，與Yang 等[8-9]用雙流體模型 (TFM) 耦合結(jié)構(gòu)曳力模型所得到的軸、徑向分布模擬結(jié)果趨勢基本相同。圖6 顯示，CPFD 模擬成功地捕捉到了擋板附近區(qū)域降低現(xiàn)象；圖7 顯示，CPFD 模擬基本上可用于描述徑向固含率分布情況。通過軸、徑向模擬結(jié)果與實驗數(shù)據(jù)的對比可知，使用CPFD 方法耦合Igci 曳力模型對兩層單旋導(dǎo)向擋板鼓泡流化床進行模擬是可行的。

圖6 有、無兩層單旋導(dǎo)向擋板時鼓泡流化床軸向時均固含率模擬值與實驗數(shù)據(jù)對比Fig.6 Comparison of simulated value and experimental data of time-averaged axial solid volume fraction of the bubbling fluidized beds with/without two louver baffles

圖7 有、無兩層單旋導(dǎo)向擋板鼓泡流化床徑向時均固含率模擬值與實驗數(shù)據(jù)對比Fig.7 Comparison of simulated value and experimental data of time-averaged radial solid volume fraction of the bubbling fluidized beds with/without two louver baffles

3.3 氣泡平均當(dāng)量直徑分布

使用Bakshi 等[23]提出的多相流三維檢測和追蹤算法 (MS3DATA) 對有、無單旋導(dǎo)向擋板時的鼓泡流化床內(nèi)的氣泡進行統(tǒng)計和分析。在模擬計算過程中，每0.01 s 記錄1 次床內(nèi)的瞬時流場狀態(tài)，取50～70 s 內(nèi)的2 000 幀數(shù)據(jù)作為樣本追蹤和統(tǒng)計氣泡，設(shè)置氣泡相的閾值為0.7 (氣相體積分?jǐn)?shù))。采用Darton[24]和Rowe 等[25]的關(guān)聯(lián)式來預(yù)測鼓泡流化床內(nèi)的氣泡當(dāng)量直徑，其表達式分別為式(11)和 (12)，使用式(13)[26]來計算最大氣泡直徑。

式中：Db為氣泡平均當(dāng)量直徑，m；ug為表觀氣速，m·s-1；umf為起始流化速度，m·s-1，實驗中所用的玻璃微珠的起始流化速度為 0.005 1m·s-1；h0為常數(shù)，取h0= 0.03 m[7]；g為重力加速度，m·s-2；Dbm為最大氣泡直徑，m；At為氣體入口截面積，m2。

圖8 為有、無兩層單旋導(dǎo)向擋板時鼓泡流化床內(nèi)的Db隨H統(tǒng)計分析結(jié)果。由圖可見，有單旋導(dǎo)向擋板時的床內(nèi)氣泡特征不同于無導(dǎo)向擋板即空筒鼓泡流化床時的氣泡特征。在空筒鼓泡流化床中，隨H增大Db不斷增大，且未出現(xiàn)氣泡破碎現(xiàn)象；而在有單旋導(dǎo)向擋板鼓泡流化床中，第一層擋板下方的Db大小以及變化規(guī)律與空筒鼓泡流化床相似，但經(jīng)過第一層擋板時，Db變小，這說明單旋導(dǎo)向擋板對氣泡有一定破碎作用；經(jīng)過第一層擋板后，氣泡又重新開始聚并，Db逐漸增大；同理，經(jīng)過第二層擋板時，Db再次變小。由圖還可見，大部分床高處的氣泡尺寸模擬結(jié)果與Darton 關(guān)聯(lián)式預(yù)測結(jié)果趨勢相吻合，但在流化床上界面附近，有、無兩層單旋導(dǎo)向擋板時鼓泡流化床內(nèi)的Db均出現(xiàn)明顯變小的趨勢，這是由于MS3DATA 算法基于體積計算氣泡當(dāng)量直徑，而在此區(qū)域氣泡已經(jīng)開始穿越流化床界面，從而導(dǎo)致該算法統(tǒng)計的氣泡體積偏小。圖8 表明，有擋板的鼓泡流化床內(nèi)Db小于空筒鼓泡流化床的，但減小程度相對有限，這可能是由于實驗中所使用的單旋導(dǎo)向擋板葉片數(shù)量較少、開孔率大，從而起到的氣固攪拌作用相對較弱，因此，在實際應(yīng)用中需要進一步對此單旋導(dǎo)向擋板的開孔率、葉片數(shù)量和角度等進行優(yōu)化。

圖8 氣泡平均當(dāng)量直徑沿床高分布圖Fig.8 Axial distribution of average equivalent bubble diameters

3.4 氣體停留時間分布

使用脈沖示蹤法研究有、無兩層單旋導(dǎo)向擋板時鼓泡流化床內(nèi)氣體返混現(xiàn)象和氣體停留時間分布(RTD)。計算中用空氣作為示蹤劑，在第20 s 至20.01 s 內(nèi)注入示蹤劑，令其流速與主風(fēng)流速保持一致。在此時間間隔內(nèi)，設(shè)置主風(fēng)流速為0，以保證示蹤劑的注入不會對床內(nèi)的氣固流動狀態(tài)產(chǎn)生影響。記錄每個時間步時流化床出口處示蹤劑濃度，每200 個時間步內(nèi)對示蹤劑濃度取一次平均，得到有、無兩層單旋導(dǎo)向擋板時鼓泡流化床內(nèi)的氣體停留時間分布曲線，如圖9 所示，圖中θ為對比時間(無因次)，E(θ)為停留時間分布密度函數(shù)。

圖9 氣體停留時間分布圖Fig.9 Profiles of gas residence time distribution

由圖9 可見，有、無擋板時鼓泡流化床內(nèi)的氣體停留時間分布趨勢基本相同，有單旋導(dǎo)向擋板鼓泡流化床內(nèi)的停留時間分布曲線相比于空筒鼓泡流化床略微向右偏移，這說明單旋導(dǎo)向擋板對大氣泡的破碎作用在一定程度上抑制了氣體的短路問題。進一步對停留時間分布的數(shù)學(xué)特征進行計算和分析，得到氣體在鼓泡流化床內(nèi)的平均停留時間和方差σ2列于表3，可見，有、無擋板的鼓泡流化床的返混程度均相對較小，與平推流理想反應(yīng)器接近，兩者的平均停留時間和方差基本相同，表明設(shè)置兩層單旋導(dǎo)向擋板后沒有明顯影響流化床內(nèi)的氣體返混。

表3 停留時間分布的數(shù)學(xué)特征Table 3 Mathematic characteristics of RTD

4 結(jié) 論

采用計算顆粒流體力學(xué) (CPFD) 方法耦合Igci 曳力模型對有、無兩層單旋導(dǎo)向擋板的鼓泡流化床進行模擬，模擬計算得到的軸向和徑向時均固含率分布與實驗數(shù)據(jù)較為吻合，表明此方法用來模擬單旋導(dǎo)向擋板鼓泡流化床的三維氣固流動可行。使用MS3DATA 算法統(tǒng)計和分析了該鼓泡流化床內(nèi)的氣泡平均當(dāng)量直徑沿床高的分布，結(jié)果表明，單旋導(dǎo)向擋板對大氣泡有一定的破碎作用，有擋板的鼓泡流化床內(nèi)氣泡平均尺寸小于空筒鼓泡流化床。使用空氣作為示蹤劑考察該鼓泡流化床內(nèi)的氣體返混和停留時間分布，結(jié)果表明，兩層單旋導(dǎo)向擋板的存在沒有明顯影響鼓泡流化床內(nèi)的氣體返混，有、無擋板的鼓泡流化床內(nèi)的氣體平均停留時間和方差基本相同，返混程度均與平推流理想反應(yīng)器接近。