基于雙尺度耦合模擬的甲烷化催化劑多級(jí)孔結(jié)構(gòu)探究

葉 銀，張亞新

新疆大學(xué)化工學(xué)院，新疆 烏魯木齊 830046

我國(guó)是天然氣消費(fèi)大國(guó)，天然氣消耗量逐年增加，對(duì)進(jìn)口天然氣的依賴(lài)度偏高。二氧化碳加氫制甲烷工藝既能實(shí)現(xiàn)溫室氣體高效綠色資源化利用，增加天然氣的供應(yīng)，又能為棄風(fēng)、棄光等不穩(wěn)定新能源電力的轉(zhuǎn)化存儲(chǔ)提供一種有效方式[1]。所以從環(huán)保、能源戰(zhàn)略安全以及能源產(chǎn)業(yè)結(jié)構(gòu)優(yōu)化的角度出發(fā)，CO2甲烷化技術(shù)具有極大的研究?jī)r(jià)值。

甲烷化過(guò)程是合成天然氣工藝中的重要一環(huán)，具有反應(yīng)迅速、熱效應(yīng)強(qiáng)的特點(diǎn)[2]，要讓反應(yīng)能夠順利進(jìn)行，并且兼顧反應(yīng)物的轉(zhuǎn)化率和甲烷的產(chǎn)率，就必須考慮催化劑的作用。這些年來(lái)，從催化劑的制備和成型角度出發(fā)，對(duì)活性中心、助劑以及異形結(jié)構(gòu)等方面已進(jìn)行了大量的研究[3-5]，但對(duì)顆粒中的擴(kuò)散限制關(guān)注較少。程永剛等[6]的研究表明，在反應(yīng)過(guò)程中，甲烷化催化劑的主要反應(yīng)區(qū)域?yàn)轭w粒的外層，催化劑內(nèi)部空間未能得到充分利用，通過(guò)一定的手段降低顆粒的內(nèi)擴(kuò)散影響將對(duì)反應(yīng)過(guò)程產(chǎn)生積極影響。早期的研究[7-8]表明，催化劑顆粒最佳孔結(jié)構(gòu)的設(shè)計(jì)非常重要。隨著實(shí)驗(yàn)技術(shù)的發(fā)展，通過(guò)合理設(shè)計(jì)并合成具有理想孔隙結(jié)構(gòu)、性能更佳的催化劑已成為可能。因此，多級(jí)孔結(jié)構(gòu)的設(shè)計(jì)已成為研究熱點(diǎn)[9]，利用大孔隙改善物質(zhì)擴(kuò)散，小孔隙提供更多的反應(yīng)區(qū)域。

就多尺度研究而言，通過(guò)實(shí)驗(yàn)手段探究催化劑的孔隙結(jié)構(gòu)對(duì)反應(yīng)的影響是較為困難的，且脫離反應(yīng)器去研究催化劑性能無(wú)法直觀得到其工業(yè)價(jià)值，目前已有一些反應(yīng)器-催化劑雙尺度耦合模擬的研究。肖瑤等[10]采用幾何映射法將甲醇制丙烯（MTP）固定床反應(yīng)器與顆粒耦合，探究了催化劑粒徑對(duì)反應(yīng)過(guò)程的影響；Shi等[11]利用反應(yīng)器-顆粒耦合模型，在反應(yīng)器尺度下優(yōu)化了重油加氫脫金屬（HDM）催化劑的孔結(jié)構(gòu)；曹軍等[12]在考慮產(chǎn)物液膜傳質(zhì)阻力的情況下，構(gòu)建了能反映反應(yīng)器與催化劑顆粒雙尺度上費(fèi)托合成產(chǎn)物分布的多場(chǎng)耦合數(shù)學(xué)模型，討論了產(chǎn)物液膜及催化劑粒徑對(duì)產(chǎn)物分布的影響。

本研究運(yùn)用雙尺度耦合方法建立了甲烷化固定床反應(yīng)器的床層-催化劑顆粒雙尺度數(shù)學(xué)模型，同時(shí)考慮了床層及顆粒內(nèi)的熱量和物質(zhì)傳遞以及化學(xué)反應(yīng)過(guò)程；在一定的操作條件下，評(píng)價(jià)了催化劑顆粒的多級(jí)孔隙結(jié)構(gòu)參數(shù)（大孔平均孔徑、介孔平均孔徑、大孔孔隙率和介孔孔隙率）對(duì)顆粒中反應(yīng)及傳質(zhì)的影響，為催化劑的孔結(jié)構(gòu)設(shè)計(jì)確定了優(yōu)化區(qū)間。

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

