基于晶格化的動力學Monte Carlo方法研究鈦酞菁配合物氧化還原反應催化能力

周 蘇, 胡 哲, 高 揚, 姜 縝

(同濟大學 a. 汽車學院； b. 中德學院， 上海 201804)

一般情況下，比較催化劑的性能，可以通過實驗獲得以下兩個參數：① 單位時間內單位活性中心所轉化的反應物分子數(Turnover Frequency，TOF)，即單位時間內在催化劑單位面積每個活性位置的轉化數；② 活性位置密度(Site Density，SD)，即催化劑單位面積上的活性位置.但是，對于未知或者具有新型結構的催化劑，很難通過實驗得到TOF值或SD值.例如，在研究一個可能的催化反應機理時，整個反應過程持續時間很短，即使超高速攝像機也只能獲得若干張宏觀變化的圖像.若想了解其微觀層級的變化并進一步推測其催化性能，并與其他催化劑進行比較，需借助其他方法或工具.

隨著計算機技術的發展，計算化學已經能夠在沒有實驗情況下給出定性或定量結論，為尋找或設計潛在的新型催化劑材料提供方向性的信息[1].例如，以鈦酞菁配合物(TiPc)為研究對象，周蘇等[2]基于密度泛函理論構建了TiPc氧化還原反應機理模型，通過計算其Gibbs自由能等指標，分析預測了TiPc的催化能力.動力學Monte Carlo方法(KMC)也是一種應用廣泛的計算化學工具[3].例如，湯星舟[4]結合動態Monte Carlo法與分子動力學研究了C面藍寶石襯底的氮化現象，分析總結了表面氮化的相關機理；Franziska Hess等[5]利用KMC研究了催化劑RuO2對于CO和HCl氧化催化反應的作用.

本文以TiPc為研究對象，基于晶格化的KMC并借助Kmos軟件，通過求解Markov鏈主方程，計算TiPc的TOF值并以其作為催化活性的參考評價指標，研究反應溫度和氫氧氣體分壓比對TiPc催化活性的影響；通過比較未知TiPc和已知Pt兩者的TOF值，預測未知催化劑TiPc的催化活性；另一方面，從計算化學工具的角度，以TiPc作為對象，旨在驗證KMC用于研究未知催化劑催化活性的可行性.

1 模型與方法

1.1 計算軟件

圖1 Kmos軟件GUI界面(氫氧反應圖)Fig.1 GUI interface of Kmos software (hydrogen and oxygen reaction)

Kmos是基于Python的一款開源軟件，其模型建立是由用戶通過Python編寫py文件作為輸入，然后由Kmos模擬給出結果.為了展示催化劑表面的微觀變化，該軟件也支持UI界面，圖1所示為氫氣和氧氣在催化劑表面某一反應時刻的截圖，圖中深色點代表氧原子，淺色點代表氫原子，白色點代表空位，基底圓形代表催化劑.

1.2 基于晶格化的KMC及應用的假設條件

基于晶格化的KMC將整個催化劑薄膜簡化為10×10個晶胞所組成的平面，同時以單個晶胞作為研究對象，以基本反應的動力學參數作為輸入，在晶胞表面采用KMC求解Markov鏈主方程[5-6,7,9]及微觀動力學模型；根據給定的反應條件(如溫度和所有反應物的分壓)，可以計算或獲取催化劑的TOF值、晶胞表面的反應物或/和生成物成分、反應鏈中的步驟、反應機理等信息[10-11].針對所研究的對象，應用晶格化的KMC時，通常給出如下4個方面的假設.

1.2.1特殊位置的吸附和晶格排布 假設所有反應物和中間物都吸附在特殊位置，即晶體表面的活性位置.在實際的KMC模擬時，只考慮這一個晶格，整個晶體表面可以看成由這個晶格重復排列組合而成.例如，低指數單晶表面中一個晶格是由多個相同的晶胞組成.圖2所示為一個面心立方排列(100)的催化劑模型，它有2個活性位置[5].

圖2 催化劑表面CO氧化反應(包含吸附、擴散、反應和解吸附過程)模型示意圖Fig.2 Catalyst surface diagram of CO oxidation reaction on catalyst surface (including adsorption, diffusion, reaction and desorption process)

1.2.2Rare-event動力學和Markov鏈主方程 假設時間演化是基于Rare-event動力學的.利用時間尺度分離，根據微觀動力學理論，假設事件的發生是相互獨立的，可以基于Markov近似來推演反應過程.因此，事件(從一種組態到另一種組態)的發生可以通過Markov鏈主方程描述：

