沈文豪, 張亞新, 趙 玲

(1. 新疆大學(xué) 化學(xué)化工學(xué)院, 新疆 烏魯木齊 830046;2. 華東理工大學(xué) 化學(xué)工程聯(lián)合國家重點實驗室, 上海 200237)

1 前 言

我國煤炭資源豐富，國家多項政策指明了煤炭行業(yè)供給側(cè)改革的發(fā)展方向是清潔高效利用[1]。大規(guī)模煤制天然氣項目已在新疆、內(nèi)蒙古等富煤地區(qū)投產(chǎn)，該工藝中甲烷化反應(yīng)段對甲烷產(chǎn)率影響最大。副反應(yīng)中甲烷吸熱會產(chǎn)生絲狀炭附著在催化劑內(nèi)孔道表面，長時間的積累會導(dǎo)致催化劑內(nèi)孔道嚴重堵塞[2-3]。張沛等[4]采用掃描電子顯微鏡(SEM)等表征方法，對商用選擇性催化還原(SCR)脫硝催化劑的失活和再生進行研究，發(fā)現(xiàn)砷和硫酸鹽嚴重堵塞了催化劑微孔道，酸洗再生能有效脫除砷和硫酸鹽，提高SCR 脫硝催化劑活性。KOO 等[5]采用不同Mg/Al 混合比控制的載體，提高了甲烷重整制合成氣催化劑的抗結(jié)焦能力。DEMICHELI 等[6-7]在實驗中測定了甲烷在Ni/Al2O3催化劑上形成積炭的反應(yīng)動力學(xué)方程，使數(shù)值模擬方法研究甲烷化反應(yīng)積炭過程成為可能。

采用 CFD 方法對固定床反應(yīng)器內(nèi)“三傳一反”(動量傳遞、質(zhì)量傳遞、熱量傳遞、化學(xué)反應(yīng)過程[8])規(guī)律進行研究能夠節(jié)省大量的實驗和人工成本，并獲得反應(yīng)器內(nèi)部更詳細的特征參數(shù)[9-10]。CHEIN 等[11]采用數(shù)值模擬方法研究了固定床管式反應(yīng)器中甲烷化反應(yīng)的影響因素，發(fā)現(xiàn)CO 甲烷化反應(yīng)對溫度最敏感，提高反應(yīng)壓力和降低空速能提高甲烷產(chǎn)率。SON 等[12]采用數(shù)值模擬方法研究了操作條件對積炭效應(yīng)的影響，但模擬也僅停留在反應(yīng)器宏觀尺度和催化劑-流體域介觀尺度。對固定床反應(yīng)器的模擬研究通常將催化劑床層簡化為多孔介質(zhì)模型或者規(guī)則排列的小球[13-14]，耦合催化劑內(nèi)和床層流體域的多尺度研究也幾乎只在國外文獻中報道，DIXON 等[15-17]在催化劑內(nèi)部反應(yīng)過程的研究中做了很大貢獻，但大都針對于蒸汽重整等反應(yīng)。掌握降低積炭效應(yīng)的技術(shù)，必需建立更趨于真實的隨機堆積固定床床層，并將介觀流體域和微觀催化劑域的多物理場耦合，研究多尺度下甲烷化反應(yīng)催化劑積炭形成過程。

本文首先使用DEM 軟件生成顆粒(催化劑)直徑dp=0.012 m、床層與顆粒直徑比N=4、高徑比H/D=5的隨機堆積球形顆粒床層，并提取顆粒位置信息，通過編程在CFD 軟件中建立床層幾何模型，將計算域分為床層流體域和催化劑顆粒域，建立“三傳一反”多物理場耦合數(shù)學(xué)模型，通過數(shù)值模擬分析了微觀顆粒域的積炭形成規(guī)律。最后討論入口溫度400、450、500 K，催化劑溫度770、810、850 K，入口流速0.005、0.010、0.015 m·s-1和入口氫氣占比0、20%、40%分別對甲烷化反應(yīng)催化劑積炭過程的影響。

2 模型的建立與驗證

2.1 幾何模型

