外齒輪形顆粒催化劑上乙烯催化氧化的反應工程計算

周繼鵬，房鼎業(yè)，李　濤,2

?

外齒輪形顆粒催化劑上乙烯催化氧化的反應工程計算

周繼鵬1，房鼎業(yè)1，李濤1,2

（1華東理工大學大型工業(yè)反應器工程教育部工程研究中心，上海 200237；2華東理工大學化學工程聯(lián)合國家重點實驗室，上海 200237）

摘要:針對乙烯氧化制環(huán)氧乙烷反應體系，對外齒輪異形催化劑建立三維反應-傳質(zhì)-傳熱模型。有效擴散系數(shù)和有效熱導率均為待求解濃度場和溫度場的函數(shù)，使得偏微分方程組模型為強非線性。采用有限元算法求解，并對模型有效性進行驗證，定量研究了催化劑幾何比外表面積和內(nèi)擴散效率因子的關(guān)系。計算結(jié)果表明，幾何比外表面積為1862 m2·m-3的外齒輪催化劑內(nèi)擴散效率因子為0.1804，而幾何比外表面積為924 m2·m-3的圓柱形催化劑內(nèi)擴散效率因子為0.0993。對單個催化劑顆粒反應-傳遞現(xiàn)象的研究能定量指導催化劑設計，并為耦合反應器流體力學和催化劑反應傳遞現(xiàn)象多尺度模擬計算奠定基礎(chǔ)。

關(guān)鍵詞:乙烯氧化；異形催化劑；傳遞過程；數(shù)值模擬；內(nèi)擴散效率因子

引 言

環(huán)氧乙烷（EO）是乙烯（ET）工業(yè)中重要產(chǎn)品，在其生產(chǎn)成本中，乙烯原料占總成本的70%，故致力于開發(fā)高性能催化劑是降低乙烯單耗、提高經(jīng)濟效益最有效的手段[1-2]。乙烯氧化合成環(huán)氧乙烷過程是一個強放熱多重反應，它的選擇性由深度氧化生成CO2副反應所決定。許多研究工作所進行的乙烯氧化過程本征動力學研究都表明:生成CO2副反應導致的溫升越高，副反應的反應速率與主反應的反應速率之比越大，選擇性越低[3]。

異形催化劑在工業(yè)生產(chǎn)中，尤其是強放熱/吸熱反應體系中有廣泛的應用[4]。裝填異形多通孔催化劑的反應器具有高空隙率，因而具有更低的壓降、更好的流動性能，更少流動死區(qū)。異形催化劑幾何比表面積大，故能減少傳遞阻力，反應物能更好擴散到催化劑顆粒中，提高反應轉(zhuǎn)化率和選擇性[5-6]，從而提高催化劑內(nèi)擴散效率因子[7-9]。

目前反應器設計更多側(cè)重于反應器內(nèi)流體力學分析，而不考慮催化劑內(nèi)擴散影響[10]。而在工業(yè)應用中，催化劑內(nèi)擴散影響往往很嚴重。房鼎業(yè)等[11-12]報道了工業(yè)催化劑的效率因子較小，僅為 0.1～0.22。且對于甲烷蒸汽重整催化劑，粒內(nèi)存在較大死區(qū)[11]。求解催化劑顆粒反應-傳質(zhì)-傳熱模型常見的方法是將其簡化為一維或者二維模型求 解[11-14]，且將模型參數(shù)（有效擴散系數(shù)和有效熱導率）簡化為常數(shù)。由于模型的簡化，使得模型計算結(jié)果精度較低。本文對異形催化劑建立三維反應-傳質(zhì)-傳熱數(shù)學模型，且模型參數(shù)為待求解濃度場和溫度場函數(shù)，使得求解模型方程組呈現(xiàn)強烈的非線性。采用有限元算法（FEM），模型收斂性很好[15-16]。本文方法為耦合反應器流體力學和催化劑反應傳遞現(xiàn)象計算奠定了基礎(chǔ)，能定量指導催化劑開發(fā)和反應器多尺度模擬計算[17-18]。

1　反應-傳質(zhì)-傳熱數(shù)學模型

1.1 數(shù)學模型