1.1 床層-催化劑顆粒耦合計(jì)算模型

本研究建立的甲烷化反應(yīng)器模擬模型可以獲得催化劑顆粒尺度和床層尺度內(nèi)的物質(zhì)傳遞狀況、溫度變化及化學(xué)反應(yīng)過(guò)程，床層-催化劑顆粒耦合模型如圖1 所示。該模型與傳統(tǒng)均質(zhì)模型的區(qū)別在于它考慮了顆粒相的存在，因此需要用不同的方程組去定義流體相和顆粒相，通過(guò)加載顆粒表面的邊界條件將顆粒中的濃度、溫度分布與流體相的相關(guān)梯度相耦合，顆粒模型與床層模型的控制方程及邊界條件如下。

圖1 物理模型示意Fig.1 Schematic diagram of physical

顆粒模型的控制方程：

針對(duì)球形催化劑顆粒，采用Wakao 和Smith 提出的宏觀-微觀孔隙模型[13-14]，該模型用大孔平均孔徑（dM）、介孔平均孔徑（dm）、大孔孔隙率（εM）和介孔孔隙率（εm）四個(gè)參數(shù)來(lái)描述催化劑顆粒的孔隙結(jié)構(gòu)。顆粒的反應(yīng)比表面積（Sg，單位質(zhì)量催化劑的表面積）和密度（ρ）定義與孔隙結(jié)構(gòu)直接相關(guān)[15]，見(jiàn)式（6）。在該模型中，同時(shí)考慮了克努森擴(kuò)散和分子擴(kuò)散，有效擴(kuò)散系數(shù)為濃度、溫度和孔徑分布的函數(shù)[16]。采用擴(kuò)散-反應(yīng)二階偏微分方程建立顆粒內(nèi)的反應(yīng)-傳質(zhì)模型，見(jiàn)式（1），各組分?jǐn)U散僅沿顆粒半徑方向，并視顆粒內(nèi)部為恒定溫度。Solsvik 等[17]研究發(fā)現(xiàn)，顆粒內(nèi)部的黏性流和壓降很小，可以忽略不計(jì)，因此在顆粒方程中不含有對(duì)流項(xiàng)和動(dòng)量平衡。

床層模型為對(duì)流-擴(kuò)散方程，模型主要考慮床層軸向的變化。保證物料守恒的同時(shí)，床層模型還計(jì)算了沿床層的溫度變化，Qr為反應(yīng)熱，Qvd表示床層與壁面的換熱量。由于物質(zhì)穿過(guò)顆粒-流體界面進(jìn)入顆粒的流量可能是由流體主體一側(cè)對(duì)質(zhì)量傳遞的阻力決定的，因此在該模型中提出了薄膜條件假設(shè)[18]，hD表示薄膜傳質(zhì)系數(shù)，可以根據(jù)流動(dòng)中對(duì)流和擴(kuò)散質(zhì)量傳遞的相對(duì)貢獻(xiàn)的舍伍德數(shù)來(lái)計(jì)算。

1.2 甲烷化反應(yīng)動(dòng)力學(xué)模型

模型的計(jì)算要求耦合化學(xué)反應(yīng)，甲烷化反應(yīng)過(guò)程較為復(fù)雜，且伴隨多個(gè)副反應(yīng)，本研究?jī)H考慮如下三個(gè)主要的反應(yīng)：

上述反應(yīng)的動(dòng)力學(xué)方程采用Xu 等[19]建立的模型，該模型是基于Ni-MgO-Al2O3催化劑的實(shí)驗(yàn)測(cè)量得到的本征動(dòng)力學(xué)模型。

式中：下角標(biāo)1，2，3 分別代表CO 甲烷化反應(yīng)、CO2甲烷化反應(yīng)以及水煤氣變換反應(yīng)。甲烷的產(chǎn)率和選擇性定義為：

1.3 模型求解方法及驗(yàn)證