在隨機堆積顆粒的固定床反應(yīng)器中，較大管徑比將導(dǎo)致床層顆粒數(shù)倍增，給數(shù)值求解帶來巨大的計算量，而小管徑比固定床反應(yīng)器在強放熱、吸熱反應(yīng)和實驗室研究階段應(yīng)用也比較廣泛，兩者徑向流動均勻性存在差異，但最終整體特征規(guī)律仍具有一定相似性。

本文建立了小管徑比固定床反應(yīng)器床層，結(jié)構(gòu)如圖1 所示，甲烷從上端進入，下端流出，床層總高為12 dp，直徑為4 dp，將床層分為入口、催化劑堆積和出口3 個區(qū)域，入口區(qū)域高2 dp，催化劑堆積區(qū)域高5 dp，出口區(qū)域高5 dp。使用DEM 軟件在催化劑堆積區(qū)域生成直徑為dp=12 mm 的球型顆粒，并提取顆粒位置信息，建立幾何模型如圖2 所示，DEM 參數(shù)設(shè)置如表1。

圖1 床層結(jié)構(gòu)示意圖Fig.1 Schematic diagram of the bed structure

圖2 床層堆積模型Fig.2 Schematic diagram of the bed packing model

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

2.2.1 控制方程

甲烷化反應(yīng)在顆粒內(nèi)放出的熱量成為積炭反應(yīng)熱源，顆粒內(nèi)由無數(shù)微孔道組成，現(xiàn)有建模手段和計算能力很難實現(xiàn)顆粒內(nèi)成千上萬微孔道的建模，積炭又覆蓋在內(nèi)孔道表面，改變了催化劑的孔隙率。因此，做出以下合理假設(shè)：

(1) 流體為理想氣體；

(2) 流體為不可壓縮、穩(wěn)態(tài)的層流；

(3) 催化劑顆粒為多孔介質(zhì)模型；

(4) 反應(yīng)所需要的熱量來自顆粒；

(5) 積炭對催化劑活性的影響非常小可忽略，主要影響孔隙率的變化。

基于以上假設(shè)，建立多物理場耦合數(shù)學(xué)模型，控制方程如下：

(1) 動量守恒方程：

表1 DEM 參數(shù)設(shè)置Table 1 Parameters for DEM

積炭會降低催化劑孔隙率，進而影響滲透率，初始滲透率采用多孔介質(zhì)Kozeny-Carman 滲透率模型，描述催化劑孔隙率和滲透率變化的表達式[18]如下：

式中：MC為積炭的摩爾質(zhì)量，0.012 kg·mol-1；ρdepos為積炭的密度，2 000 kg·m-3。

2.2.3 邊界條件與計算方法

入口初始甲烷摩爾分數(shù)為1，溫度為500 K，平均流速為0.015 m·s-1，出口壓力為101 325 Pa，將反應(yīng)源添加在催化劑顆粒域，壁面設(shè)置為無滑移和熱絕緣邊界條件，催化劑區(qū)域設(shè)置為熱源，溫度為850 K。催化劑初始孔隙率采用Dixon 實驗測得的Ni/Al2O3球型催化劑孔隙率0.44[19]。

考慮到多物理場耦合中非線性方程組求解的難度，首先將物理場解耦，單獨計算穩(wěn)態(tài)流場，最后將穩(wěn)態(tài)流場與其他物理場耦合計算，采用直接迭代PARDISO 算法求解，瞬態(tài)求解8 000 s，步長為60 s，使用默認松弛因子值。

2.3 驗證

2.3.1 幾何模型準確性驗證

早期學(xué)者提出球型固定床床層整體空隙率[20]與徑向空隙率分布[21]的經(jīng)驗公式，并與實驗數(shù)據(jù)吻合較好，得到學(xué)者們的廣泛認可。整體空隙率經(jīng)驗公式如下：

經(jīng)計算，本文中DEM 堆積的床層整體空隙率為0.47，擬合公式計算空隙率為0.44，相差3%，誤差主要來源于堆積的隨機性以及計算空隙率時床層高度的取值。

徑向空隙率經(jīng)驗公式如下：

圖3 徑向空隙率截取示意圖Fig.3 Schematic diagram of the radial void fraction

圖4 徑向空隙率與文獻對比Fig.4 Comparison of the radial void fraction results with literature

2.3.2 網(wǎng)格無關(guān)性驗證