Ag/Al2O3催化劑上，乙烯與氧發(fā)生如下3個不可逆反應。由原子系數(shù)矩陣法可知該體系有兩個獨立反應，本模型取反應（1）和反應（2）為獨立反應。

傳遞現(xiàn)象會影響反應進程，同時反應會產(chǎn)生溫度分布和濃度分布，影響擴散系數(shù)和熱導率，進而影響傳遞現(xiàn)象。為了準確計算模擬催化劑顆粒的反應過程，必須建立催化劑顆粒內(nèi)反應-傳質(zhì)-傳熱偏微分方程組并求解。選取ET和CO2作為關(guān)鍵組分。

反應-傳質(zhì)方程

傳熱方程

其中

式中，1r為反應1中乙烯（ET）的反應速率；r2為反應2中乙烯（ET）的反應速率； r3為反應3中環(huán)氧乙烷（EO）的反應速率。

模型中，有效擴散系數(shù)為濃度分布、溫度分布、孔徑分布的函數(shù)。有效擴散系數(shù)由分子擴散系數(shù)和Knudsen擴散系數(shù)組成[19]。

分子擴散系數(shù)

式中，i=ET、CO2，j=ET、CO2、EO、H2O、O2、N2，且ji≠；iy為催化劑內(nèi)組分的摩爾分數(shù)。

組分i，j二元擴散系數(shù)

式中，ε為催化劑孔隙率；τ為曲折因子；pt為系統(tǒng)總壓力；Mi為相對分子質(zhì)量；Vt為擴散體積[19]。

Knudsen擴散系數(shù):

式中，dpore為催化劑的平均孔徑。

催化劑顆粒的有效熱導率主要受催化劑骨架結(jié)構(gòu)和孔內(nèi)氣體熱導率影響。式（14）可用于計算催化劑顆粒的有效熱導率[20]

式中，c為與催化劑顆粒有關(guān)的參數(shù)；N為催化劑孔內(nèi)氣體組分數(shù)；λs為催化劑骨架熱導率，對Al2O3，熱導率與溫度呈線性關(guān)系[20]。

催化劑孔內(nèi)單組分氣體熱導率λfi

式中，λgi為純組分i熱導率，可關(guān)聯(lián)為溫度的函數(shù)關(guān)系式[21]。L為分子平均自由程

式中，k為 Boltzmann常數(shù)；de為分子有效直徑。

反應焓變?yōu)闇囟鹊暮瘮?shù)[20]。

反應中加入10-6級別二氯乙烷（DEC）作為抑制劑，催化劑反應動力學方程如下[22]

模型參數(shù)見文獻[22]。

1.2 物理模型及網(wǎng)格劃分

外齒輪催化劑尺寸為?6.5 mm×6.5 mm，中孔直徑1.4 mm，齒輪曲率半徑1.1 mm。

采用有限元軟件COMSOL Multiphysics。利用對稱性，對八分之一模型求解，再將計算結(jié)果映射到整個催化劑上，可以在不影響求解精度的情況下極大減少計算量。催化劑高度方向上因變量梯度在表面大，內(nèi)部小，所以網(wǎng)格采用掃掠方法，在表面網(wǎng)格密集，內(nèi)部稀疏。模型棱柱單元網(wǎng)格數(shù)為16170，求解自由度為86791。經(jīng)驗證，再細化網(wǎng)格不影響模型計算精度，網(wǎng)格獨立性很好。圖1為八分之一催化劑顆粒網(wǎng)格剖分示意圖。

圖1　八分之一催化劑網(wǎng)格示意圖Fig.1 Grid map of catalyst

1.3 邊界條件及求解方法

模擬條件為工業(yè)反應中操作條件。操作壓力為2.2 MPa，溫度為483.15 K，初始組分含量ET 28.5%、O28.5%、CO28.0%、N255.0%。

邊界條件為催化劑顆粒外表面濃度和溫度與氣相主體濃度和溫度相同。

圖1中，在催化劑外表面選擇濃度和溫度邊界條件。在對稱面上選取對稱邊界條件。采用直接迭代求解法PARDISO，收斂相對容差為0.1%。由于模型為高度非線性方程，所以先求解反應-傳質(zhì)方程式（4）、式（5）。然后將得到的解作為初值，求解反應-傳質(zhì)-傳熱方程式（4）～式（6），模型收斂性很好。