上述數(shù)學(xué)模型采用COMSOL 5.4 多物理場(chǎng)耦合仿真軟件求解。通過(guò)COMSOL 軟件中的組件耦合模塊，使得床層尺度和顆粒尺度的實(shí)時(shí)數(shù)據(jù)在兩個(gè)維度間進(jìn)行數(shù)據(jù)投影映射，從而實(shí)現(xiàn)催化劑與床層模型的耦合。計(jì)算過(guò)程中，提取床層中某點(diǎn)的各項(xiàng)物理場(chǎng)數(shù)據(jù)，以此作為該處催化劑顆粒模型的邊界條件，同時(shí)把床層模型得到的物質(zhì)濃度數(shù)據(jù)傳遞給顆粒方程，在顆粒維度對(duì)催化劑顆粒中的反應(yīng)速率、物質(zhì)濃度進(jìn)行求解，得到相應(yīng)數(shù)據(jù)；同時(shí)，在薄膜條件的假設(shè)下，顆粒內(nèi)部的物質(zhì)濃度數(shù)據(jù)經(jīng)過(guò)處理?yè)Q算為顆粒表面的物質(zhì)通量，從而將相應(yīng)數(shù)據(jù)傳遞至床層尺度，以此確定床層方程中反應(yīng)速率；各物理場(chǎng)數(shù)據(jù)通過(guò)耦合變量在兩個(gè)尺度間相互換算，兩尺度間的數(shù)據(jù)傳遞與模型求解交替進(jìn)行，直至模型收斂求得反應(yīng)器內(nèi)各物理場(chǎng)分布。

本研究為瞬態(tài)模擬，模擬時(shí)長(zhǎng)定為1 000 s。由于模型高度對(duì)稱(chēng)，只對(duì)1/8 反應(yīng)器空間進(jìn)行模擬計(jì)算。顆粒維度采用立方根序列進(jìn)行離散，單元數(shù)定為6。三維反應(yīng)器網(wǎng)格的劃分采用掃掠劃分，沿徑向作等差分布，經(jīng)網(wǎng)格無(wú)關(guān)性驗(yàn)證，最終確定網(wǎng)格數(shù)為23 860。

為檢驗(yàn)本工作所選用動(dòng)力學(xué)模型的可靠性，用上述動(dòng)力學(xué)模型計(jì)算得到的CO 甲烷化反應(yīng)和CO2甲烷化反應(yīng)的模擬數(shù)據(jù)與文獻(xiàn)數(shù)據(jù)[20-21]進(jìn)行了比較，結(jié)果如圖2 所示。由圖2 可知，模擬數(shù)據(jù)與實(shí)驗(yàn)結(jié)果基本吻合，證明模擬準(zhǔn)確可靠。

圖2 模擬值與實(shí)驗(yàn)值對(duì)比Fig.2 Comparison of simulated values and experimental data

2 結(jié)果與討論

本研究模擬的反應(yīng)器尺寸參數(shù)參照Hwang 等[21]進(jìn)行CO2甲烷化實(shí)驗(yàn)的固定床反應(yīng)器的，催化劑粒徑（dp）為5.4 mm 的球形顆粒，床層總體孔隙率（εbed）用經(jīng)驗(yàn)公式（21）[22]進(jìn)行計(jì)算，其余初始條件如表1 所示。

表1 初始條件Table 1 Initial conditions

2.1 參考模型的模擬結(jié)果

本研究需先確定一個(gè)參考模型，以此作為后期多級(jí)孔結(jié)構(gòu)優(yōu)化設(shè)計(jì)的對(duì)比基準(zhǔn)。參考模型的顆?？傮w孔隙率（εtotal）設(shè)定為0.5，顆粒大孔平均孔徑（dM）為100 nm、大孔孔隙率（εM）為0.1；顆粒介孔平均孔徑（dm）為5 nm、介孔孔隙率（εm）為0.4。計(jì)算所得反應(yīng)器出口處CO2轉(zhuǎn)化率如圖3 所示。模擬至900 s 左右時(shí)反應(yīng)器內(nèi)各個(gè)物理場(chǎng)趨于穩(wěn)態(tài)。提取1 000 s 時(shí)，反應(yīng)各組分在床層尺度的濃度分布及溫度分布如圖4 所示。

圖3 CO2 轉(zhuǎn)化率隨反應(yīng)器運(yùn)行時(shí)間變化Fig.3 Change of CO2 conversion with running time

圖4 1 000 s 時(shí)床層中各物質(zhì)濃度及溫度分布Fig.4 The concentration and temperature distributions of various substances in the bed at 1 000 s

計(jì)算得到反應(yīng)各組分在反應(yīng)器不同位置處催化劑顆粒內(nèi)部的濃度分布如圖5 所示，由床層中心入口處（R=0 mm，L=15.9 mm）顆粒內(nèi)組分濃度分布可以看出，CO2和H2進(jìn)入顆粒后迅速反應(yīng)，濃度快速下降；由于存在微量水汽變換反應(yīng)和CO 甲烷化反應(yīng)，CO 濃度先小幅升高，而后降低；在顆粒中心處，各組分濃度變化趨于平緩。床層中部（R=0 mm，L=30.9 mm）和出口端（R=0 mm，L=45.9 mm）的催化劑顆粒內(nèi)組分濃度分布與床層入口端處（R=0 mm，L=15.9 mm）相比，濃度變化漸漸平緩，顆粒內(nèi)的反應(yīng)速率明顯降低。