由于反應(yīng)主要發(fā)生在催化劑堆積區(qū)，為減少計算量將進出口區(qū)域劃分為六面體網(wǎng)格，催化劑堆積區(qū)均劃分為四面體網(wǎng)格，筒體和顆粒邊界劃分3 層邊界層。反應(yīng)、擴散等行為在顆粒域更為復(fù)雜，對其網(wǎng)格進行加密處理，由于在顆粒中心施加溫度邊界，為更精確計算邊界處參數(shù)，將顆粒中心網(wǎng)格再次細化，網(wǎng)格劃分如圖5 所示。直接劃分網(wǎng)格會導(dǎo)致壁面接觸點附近網(wǎng)格質(zhì)量急劇降低，將所有顆粒直徑減小1%，床層空隙率只下降1.2%，對結(jié)果影響較小。采用同樣的網(wǎng)格劃分方法劃分1 040 000、1 570 000 和2 060 000網(wǎng)格，比較3 種網(wǎng)格床層壓降變化，如圖6 所示，結(jié)果表明，繼續(xù)細化網(wǎng)格對結(jié)果并無較大影響，為降低計算成本，將網(wǎng)格數(shù)量確定為1 570 000。

圖5 網(wǎng)格劃分Fig.5 Mesh generation

圖6 網(wǎng)格數(shù)量對壓力變化的影響Fig.6 Effects of grid number on pressure variation

3 結(jié)果與討論

3.1 模擬有效性驗證

KUVSHINOV 等[23]在填裝Ni/Al2O3催化劑，外部持續(xù)加熱維持等溫狀態(tài)的管式反應(yīng)器中進行甲烷積碳反應(yīng)實驗，測得H2占比為0、0.15、0.30、0.40 時，溫度從760 K 增加到860 K 時，反應(yīng)速率的變化規(guī)律。本文采用CFD 方法同樣進行等溫反應(yīng)的模擬，將模擬值與文獻值對比，如圖7 所示，在氫氣質(zhì)量分數(shù)為0 時，模擬值與試驗值有很小的誤差，其余參數(shù)下二者結(jié)果基本吻合，證明模擬方法準確可靠。

圖7 模擬值與實驗值對比Fig.7 Comparison of simulation with experimental values

圖8 多物理場分布Fig.8 Multiple physical field distribution

3.2 多物理場分布

將流體域與顆粒域多物理場耦合計算，得到如圖8 所示的多物理場分布示意圖。從速度流線圖可以看出，流體流動在出入口區(qū)很穩(wěn)定，在床層堆積區(qū)會沿顆粒表面無規(guī)則流動，局部速度會增大到入口流速的10 倍，顆粒尾部還會出現(xiàn)尾流和回流現(xiàn)象，流體流經(jīng)顆粒多孔介質(zhì)區(qū)域，在剛接觸顆粒、離開顆粒和顆粒內(nèi)的速度幾乎為0。反應(yīng)器壁面(右)是壓力分布，壓降發(fā)生在床層堆積區(qū)，進出口壓降為0.104 Pa。反應(yīng)器壁面(左)是溫度分布，入口甲烷溫度為500 K，隨著距入口距離的增加，床層區(qū)溫度逐漸上升，第1 層顆粒溫度從500 K 迅速上升至850 K，到第2 層顆粒時溫度基本達到850 K。顆粒表面是反應(yīng)剛開始時甲烷濃度分布，甲烷在顆粒上半部分濃度高，下半部分濃度低，在入口區(qū)域濃度為14.4 mol·m-3，底部最小濃度能達到 14 mol·m-3。

3.3 場分布隨時間的變化

3.3.1 熱態(tài)場分布

分別提取了不同時刻床層中心截面溫度、甲烷摩爾分數(shù)、反應(yīng)速率和顆粒中心截面孔隙率的分布。圖9 為不同時刻床層中心截面溫度分布，當(dāng)t=0 s 時，入口溫度為500 K，顆粒中心溫度為850 K，其他區(qū)域溫度均為293.15 K。當(dāng)t=60 s 時，各區(qū)域溫度基本達到穩(wěn)定，床層堆積區(qū)和流體區(qū)溫度迅速上升至850 K，但由于入口區(qū)不斷有500 K 的甲烷流入，使入口區(qū)和床層堆積區(qū)產(chǎn)生500～850 K 的溫度梯度。