圖2　COMSOL計算值和文獻值[14]比較Fig.2 Concentration distribution of ET calculated by COMSOL comparing with literature values

2　結(jié)果與討論

2.1 模型驗證

用COMSOL軟件計算了?6.5 mm×6.5 mm，中孔直徑 1.4 mm環(huán)柱狀催化劑顆粒。圖 2是COMSOL與正交配置法計算ET濃度分布比較[12]。可見COMSOL算法與正交配置法計算結(jié)果吻合。由于文獻中正交配置法只選取 10個點對求解方程組進行離散，在對方程組離散化過程中存在較大的舍入誤差。而COMSOL可以根據(jù)求解精度對方程組進行任意階次離散，舍入誤差較文獻采用10點離散正交配置法低。故COMSOL理論上相比于文獻中采用10點離散正交配置法精度更高。

2.2 催化劑顆粒有效擴散系數(shù)與有效熱導率

由式（9）～式（16）可知，有效擴散系數(shù)和有效熱導率是待求解的濃度場和溫度場的函數(shù)。圖3(a)為ET有效擴散系數(shù)。催化劑顆粒內(nèi)ET平均有效擴散系數(shù)為1.97×10-7m2·s-1，CO2平均有效擴散系數(shù)為1.85×10-7m2·s-1，且有效擴散系數(shù)在催化劑外表面達到最大。有效擴散系數(shù)基本不變。分析可知，有效擴散系數(shù)與溫度和組分濃度有關(guān)，催化劑內(nèi)部溫度升高，會導致有效擴散系數(shù)變大，由于存在內(nèi)擴散阻力，導致有效擴散系數(shù)變小。由于溫升很小，內(nèi)擴散阻力引起的有效擴散系數(shù)變低占主要因素。所以催化劑外表面有效擴散系數(shù)最大，內(nèi)部變低。

圖3　有效擴散系數(shù)與有效熱導率Fig.3 ET effective diffusivity and effective thermal conductivity

圖3(b)為有效熱導率分布。平均有效熱導率為0.45 W·m-1·K-1。有效熱導率主要受溫度影響，內(nèi)部溫升大，有效熱導率在內(nèi)部達到最大。由于催化劑顆粒溫升不大所以有效熱導率變化不大。

2.3 ET濃度與催化劑顆粒溫度分布圖

圖4為關(guān)鍵組分ET濃度分布。由圖可知催化劑顆粒在外表面濃度降低最大，反應主要發(fā)生在外表面。對整個催化劑顆粒進行積分運算可知，ET在催化劑顆粒內(nèi)的平均濃度為 149.13 mol·m-3，CO2平均濃度為45.91 mol·m-3。

圖4　ET濃度分布Fig.4 ET concentration distribution

圖5為溫度分布。反應為放熱反應，催化劑顆粒溫度升高，最高溫度為484.58 K，溫升為1.43 K。溫升較小原因是模型求解中選取的邊界條件是溫度邊界條件。由牛頓冷卻定律，傳熱量為對流傳熱系數(shù)與溫度差的乘積，催化劑與主流體傳熱量為一有限常數(shù)，而催化劑與主流體之間溫差趨于 0（選取了溫度邊界條件），所以對流傳熱系數(shù)趨于無窮大，故相當于催化劑與主流體區(qū)域傳熱性能最好。

2.4 EO生成速率與選擇性

圖6為EO反應速率和選擇性分布。由圖6(a)可知，催化劑外表面反應速率達到最大，經(jīng)過一個很短的距離，大約0.2 mm左右，催化劑內(nèi)反應速率降低為外表面反應速率1/10，而且維持不變，可知催化劑有效反應區(qū)域為催化劑外表面。

由圖6(b)可知，催化劑選擇性在外表面達到最大。由于催化劑內(nèi)部溫度比外表面大，催化劑內(nèi)選擇性降低，提高催化劑內(nèi)擴散效率因子，降低催化劑溫升，可有效提高催化劑選擇性。

圖5　溫度分布Fig.5 Temperature distribution/K

圖6　EO反應速率和選擇性Fig.6 Reaction rate and selectivity of EO

2.5 催化劑效率因子與催化劑設計