(1)

1.2.3基本反應的確定 假設前文提到的一個晶格中組態變化只考慮1～10個位置.因此，可以根據局部產物Ea的晶格狀態和變化后的局部產物Pa的晶格狀態以及相應的伴隨速率Ka來確定一個基本反應，即

(2)

1.2.4狀態轉移矩陣的大小和結構 式(1)中w可分解成一系列基本反應：

(3)

(4)

式中：ka為基本反應a的反應速率.

以上假設用于簡化Markov鏈主方程的求解.

1.3 模型

1.3.1催化反應速率計算 KMC模擬的主要特點之一就是可以計算各種催化劑的反應速率，包括TOF值.如果基本反應只有一個最終生成物，那么TOF的計算只考慮這一種生成物.

對于基本反應a，在任意時刻t，

TOFa=

(5)

1.3.2參數計算 基于密度泛函理論(DFT)，利用Amsterdan Density Functional(ADF)等軟件，根據構造的TiPc在氧化還原過程中形成的各粒子能量變化，可計算TiPc的氧吸附能：

EbO2=EML-O2-(EO2+EML)

(6)

式中：EML-O2為氧分子和吸附粒子的總能量；EO2為孤立的氧分子的總能量；EML為吸附粒子的能量.此式也可用于計算Pt的氧吸附能.

同理，計算TiPc的氫吸附能：

EbH2=EML-H2-(EH2+EML)

(7)

式中：EML-H2為氫分子和吸附粒子的總能量；EH2為孤立的氫分子的總能量；EML為吸附粒子的能量.此式也用于計算Pt的氫吸附能.

2 結果與討論

2.1 TiPc的TOF值

在模擬TiPc表面氫氣和氧氣反應時，反應溫度T和TiPc表面的氫氣分壓與氧氣分壓的比值pH2/pO2是兩個重要的參數.在Kmos軟件運行環境下，取T在470～570 K(間隔5 K)之間的21個值和pH2/pO2在 0.1～10之間的40個值，可以計算得到相應的TOF值.以溫度為470、520和570 K的情況為例，分析反應溫度T和氫氧氣體分壓比值pH2/pO2對TiPc催化劑TOF值的影響.

2.1.1TOF值隨反應溫度T的變化 由化學反應基本原理可知，反應溫度T是影響催化劑活性(TOF值)的重要因素.溫度較低時，催化劑活性較小，對于反應促進效果較差；隨著溫度上升，催化效果逐漸增強.但是，大多數催化劑都有其催化活性適宜溫度范圍.溫度過高，會使催化劑發生燒結，從而破壞其催化活性.為了分析反應溫度對TOF值的影響，圖3所示為3種溫度條件下不同氫氧氣體分壓比值對應的TOF值(實線部分).與溫度條件520 K及570 K相比，溫度470 K下的TOF值過低，所以圖3的虛線部分也給出了lg(TOF)值，以便于對比分析.

圖3 不同溫度下pH2/pO2對應的TiPc催化劑TOF及lg(TOF)值Fig.3 TOF and lg(TOF) values of TiPc under different pH2/pO2 at 470, 520, 570 K

由圖3可明顯看出：① 隨著反應溫度上升，同一氣體分壓比值對應的TOF值也依次增高，這符合催化劑基本特性.但反應溫度超過催化劑適宜反應溫度區間后，催化劑催化性能下降.Choi等[14]及陳維民等[15]通過各自的研究發現，催化劑會隨著溫度上升而發生聚集，超過一定程度會降低催化劑催化活性；② 與溫度區間520～570 K相比，在溫度區間470～520 K內TOF值的增長幅度較小，這說明TiPc的活性雖然與溫度正相關，但并不是線性相關；③ lg(TOF)值變化隨溫度變化近似呈線性關系.

2.1.2TOF值隨氫氧氣體分壓比值pH2/pO2的變化 除反應溫度T外，氫氧氣體分壓比值pH2/pO2同樣會影響催化劑活性(即TOF值).為了研究相關規律，將3種溫度條件下TOF值(圖3實線部分)在橫坐標 0.1～0.33 區間內的曲線放大得到圖4.

圖4 不同溫度下pH2/pO2在 0.1～0.33 區間時對應的TiPc催化劑TOF值Fig.4 TOF values of TiPc within pH2/pO2 0.1-0.33 at 470, 520, 570 K