圖9 溫度隨時間的變化Fig.9 Temperature profiles as a function of time

圖10 甲烷摩爾分數(shù)隨時間的變化Fig.10 Profiles of methane mole fraction as a function of time

圖10 為甲烷摩爾分數(shù)隨時間的變化，當(dāng)t=60 s時，顆粒內(nèi)少部分甲烷已經(jīng)開始反應(yīng)。當(dāng)t=2 000 s時，由于甲烷向催化劑內(nèi)層擴散，擴散阻力增大，部分甲烷反應(yīng)生成了碳和氫氣，甲烷摩爾分數(shù)開始減小，最小為0.93；氣流方向至上而下，單個顆粒上半部分甲烷供給更快，導(dǎo)致甲烷在顆粒上半部分多、下半部分少。當(dāng)t=4 000 s 時，顆粒內(nèi)層甲烷摩爾分數(shù)最小為0.90，此時外層孔道已經(jīng)嚴重堵塞，顆粒內(nèi)部傳質(zhì)受外部氣流影響減小，導(dǎo)致甲烷在顆粒內(nèi)近似呈圓環(huán)狀減小。當(dāng)t=8 000 s 時，顆粒幾乎完全堵塞，甲烷在顆粒內(nèi)呈圓環(huán)狀減小的趨勢更明顯，不同的是顆粒內(nèi)層甲烷摩爾分數(shù)最小已下降至0.85。

圖11 為反應(yīng)速率隨時間的變化，當(dāng)t=60 s 時，最小反應(yīng)速率為 0.006 1 mol·m-3·s-1；當(dāng) t=2 000 s時，顆粒內(nèi)甲烷被消耗減少，最小反應(yīng)速率減小至0.005 3 mol·m-3·s-1；當(dāng) t=4 000 、8 000 s 時，最小反應(yīng)速率分別為 0.004 8、0.004 0 mol·m-3·s-1，分布規(guī)律同甲烷摩爾分數(shù)分布規(guī)律相同，反應(yīng)速率在顆粒內(nèi)也近乎呈圓環(huán)狀減小。

圖11 反應(yīng)速率隨時間的變化Fig.11 Variation of reaction rate as a function of time

3.3.2 孔隙率分布

如圖12 所示，反應(yīng)從60 s 進行到8 000 s，顆粒內(nèi)平均孔隙率值從0.44 下降到0.04左右。隨反應(yīng)的進行，顆粒外層比內(nèi)層積炭更嚴重，當(dāng)t=4 000 s 時，顆粒內(nèi)、外層孔隙率值最大相差0.05，對比圖12 和11，發(fā)現(xiàn)反應(yīng)速率越快的區(qū)域，孔隙率越小。提取了相同時間間隔(1 000 s)床層中心軸線上顆粒內(nèi)的孔隙率分布，如圖 13 所示，表明反應(yīng)越充分，顆粒內(nèi)孔隙率分布越不均勻，呈現(xiàn)出外層低、內(nèi)層高的規(guī)律。提取前 8 000 s每1 000 s 顆粒內(nèi)平均孔隙率下降值，如圖14 所示，表明孔隙率下降程度隨時間的增加在不斷減緩，4 000 s 前，催化劑未完全堵死，反應(yīng)劇烈，平均每1 000 s 顆粒內(nèi)平均孔隙率值下降0.08；4 000 s 之后，顆粒內(nèi)平均孔隙率值下降到0.1 以下，積炭反應(yīng)減慢，平均每1 000 s 顆粒內(nèi)平均孔隙率值下降0.025。

圖12 中心截面孔隙率隨時間的變化Fig.12 Variation of central section porosity as a function of time

圖13 中心截線孔隙率隨時間的變化Fig.13 Variation of central line porosity as a function of time

圖14 相同時間間隔孔隙率下降趨勢Fig.14 Decrease of porosity under same time intervals

3.4 穩(wěn)態(tài)時場分布規(guī)律