圖5 反應(yīng)器不同位置處顆粒內(nèi)的物質(zhì)濃度分布Fig.5 The concentration distributions of components in particles at different positions of reactor

圖6 為床層入口處催化劑顆粒內(nèi)各組分反應(yīng)速率的徑向分布。圖7 為反應(yīng)器中心處床層濃度及顆粒濃度的軸向分布?？梢钥吹剑w粒內(nèi)部與床層間存在較大的濃度梯度。反應(yīng)可能發(fā)生在催化劑顆粒近外壁面區(qū)域，約為顆?？v深20%的區(qū)域內(nèi)，催化劑中心處反應(yīng)速率較低，催化劑利用率較低，因而甲烷化反應(yīng)是受內(nèi)擴(kuò)散嚴(yán)重影響的過(guò)程。從催化劑制備的角度來(lái)說(shuō)，大孔孔隙將會(huì)降低內(nèi)擴(kuò)散阻力，使各物質(zhì)更容易在顆粒內(nèi)部進(jìn)行輸運(yùn)，從而更好地利用顆粒內(nèi)部空間，所以研究催化劑顆粒的多級(jí)孔結(jié)構(gòu)是有必要的。

圖6 顆粒內(nèi)各物質(zhì)反應(yīng)速率的徑向分布Fig.6 The reaction rate distribution of components in particles along radial direction of the reactor

圖7 反應(yīng)器內(nèi)床層濃度及顆粒濃度的軸向分布Fig.7 Axial concentration distributions in the bed and in particles at the centerline of the reactor

2.2 多級(jí)孔結(jié)構(gòu)探究

球形催化劑顆粒的宏觀-微觀孔隙模型[13-14]中，模型參數(shù)有大孔平均孔徑、介孔平均孔徑、大孔孔隙率和介孔孔隙率，本研究在設(shè)定顆?？傮w孔隙率為0.5 的情況下，考察了顆粒結(jié)構(gòu)參數(shù)對(duì)催化劑活性的影響，尋求更有效的催化劑顆粒孔結(jié)構(gòu)。

保持大孔孔徑為100 nm，介孔孔徑為5 nm，模擬了床層入口處顆粒內(nèi)CO2徑向濃度分布隨大孔孔隙率的變化情況，結(jié)果見(jiàn)圖8。由圖8 可知，大孔孔隙率的增大改善了顆粒內(nèi)部的擴(kuò)散阻力，提高了反應(yīng)程度；但是，當(dāng)大孔孔隙率大到一定程度時(shí)（εM=0.5），由于反應(yīng)比表面積減少很多，反應(yīng)程度反而會(huì)降低。不同大孔孔隙率時(shí)固定床出口的CO2轉(zhuǎn)化率、CH4選擇性和產(chǎn)率的模擬結(jié)果見(jiàn)圖9。由圖9 可知，反應(yīng)器性能隨著大孔孔隙率的增加先逐漸變好而后又急劇變差。當(dāng)大孔孔隙率為0.4 時(shí)，甲烷的產(chǎn)率達(dá)到81.1%，相較于參考模型提高了約16%。這是由于大孔孔隙率控制著大孔和介孔的分布，影響反應(yīng)的擴(kuò)散和比表面積，如果顆?？紫堵适怯奢^高比例的大孔隙組成的，擴(kuò)散速度會(huì)更快，但反應(yīng)表面活性較低。因此如果大孔孔隙率過(guò)大，反應(yīng)組分在催化劑顆粒內(nèi)的擴(kuò)散速度加快，導(dǎo)致反應(yīng)速率會(huì)先增大，但在某個(gè)值上活性表面的損失太高，轉(zhuǎn)化率將會(huì)降低。

圖8 大孔孔隙率對(duì)顆粒內(nèi)CO2 濃度分布的影響Fig.8 Effect of macroporous porosity on the distribution of CO2 concentration in particles

圖9 大孔孔隙率對(duì)反應(yīng)結(jié)果的影響Fig.9 Effect of the macroporous porosity on the reaction results

圖10 為顆粒的比表面積（Sg）和CO2有效擴(kuò)散系數(shù)（DCO2eff）隨大孔孔隙率的變化。如圖10 所示，隨大孔孔隙率的增大，比表面積呈線性減小，擴(kuò)散系數(shù)呈二次增大。因此，在活性成分沒(méi)有改變的情況下，若在造粒過(guò)程中對(duì)顆?？紫督Y(jié)構(gòu)進(jìn)行優(yōu)化，可制造出性能優(yōu)異的催化劑。