內(nèi)擴散效率因子定義

式中，R(V)為催化劑顆粒內(nèi)反應速率，為空間坐標函數(shù)。R?為催化劑外表面反應速率。計算可知，催化劑內(nèi)擴散效率因子為0.1804，還有很大提升空間。提高催化劑效率因子，能增加催化劑利用率，從而提高反應轉(zhuǎn)化率。在生產(chǎn)規(guī)模不變條件下可以減小反應器設備尺寸，降低投資。

EO氧化反應屬于內(nèi)擴散控制，反應速率經(jīng)過外表面很短距離迅速降低90%，要提高催化劑效率因子，必須增加催化劑幾何比外表面積。圖7為催化劑選擇性和內(nèi)擴散效率因子與幾何比外表面積關(guān)系。計算了圓柱形、環(huán)柱狀、外齒輪形催化劑，得出催化劑選擇性、效率因子與幾何比外表面積的定量關(guān)系。其中，圓柱形催化劑尺寸為?6.5 mm×6.5 mm。環(huán)柱狀催化劑尺寸見2.1節(jié)。將外齒輪異形催化劑等效為等比外表面積球形催化劑顆粒時，經(jīng)計算得出兩者內(nèi)擴散效率因子相同。

圖7　催化劑選擇性和內(nèi)擴散效率因子與幾何比外表面積關(guān)系Fig.7 Selectivity and internal effectiveness factor of catalyst with geometrical specific surface area

3 結(jié) 論

（1）乙烯環(huán)氧化反應屬于內(nèi)擴散控制，且擴散阻力區(qū)在距離催化劑外表面約0.2 mm左右。在擴散阻力區(qū)反應速率下降90%，越過擴散阻力區(qū)，催化劑內(nèi)部反應速率基本維持不變。

（2）幾何比外表面積為1862 m2·m-3的外齒輪催化劑內(nèi)擴散效率因子為 0.1804，EO選擇性為86.00%，而幾何比外表面積為924 m2·m-3的圓柱形催化劑內(nèi)擴散效率因子為 0.0993，EO選擇性為84.12%。

（3）催化劑幾何比外表面積能作為評價催化劑內(nèi)擴散影響的重要指標。且催化劑內(nèi)擴散效率因子幾乎隨催化劑幾何比外表面積呈線性遞增關(guān)系。

符 號 說 明

ci——組分i濃度，i=ET,CO2，mol·m-3

ci,b——組分i在主流體濃度，mol·m-3

ci,??——組分i在催化劑表面濃度，mol·m-3

Deff(CO2) ——CO2有效擴散系數(shù)，m2·s-1

Deff(ET) ——ET有效擴散系數(shù)，m2·s-1

Deff,k,i——Knudsen擴散系數(shù)，m2·s-1

Deff,m,i——分子擴散系數(shù)，m2·s-1

dpore——催化劑平均孔徑，nm

ΔHi——反應焓變，i=1,2,3，J·mol-1

Ki——吸附常數(shù)

k ——Boltzmann常數(shù)，k =1.38×10-23J·K-1

ki——反應速率常數(shù)

Mi——相對分子質(zhì)量

pt——總壓，MPa

ri——反應i的反應速率，i=1,2,3，kmol·kg-1·h-1

r(CO2) ——CO2總反應速率，mol·kg-1·s-1

r(ET) ——ET總反應速率，mol·kg-1·s-1

T ——催化劑溫度，K

Tb——主流體溫度，K

T??——催化劑外表面溫度，K

Vi——擴散體積

yi——組分i的摩爾分數(shù)，i=ET, CO2, EO, H2O, O2, N2

ε ——催化劑孔隙率

η ——催化劑內(nèi)擴散效率

因子

λeff——有效熱導率，W·m-1·

K-1

λfi——催化劑內(nèi)單組分氣體熱導率，W·m-1·K-1

λgi——純組分氣體熱導率，W·m-1·K-1

ρcat——催化劑密度，kg·m-3

τ ——催化劑曲折因子

References

[1] 李濤, 樊蓉蓉, 朱炳辰. 薄壁多通孔環(huán)氧乙烷合成銀催化劑工程研究-外冷管式反應器的數(shù)學模擬及比較 [J]. 石油化工, 2008, 34: 817-818.