圖5 溫度470至570 K條件下不同氧氫氣體分壓比值對應的TiPc催化劑lg(TOF)值Fig.5 lg(TOF) values of TiPc under different O2 and H2 partial pressure ratios between 470 K and 570 K

結合圖3和4的信息可知，3種不同溫度條件下TOF值的最大值出現在氫氧氣體分壓比值0.289～0.326區間內.在氫氧氣體分壓比值0.1～0.2之間，TOF值單調緩慢上升；在氫氧氣體分壓比值0.2～2區間內，TOF值先急速上升至最高點，然后單調急速下降；從氫氧氣體分壓比2開始，TOF值單調平穩下降至接近0值.催化反應主要經歷以下3個步驟[12-13]：① 反應物被催化劑吸附在表面；② 反應物在催化劑表面進行反應；③ 反應生成物從催化劑表面脫附.其中，步驟①對應前文所計算吸附能，步驟③保證反應生成物能夠及時富集及保證催化劑再生.根據 1.2 中的假設，所有的反應物和中間物都吸附在表層催化活性位置.因此，根據化學平衡移動原理，當氫氧分壓比從零開始升高時，催化劑表面反應物氫氣的濃度增加，反應傾向于向正反應方向移動，對于催化劑表面的活化中心具有很好的活化作用.但是，氫氧分壓比進一步增大，使得氫氣的覆蓋度過大，會出現“氫中毒”現象，反而使催化劑活性下降，抑制反應的進行.本文所用方法不能從定量上解釋TOF值急速變化的原因，但是，從上述的分析中可以推測，TiPc催化劑的活性區域在氫氧分壓比值 0.2 ～2之間，此區間內氫氣和氧氣在催化劑表面以較為合適正向反應的比例共存.

2.1.3采用2D Map圖分析TOF 為了能夠同時研究反應氫氧氣體分壓比及反應溫度對TiPc催化活性或反應速率的影響，利用線性插值，將上述21組數據繪制成2D Map圖；為了分析方便，取氫氧氣體分壓比值的倒數即氧氫氣體分壓比值為橫坐標，縱坐標為反應溫度(見圖5)；二維圖中區域的顏色深淺則對應TOF值(反應速率)，顏色越淺代表TOF值越大.從TiPc的計算數據及上文的結論可知，在氫氧氣體分壓比值 0.1～0.2 區間內，TiPc的TOF值非常小，接近于零.為了增強這種區別，取TOF值的常用對數.圖5所示為反應溫度 470～570 K 下不同氧氫氣體氣體分壓比值對應的 lg(TOF) 值.

圖5中暗部區域對應的橫坐標范圍為氧氫氣體分壓比值 3.8～10，對應于氫氧氣體分壓比值 0.1～0.26，基本符合 2.1.1 和 2.1.2 得出的結論，該區域內TiPc反應活性遠低于左側亮部區域，且亮部區域和暗部區域突變位置即為TOF值陡降陡升位置.

2.2 Pt的TOF值

采用上述的相同條件，在Kmos軟件運行環境下計算得到相應的Pt催化劑TOF值，用于分析反應溫度T及氫氧氣體分壓比值pH2/pO2對TOF值的影響，從而進一步比較TiPc與Pt兩者的催化活性.

2.2.1TOF值隨反應溫度T的變化 由Pt催化劑的計算數據可知，Pt催化劑TOF值在氫氧氣體分壓比值 0.8～10 區間內幾乎不變，為了便于分析，只選用氫氧氣體分壓比值 0～0.8 區間內的計算數據.圖6所示為3種溫度條件下不同氫氧氣體分壓比值對應的TOF值(實線部分)，且與TiPc催化劑情況相同，為突出470 K與520 K溫度下TOF值曲線間的差別，圖6的虛線部分也給出了 lg(TOF)值.

圖6 不同溫度下pH2/pO2對應的Pt催化劑TOF及lg(TOF)值Fig.6 TOF and lg(TOF) values of Pt under pH2/pO2 at 470, 520, 570 K

由圖6可看出，不同溫度條件下，Pt催化劑表現出與TiPc催化劑相似的特性：① 隨著反應溫度上升，同一氣體分壓比值對應的TOF值也依次增高；② 雖然TOF值與溫度不是線性正相關，但是 lg(TOF) 值隨溫度變化近似呈線性關系.

2.2.2TOF值隨氫氧氣體分壓比值pH2/pO2的變化 為了研究氫氧氣體分壓比值pH2/pO2對催化劑活性(即TOF值)的影響關系，將3種溫度條件下TOF值(圖6實線部分)在橫坐標 0.1～0.26 區間內的曲線放大得到圖7.