圖10 反應(yīng)比表面積和CO2 有效擴(kuò)散系數(shù)隨εM 的變化情況Fig.10 The variations of reaction specific surface area and CO2 effective diffusion coefficient with εM

在孔隙率分布固定的情況下（εM=0.4，εm=0.1），模擬研究了大孔平均孔徑和介孔平均孔徑對(duì)反應(yīng)的影響。CO2轉(zhuǎn)化率、CH4產(chǎn)率、CH4選擇性和熱點(diǎn)溫度隨大孔孔徑的變化如圖11 所示。由圖11 可知，隨著dM的增大，轉(zhuǎn)化率和選擇性有所增加，當(dāng)dM≥120 nm 時(shí)，產(chǎn)率增幅趨于平穩(wěn)。因?yàn)榭伺瓟U(kuò)散系數(shù)與孔徑為線性關(guān)系，隨著大孔孔徑的增大，有效擴(kuò)散系數(shù)變大，從而導(dǎo)致產(chǎn)率有所增加，但反應(yīng)比表面積隨dM的增大而減小，則會(huì)抑制反應(yīng)，而產(chǎn)率趨于平衡可能是由于在所探究的大孔孔徑范圍內(nèi)，擴(kuò)散的影響占主導(dǎo)，彌補(bǔ)了反應(yīng)比表面積減小的影響。

圖11 大孔孔徑對(duì)反應(yīng)結(jié)果的影響Fig.11 Effect of the macroporous pore diameter on the reaction results

介孔孔徑對(duì)反應(yīng)結(jié)果的影響如圖12 所示。由圖12 可見(jiàn)，CO2轉(zhuǎn)化率以及CH4產(chǎn)率都隨介孔孔徑的增大而降低。雖然顆粒的介孔孔徑越大，其比表面積越大[見(jiàn)式（6）]，但擴(kuò)散系數(shù)隨dm的增大而增大[見(jiàn)式（3）]，而擴(kuò)散的影響起了主導(dǎo)作用。同時(shí)，考慮到熱點(diǎn)溫度不超過(guò)積炭溫度（823 K）[23]，介孔直徑為4.5 nm 左右比較合適。

圖12 介孔孔徑對(duì)反應(yīng)結(jié)果的影響Fig.12 Effect of the mesoporous pore diameter on the reaction results

3 結(jié)論

a）利用多物理場(chǎng)耦合計(jì)算方法所建立的床層-顆粒雙尺度耦合模型能直觀地反映床層及顆粒內(nèi)的傳遞行為以及在不同尺度上的化學(xué)反應(yīng)過(guò)程。

b）床層-顆粒雙尺度耦合模擬結(jié)果表明：在未經(jīng)優(yōu)化的球型顆粒（dp=5.4 mm）中，甲烷化反應(yīng)受催化劑中內(nèi)擴(kuò)散阻力的影響，主要發(fā)生在催化劑顆粒近壁面區(qū)域，有效反應(yīng)區(qū)域約為顆粒外層的20%區(qū)域，催化劑整體利用率較??；且顆粒內(nèi)與床層間存在較大濃度梯度。

c）在給定的操作條件下，基于雙尺度耦合模擬對(duì)球形甲烷化催化劑（dp=5.4 mm，εtotal=0.5）的孔結(jié)構(gòu)探究發(fā)現(xiàn)：反應(yīng)器甲烷產(chǎn)率隨顆粒大孔孔隙率增大呈現(xiàn)先增大后減小的趨勢(shì)，大孔孔隙約為0.4 時(shí)甲烷產(chǎn)率達(dá)到最大；甲烷產(chǎn)率隨大孔平均孔徑的變化波動(dòng)較小，當(dāng)dM≥120 nm 甲烷產(chǎn)率趨于平穩(wěn)；介孔平均孔徑越小，甲烷產(chǎn)率越大，當(dāng)dm減小至4.5 nm 時(shí)產(chǎn)率變化變緩。εM和dm對(duì)反應(yīng)結(jié)果的影響較大，采用多級(jí)孔結(jié)構(gòu)催化劑有助于實(shí)現(xiàn)甲烷產(chǎn)率的提升；當(dāng)εM為0.35～0.4，dM≥120 nm，dm為4～5 nm 時(shí)，甲烷產(chǎn)率為80%左右，較參考模型提高約16%。

符號(hào)說(shuō)明