根據(jù)3.3 節(jié)描述，發(fā)現(xiàn)時間進行到8 000 s 時反應(yīng)程度已經(jīng)很小，為增加結(jié)論的準確性，將計算時間增加到12 000 s，提取穩(wěn)態(tài)時甲烷、氫氣摩爾分數(shù)和顆粒內(nèi)孔隙率分布云圖。如圖15 所示，當(dāng)t=12 000 s 時，流體區(qū)甲烷摩爾分數(shù)高達 99%，氫氣含量很少，從顆粒外層到內(nèi)層，甲烷含量越來越低，氫氣含量越來越高。如圖16 所示，當(dāng)t=12 000 s 時顆粒內(nèi)平均孔隙率值為0.015，此時顆粒內(nèi)、外層孔隙率值相差0.02，顆粒外層積炭仍然更嚴重。這種現(xiàn)象是顆粒內(nèi)擴散和反應(yīng)同時造成的，顆粒外層反應(yīng)速率大，生成的氫氣和碳也越多，同時大量甲烷也會稀釋氫氣，使顆粒外層的氫氣減少；顆粒內(nèi)層甲烷內(nèi)擴散的同時也不斷地反應(yīng)產(chǎn)生氫氣，最后集中在顆粒中心，內(nèi)層反應(yīng)程度反而越來越小，積炭程度也越來越小。

圖15 甲烷和氫氣摩爾分數(shù)對比（t = 12 000 s）Fig.15 Comparison of mole fraction of methane and hydrogen (t = 12 000 s)

圖16 孔隙率分布(t = 12 000 s)Fig.16 Porosity distribution (t = 12 000 s)

3.5 不同操作參數(shù)對積炭的影響

實際生產(chǎn)中，最容易控制的是操作參數(shù)，且通過改變操作參數(shù)來提高生產(chǎn)效率經(jīng)濟成本最低。由圖16 可知，每一個顆粒內(nèi)孔隙率分布都表現(xiàn)出對稱性，表明顆粒中心截線上孔隙率的變化能夠反映顆粒內(nèi)平均孔隙率的變化規(guī)律。因此，提取當(dāng)t=12 000 s，穩(wěn)態(tài)時，床層中心軸線上顆粒內(nèi)孔隙率的變化，研究了不同入口氫氣占比、溫度、流速和催化劑溫度對積炭效應(yīng)的影響，如圖17～20 所示。

如圖17 所示，入口組分無氫氣時，顆粒內(nèi)平均孔隙率值為0.02，積碳非常嚴重，加入氫氣，由于增加了生成物使促使反應(yīng)逆向進行，顆粒內(nèi)平均孔隙率值明顯提高；當(dāng)氫氣占比 40%時，顆粒內(nèi)平均孔隙率值增大到0.29，同時，加入氫氣會使顆粒內(nèi)、外層甲烷濃度梯度減小，積炭分布越均勻。如圖18 所示，顆粒內(nèi)溫度為770 K 時，顆粒內(nèi)平均孔隙率值為0.23；溫度從770 K 上升到810 K，顆粒內(nèi)平均孔隙率下降約0.15；溫度從810 K 升高到850 K，顆粒內(nèi)平均孔隙率值下降約0.06，表明催化劑升溫會加速積炭的形成，繼續(xù)升溫對積炭的影響越來越小。如圖19 所示，入口速度從0.005 m·s-1增加到0.010、0.015 m·s-1對顆粒內(nèi)平均孔隙率幾乎無影響，但顆粒中心處孔隙率值會下降0.001 左右，這是由于增大流速會增大反應(yīng)器軸向壓差，在更大壓力驅(qū)動下相同時間甲烷內(nèi)擴散越多，生成的積炭也越多。如圖20 所示，入口溫度的升高只導(dǎo)致第一層顆粒中心的孔隙率減小，因為入口溫度越高，進入床層的氣體能夠更快達到反應(yīng)區(qū)溫度，更快生成積炭。這與溫度場分布規(guī)律一致，由圖10 可知，第1 層顆粒溫度變化梯度較大，第2 層溫度基本達到穩(wěn)定。

圖17 入口氫氣含量對孔隙率的影響Fig.17 Effects of inlet hydrogen contents on porosity

圖18 催化劑溫度對孔隙率的影響Fig.18 Effects of catalyst temperature on porosity