LI T, FAN R R, ZHU B C. Engineering design of thin-walled silver catalyst of ethylene to ethylene oxide: mathematical modeling and compare of external cooling tube reactor [J]. Petrochemical Technology, 2008, 34: 817-818.

[2] SALMI T, CARUCCI J H, ROCHE M, et al. Microreactors as tools in kinetic investigations: ethylene oxide formation on silver catalyst [J]. Chem. Eng. Sci., 2013, 87: 306-314.

[3] 李紅權(quán). 乙烯氧化制環(huán)氧乙烷反應的熱力學分析及工藝改進方向[J]. 工業(yè)催化, 1999, 7 (1): 33-37.

LI H Q. Thermodynamics for synthesis of ethylene oxide through ethylene oxidation and the trends in improvement on its process condition [J]. Industrial Catalysis, 1999, 7 (1): 33-37.

[4] DEROUANE E G, PARMON V, LEMOS F. Sustainable Strategies for the Upgrading of Natural Gas in Sustainable Strategies for the Upgrading of Natural Gas: Fundamentals, Challenges and Opportunities [M]. Netherlands: Springer, 2005: 413-420.

[5] SIE S T, KRISHNA R. Process development and scale up (Ⅱ): Catalyst design strategy [J]. Rev. Chem. Eng., 1998, 14 (3): 159-202.

[6] DIXON A G, NIJEMEISLAND M, STITT E H. Packed tubular reactor modeling and catalyst design using computational fluid dynamics [J]. Advances in Chemical Engineering, 2006, 31: 307-389.

[7] BRUNO S P, BARRETO G F, GONZALEZ M G. Effect of the geometric characteristics of commercial catalysts for steam reforming [J]. Chem. Eng. J., 1988, 39 (3): 147-156.

[8] BAEK S M, KANG J H, LEE K J, et al. A numerical study of the effectiveness factors of nickel catalyst pellets used in steam methane reforming for residential fuel cell applications [J]. Int. J. Hydrogen Energy, 2014, 39 (17): 9180-9192.

[9] ALBERTON A L, SCHWAAB M, FONTES C E, et al. Hybrid modeling of methane reformers (Ⅰ): A metamodel for the effectiveness factor of a catalyst pellet with complex geometry [J]. Ind. Eng. Chem. Res., 2009, 48 (21): 9369-9375.

[10] DIXON A G. CFD study of effect of inclination angle on transport and reaction in hollow cylinder catalysts [J]. Chem. Eng. Res. Des., 2014, 92 (7): 1279-1295.

[11] 房鼎業(yè), 孫啟文, 朱炳辰. 非等溫甲烷蒸汽轉(zhuǎn)化反應催化劑的效率因子 [J]. 高?；瘜W工程學報, 1991, 3 (5): 225-231.

FANG D Y, SUN Q W, ZHU B C. The effectiveness factors for non-isothermal methane steam-reforming catalyst [J]. Journal of Chemical Engineering of Chinese Universities, 1991, 3 (5): 225-231.

[12] 高崇, 潘銀珍, 朱炳辰. 環(huán)柱狀催化劑內(nèi)強放熱復合反應-傳質(zhì)-傳熱耦合過程研究（Ⅲ）: 數(shù)學模型的求解及實驗驗證 [J]. 化工學報, 1998, 49 (5): 617-623.

GAO C, PAN Y Z, ZHU B C. Study on coupling process of multiple reaction-mass transfer-heat transfer in hollow cylindrical catalysts with strong heat effect (Ⅲ): Solution and experimental verification of mathematical model [J]. Journal of Chemical Industry and Engineering (China), 1998, 49 (5): 617-623.

[13] LI T, XU M S, ZHU B C, et al. Reaction-diffusion model for irregularly shaped ammonia synthesis catalyst and its verification under high pressure [J]. Ind. Eng. Chem. Res., 2009, 48 (19): 8926-8933.

[14] PAN T S, ZHU B C. Study on diffusion-reaction process inside a cylindrical catalyst pellet [J]. Chem. Eng. Sci., 1998, 53 (5): 933-946. [15] HAYES R E, KOLACZKOWSKI S T, THOMAS W J. Finite element model for a catalytic monolith reactor [J]. Computers & Chemical Engineering, 1992, 16 (7): 645-657.

