葉 銀，張亞新

新疆大學化工學院，新疆 烏魯木齊 830046

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

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

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

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

1 數(shù)學模型

1.1 床層-催化劑顆粒耦合計算模型

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

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

顆粒模型的控制方程：

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

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

1.2 甲烷化反應動力學模型

模型的計算要求耦合化學反應，甲烷化反應過程較為復雜，且伴隨多個副反應，本研究僅考慮如下三個主要的反應：

上述反應的動力學方程采用Xu 等[19]建立的模型，該模型是基于Ni-MgO-Al2O3催化劑的實驗測量得到的本征動力學模型。

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

1.3 模型求解方法及驗證

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

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

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

圖2 模擬值與實驗值對比Fig.2 Comparison of simulated values and experimental data

2 結(jié)果與討論

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

表1 初始條件Table 1 Initial conditions

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

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

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

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

計算得到反應各組分在反應器不同位置處催化劑顆粒內(nèi)部的濃度分布如圖5 所示，由床層中心入口處（R=0 mm，L=15.9 mm）顆粒內(nèi)組分濃度分布可以看出，CO2和H2進入顆粒后迅速反應，濃度快速下降；由于存在微量水汽變換反應和CO 甲烷化反應，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)的反應速率明顯降低。

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

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

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

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

2.2 多級孔結(jié)構探究

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

符號說明