圖19 入口速度對孔隙率的影響Fig.19 Effects of inlet velocities on porosity

圖20 入口溫度對孔隙率的影響Fig.20 Effects of inlet temperature on porosity

圖21 更直觀的反應(yīng)了4 種操作參數(shù)對顆粒內(nèi)平均孔隙率的影響，入口氫氣從0 增加到40%，顆粒內(nèi)平均孔隙率值從0.02 上升至0.29，增加了0.27。催化劑溫度從770 K 增加到850 K，顆粒內(nèi)平均孔隙率值從0.24變?yōu)?.02，減小了0.22。當(dāng)入口流速從0.005 m·s-1增加到0.015 m·s-1時，顆粒中心孔隙率減小 0.01，平均孔隙率幾乎無變化。當(dāng)入口溫度從400 K 增加到500 K，第1 層顆粒中心孔隙率減小0.005，其他顆粒中心、平均孔隙率均無變化。

圖21 孔隙率的影響對比(t = 12 000 s)Fig.21 Comparison of effects of different factors on porosity (t = 12 000 s)

4 結(jié) 論

建立隨機堆積球形催化劑固定床床層幾何模型和甲烷化反應(yīng)積炭過程的多物理場有限元模型，模擬了床層介尺度流體域和微尺度顆粒域的“三傳一反”現(xiàn)象，解析了催化劑內(nèi)部積炭形成過程，并討論不同操作參數(shù)對積炭效應(yīng)的影響，得到以下結(jié)論：

(1) 以模擬的方法揭示了積炭是由催化劑外層到內(nèi)層逐漸形成的，積碳主要發(fā)生在反應(yīng)初期，最先產(chǎn)生在催化劑外表面，積炭發(fā)生在內(nèi)孔道表面相對滯后，導(dǎo)致顆粒內(nèi)孔隙率從外到內(nèi)呈圓環(huán)狀減?。?/p>

① 當(dāng)t=4 000 s 時，顆粒內(nèi)、外層孔隙率值最大相差0.05；

② 當(dāng)t=0 s 時，催化劑孔隙率值為0.44，反應(yīng)進行4 000 s 時減小至0.1，8 000 s 減小至0.025。

(2) 研究了不同操作參數(shù)對積炭效應(yīng)的影響，發(fā)現(xiàn)增加入口氫氣含量、降低催化劑溫度能明顯降低積炭效應(yīng)，降低入口流速和入口溫度能輕微減緩催化劑中心積炭程度：

① 當(dāng)入口氫氣從0 增加到40%時，顆粒內(nèi)平均孔隙率值增大0.27；

② 當(dāng)催化劑溫度從770 K 增加到850 K 時，顆粒內(nèi)平均孔隙率值減小0.22；

③ 當(dāng)入口流速從0.005 m·s-1增加到0.015 m·s-1時，每層催化劑中心孔隙率值減小0.01；

④ 當(dāng)入口溫度從400 K 增加到500 K 時，第1 層催化劑中心孔隙率值減小0.005。

符號說明：

ci— 濃度，mol·m3

Cp— 恒壓比熱容，J·kg-1·K-1

dp— 顆粒直徑，m

D — 床層直徑，m

Dij— 擴散系數(shù)，m2·s-1

kpor— 無量綱系數(shù)

k — 速率常數(shù)，mol·m-3·s-1·Pa-1

Kp— 平衡常數(shù)，Pa

KH— 吸附平衡常數(shù)，Pa-0.5

M — 摩爾質(zhì)量，kg·mol-1

N — 床層直徑與顆粒直徑比

p — 壓力，Pa

Q — 熱源，W·m3

r — 反應(yīng)速率，mol·m-3·s-1

R — 氣體常數(shù)，J·mol-1·K-1

T — 溫度，K

u — 速度，m·s-1

w — 質(zhì)量分數(shù)

xi— 摩爾分數(shù)

ε — 催化劑孔隙率

ρ — 混合物密度，kg·m-3

κ — 滲透率，m2

λe— 有效導(dǎo)熱率，W·m-1·K-1

ε0— 初始孔隙率

κ0— 初始滲透率，m-2

下標

ave — 平均

cat — 催化劑

depos — 積炭

i,j — 各物質(zhì)名稱

mix — 混合物