[16] ROY S, HEIBEL A K, LIU W, et al. Design of monolithic catalysts for multiphase reactions [J]. Chem. Eng. Sci., 2004, 59 (5): 957-966.

[17] KO?í P, NOVáK V, ?TěPáNEK F, et al. Multi-scale modelling of reaction and transport in porous catalysts [J]. Chem. Eng. Sci., 2010, 65 (1): 412-419.

[18] KULKARNI K, MOON J, ZHANG L B, et al. Multi-scale modeling and solution multiplicity in catalytic pellet reactors [J]. Ind. Eng. Chem. Res., 2008, 47 (22): 8572-8581.

[19] FOGLER H S. Elements of Chemical Reaction Engineering [M]. 4th ed. Pearson Education, 2006: 813-827.

[20] 高崇, 朱英, 李樹森, 等. 環(huán)柱狀催化劑內(nèi)強放熱復合反應-傳質(zhì)-傳熱耦合過程研究（Ⅰ）: 催化劑曲折因子有效熱導率的測定 [J].化工學報, 1998, 49 (5): 601-609.

GAO C, ZHU Y, LI S S, et al. Study on coupling process of multiple reaction-mass transfer-heat transfer in hollow cylindrical catalysts with strong heat effect (Ⅰ): Experimental measurement of tortuosity factors and effective thermal conductivity of catalyst pellets [J].Journal of Chemical Industry and Engineering (China), 1998, 49 (5): 601-609.

[21] 時鈞, 汪家鼎, 余國琮, 等. 化學工程手冊 [M]. 2版. 北京: 化學工業(yè)出版社, 1996: 1-141-1-145.

SHI J, WANG J D, YU K T, et al. Chemical Engineers Handbook [M]. 2nd ed. Beijing: Chemical Industry Press, 1996: 1-141-1-145.

[22] 高崇, 潘銀珍, 朱炳辰. 環(huán)柱狀催化劑內(nèi)強放熱復合反應-傳質(zhì)-傳熱耦合過程研究（Ⅱ）: 本征反應動力學及反應-傳質(zhì)-傳熱耦合過程數(shù)學模型 [J]. 化工學報, 1998, 49 (5): 610-616.

GAO C, PAN Y Z, ZHU B C. Study on coupling process of multiple reaction-mass transfer-heat transfer in hollow cylindrical catalysts with strong heat effect (Ⅱ): Intrinsic kinetics and mathematical model of coupling process of reaction-mass transfer-heat transfer [J]. Journal of Chemical Industry and Engineering (China), 1998, 49 (5): 610-616.

2016-01-19收到初稿，2016-03-08收到修改稿。

聯(lián)系人:李濤。第一作者:周繼鵬（1990—），男，碩士研究生。

Received date: 2016-01-19.

中圖分類號:TQ 028.8

文獻標志碼:A

文章編號:0438—1157（2016）07—2808—07

DOI:10.11949/j.issn.0438-1157.20160081

Corresponding author:Prof. LI Tao, tli@ecust.edu.cn

Reaction engineering calculation for ethylene catalytic oxidation over gear-shaped catalysts

Abstract:For the reaction of ethylene catalytic oxidation to ethylene oxide over a gear-shaped catalyst particle, a three-dimensional reaction-mass transfer-heat transfer model is developed to calculate effectiveness factor of the catalyst. Effective diffusivities and effective thermal conductivity are the function of temperature and concentration distribution in catalyst particle. The finite element method is employed to solve the set of highly nonlinear partial differential equations, and the values obtained agree well with those given by literatures. Geometric external surface area is a key parameter for the effectiveness factor. The effectiveness factor is 0.1804 for gear-shaped catalyst with geometric external surface area of 1862 m2·m-3, and 0.0993 for cylindrical catalyst with geometric external surface area of 924 m2·m-3. The quantitative calculation for reaction-transport phenomena in a single catalyst particle could be a guide for catalyst design and a basis for multi-scale simulation of coupling hydromechanics and catalytic reaction-transport phenomena.

Key words:ethylene oxidation; irregular shape catalyst; transport processes; numerical simulation; internal effectiveness factor