由圖7可以看出，3種不同溫度條件下TOF值的最大值均出現在氫氧氣體分壓比值0.229處，此時Pt催化劑的催化活性最高；與TiPc催化劑相比，Pt催化劑TOF值的最大值明顯小幾個數量級，所以也可以預測TiPc催化劑的活性高于Pt催化劑.

圖7 不同溫度下pH2/pO2 在 0.1～0.26 區間時對應的Pt催化劑TOF值Fig.7 TOF values of Pt within pH2/pO2 0.1-0.26 at 470, 520, 570 K

2.2.3采用2D Map圖分析TOF 將Pt催化劑21組原始數據利用線性插值繪制成2D Map圖.由圖6實線部分可知，在氫氧氣體分壓比值 0.1～0.3 區間內，Pt催化劑TOF值遠大于其他區域.為了增強這種區別，取TOF值的常用對數.圖8所示為反應溫度470～570 K下不同氧氫氣體分壓比值對應的lg(TOF)值.

圖8中亮部區域對應的橫坐標范圍為氧氫氣體分壓比值 3.9～10，對應氫氧氣體分壓比值 0.1～0.256，基本符合本節 2.2.1 和 2.2.2 所得結論，該區域內Pt催化劑反應活性遠高于左側暗部區域，且亮部區域和暗部區域突變位置即為TOF值陡降陡升位置.

圖8 470～570 K條件下不同氧氫氣體分壓比值對應的Pt催化劑lg(TOF)值Fig.8 lg(TOF) values of Pt under different O2 and H2 partial pressure ratios between 470 K and 570 K

2.3 分析與討論

TiPc催化劑TOF值與反應溫度T及氫氧氣體分壓比值pH2/pO2存在一定復雜關系.TiPc催化劑TOF值會隨溫度變化而發生改變，隨著反應溫度由470 K增長到570 K，TOF值同樣呈現出遞增的趨勢；氫氧氣體分壓比值同樣會影響催化劑活性(TOF值)，TiPc催化劑TOF值的最大值出現在氫氧氣體分壓比值 0.289～0.326 區間內，氫氧氣體分壓比值在 0.1 至該區間，TOF值首先接近于零，之后單調遞增；氫氧氣體分壓比值在該區間至10，TOF值單調遞減直至接近0值.此外，TiPc催化劑TOF值的最大值出現在反應溫度570 K和氫氧氣體分壓比值 0.326 處，其值為 2.034×107.

與TiPc催化劑不同，Pt催化劑TOF值在氫氧氣體分壓比值 0.1～0.229 之間單調遞增，在氫氧氣體分壓比值 0.229～0.26 之間單調遞減，可以推測氫氧氣體分壓比值 0.1～0.26 區間是Pt催化劑的活性區域；當氫氧氣體分壓比值處在 0.26～10區間內，TOF值接近于零且單調遞減.Pt催化劑TOF值的最大值出現在反應溫度570 K和氫氧氣體分壓比值 0.229 處，其值約為82.

與Pt催化劑相比，TiPc催化劑TOF值的最大值以及TiPc催化劑的2D Map圖中亮部區域對應氫氧氣體分壓比值區間都明顯更大，可以推測TiPc催化劑的活性區域更廣且其催化活性更高.

3 結語

基于晶格化的動力學Monte Carlo方法，采用TOF值作為參考評價指標，評價TiPc催化劑在反應中的催化活性，模擬氣-固界面氫、氧分子在催化劑表面的反應.討論了TiPc在不同反應溫度及不同氫氧氣體分壓比值的反應條件下的催化活性.TiPc催化劑TOF值隨反應溫度上升而遞增，同時也隨氫氧氣體分壓比值上升而首先遞增之后遞減.TiPc催化劑TOF值的最大值出現在反應溫度為570 K和氫氧氣體分壓比值 0.326 處.采用與TiPc催化劑相同的反應條件，在Kmos軟件運行環境下計算得到相應的Pt催化劑TOF值，并與TiPc催化劑TOF值進行比較，得到了與文獻[2]相同的預測結果，即TiPc催化劑催化活性強于Pt催化劑.

因計算方法和計算能力的局限性，本文僅研究了TiPc在相對較簡單的氣-固界面上反應的催化性能，沒有涉及質子交換膜燃料電池中催化劑在固-液界面的真實反應環境.