Nd摻雜CeO2（111）表面穩定結構的熱力學分析實驗設計

朱后禹，劉東源，欒 波，趙 文，郭文躍

（1.中國石油大學（華東）材料科學與工程學院，山東青島 266580；2.山東京博控股集團有限公司，山東濱州 262500）

隨著科學技術和計算機的發展，科學研究的體系越來越復雜，傳統的解析推導方法已不復應用，而計算機運算能力不斷提高，為復雜體系的研究提供了新的科學手段—計算材料科學，并迅速得到發展。通過計算機模擬，可以初步了解不同溫度、壓強條件下材料的結構穩定性和使用性能等，基于所計算的材料性能演變規律和機理，可以為進一步改善材料性能提供理論指導。目前基于密度泛函理論（Density functional theory，DFT）[1]的Vienna ab initio simulation package （VASP），Material Studio 等軟件所做的結構優化都是在0 K 溫度下求解Kohn-Sham 方程[2]（即電子密度泛函理論中的單電子方程）來求得體系的基態電荷密度，從而達到收斂標準。在很多應用上，比如準確地計算開殼層體系、具有高自旋多重度的體系、激發態體系等，還存在一定的局限性。因此，通過結構優化而得到的固體結構模型通常是滿足化學計量數的無缺陷模型，而不是應用中某個溫度壓強下所處的一個實際狀態。而用第一性原理計算熱力學性質主要目的是預測固體材料的結構模型在某一溫度和壓強下處于一個什么樣的實際狀態。根據DFT計算得到結構模型的能量和熵，求出兩個表面結構反應變化前后的吉布斯自由能變，再結合實際應用條件給定溫度和壓強數值，來判斷這個結構處于什么樣的相態。實驗采用了基于密度泛函理論的第一性原理計算方法，依托于超算中心集群資源，通過模擬計算預測了摻雜CeO2材料在一定的溫度和壓強范圍內表面氧空位的自發形成情況，從而確定實際反應條件下摻雜CeO2表面的穩定結構。

1 實驗設計依據

CeO2是一種重要的稀土氧化物材料，在固體氧化物燃料電池、尾氣處理和生物醫藥等方面都有重要的實際工業應用。大多數關于CeO2材料的研究都是利用其優越的氧化還原能力來做表面催化反應，具體表現在CeO2表面比較容易發生晶格氧原子的釋放和遷移。在某一溫度下，由于熱漲落效應，氧原子會脫離晶格結構形成空位點缺陷，即氧空位。脫離的氧原子通常可以兩兩結合成為O2分子，或者直接參與表面催化氧化反應氧化吸附在表面的反應物。另外，晶格氧原子脫離形成一個氧空位后會在表面釋放2個電子，由于Ce-4f 軌道的強局域性，電子會將近鄰的Ce4+離子還原為Ce3+。由此可知，氧空位形成的數量具有飽和性，隨著氧空位形成數量的增多，后續氧空位的形成越困難，CeO2對表面吸附物種的氧化能力會降低。因此，在做有關CeO2表面催化的模擬研究之前，有必要結合熱力學計算預測實際反應條件下CeO2表面的氧空位形成情況。

近幾年的研究都致力于對純CeO2表面進行金屬團簇修飾或金屬原子摻雜來降低氧空位形成能，加速氧空位形成速率，從而提高CeO2表面的催化氧化反應速率。而實際應用較廣的摻雜CeO2的氧空位形成能可調控至1～2eV。因此在中低溫（300～800K）條件下，CeO2表面自發形成氧空位的可能性較大。在進行CeO2表面催化的模擬研究時，直接使用無缺陷的完整CeO2表面是不嚴謹的。本實驗應用VASP 軟件，以稀土金屬元素Nd 摻雜CeO2表面為例，基于周期性密度泛函理論對Nd 摻雜CeO2進行結構優化，并通過熱力學計算預測在高溫還原性氣氛條件下的表面氧空位形成情況，最終繪制Nd 摻雜CeO2表面穩定結構相圖。

2 模型構建及參數設置

采用Vienna Ab-into Simulation Package（VASP）計算軟件，并基于超算中心集群開展計算工作。

2.1 模型構建可分為如下幾個步驟

（1）從Material Studio 軟件自帶的結構庫中導入CeO2的晶胞結構；

（2）利用建模工具進行切面操作得到CeO2（111）面，建立超晶胞，設置真空層高度；

（3）用Nd 原子替換CeO2（111）表面的一個Ce 原子，期間需充分考慮Nd 原子在表面上各種可能的替位摻雜位置，并通過結構優化比較不同摻雜結構的總能量，確定最穩定的Nd 摻雜CeO2（111）表面結構。

如圖1所示，Nd 摻雜CeO2（111）表面模型為六層原子組成，每層含有9個原子（3×3）的超晶胞模型。在表面垂直方向添加了15nm 的真空層，以防止在該方向上出現周期性結構的相互作用。在模擬過程中，對Nd 摻雜CeO2（111）表面模型的上三層原子進行無對稱約束的結構優化，并固定下三層原子于計算的體相點陣位置，以模擬固體內部幾乎不受表面弛豫影響的情況。

圖1 Nd摻雜CeO2（111）表面模型的俯視圖（左）和側視圖（右）。為了加以區分，最頂層的兩層原子用大號球棍模型顯示。紅色和淺黃色表示O原子和Ce原子，其中一個Ce原子被Nd原子取代并以綠色標示。圖標O1T-O3T和O1S-O3S分別表示表層和次表層可能形成的氧空位位置。

2.2 計算參數設置如下

所有計算都是在以DFT 理論為基礎的VASP 軟件中進行的。離子-電子間相互作用采用平面波增強方法（PAW）處理，電子間交換關聯勢能采用廣義梯度近似（GGA）處理。在靜態自洽計算中，Ce（5s，5p，6s，4f，5d），Nd（5s，5p，6s，4f，5d），和O（2s，2p）作為價電子處理，其余電子與原子核視為離子并用PBE 贗勢代替。在所有結構優化中設置截斷能為400eV，離子位移收斂標準為0.03eV/?，電子能量收斂標準為10-6eV。所有計算均采用自旋極化。布里淵區內的積分采用3×3×1的Monkhorst-Pack k 點設置。采用DFT+U 方法處理Ce-4f 軌道的強關聯相互作用，其中Hubbard U修正數值U-J設置為5eV。

氧空位形成能（Evac）計算公式是：

其中，En·Vo-surf，EO2和Efull-surf分別表示生成n個氧空位后的Nd 摻雜CeO2（111）表面的總能量，真空中一個氧氣分子的能量以及無缺陷表面的總能量。從頭算熱力學的基本原理是把Nd 摻雜CeO2（111）表面生成一個氧空位的過程等價于一個固體表面被還原性氣體還原的具體反應過程，從而與溫度和還原性氣體的分壓建立函數關系，計算該氧化還原反應的吉布斯自由能變。具體反應表達式可用如下式子表示

其中n表示在一個反應周期內有多少H2分子被Nd 摻雜CeO2（111）表面氧化，同時生成了n個氧空位以及最終產物H2O。則該反應過程的吉布斯自由能變可用如下式子表示：

GH2OandGH2表示氣相H2O 分子和H2分子的吉布斯自由能。單位摩爾的氣相分子H2O 或H2的吉布斯自由能可進一步由兩者的化學勢μH2O(T,P)和μH2(T,P)表示：

式中，下標Y表示氣相分子H2O 或H2，R為熱力學理想氣體常數8.314J/（mol·K），標準熱力學數據焓HY(T,Pθ)和熵SY(T,Pθ)可從NIST 官方數據庫中查詢。結合式（3）、（4）可知，當ΔG=0，n=1時，可得到壓強pY與溫度T的函數關系曲線，此曲線的物理意義即為無缺陷的Nd 摻雜CeO2（111）表面和生成一個最近鄰氧空位表面的臨界狀態。以此類推，代入n=2,3…時可得到形成后續氧空位的臨界相態曲線。以壓強pY與溫度T為坐標，整合n=1,2,3…代入后得到的臨界相態曲線即可。

3 實驗結果及討論

純CeO2（111）表面的表層氧空位形成能約為2.7eV，比貴金屬氧化物更易于形成氧空位，但仍不能滿足固體氧化物燃料電池啟動溫度低溫化的需求，所以本研究中也采用當前研究較為廣泛的稀土金屬Nd 摻雜CeO2模型，來進一步調控氧空位形成能。圖2為結構優化并最終達到收斂標準的Nd 摻雜CeO2（111）表面模型，其中一個表層金屬Ce 原子替換為Nd原子。首先，由于摻雜原子引入，導致了一定程度上的晶格畸變。相比于純CeO2（111）表面中的Ce-OT，Ce-OS 鍵長2.355?和2.356?，Nd 原子與表層晶格氧原子OT之間的鍵長拉伸至為2.362?，與次表層氧原子OS之間的鍵長縮短為2.348?。Nd摻雜破壞了晶格結構原有的對稱性，利于表面氧空位的形成。

圖2 Nd摻雜CeO2（111）表面的優化結構及相關鍵長參數

根據周圍晶格氧原子與摻雜原子距離的不同，需要考慮6個可能的氧空位形成位置。計算其氧空位形成能（表1）后發現，相對于摻雜原子，最近鄰氧原子脫離晶格形成氧空位所需能量最低，只有2.253eV。而最近鄰次表層氧空位形成能為2.483eV。這一結果與上文討論的鍵長數據相吻合，氧空位形成能的大小與金屬-氧原子鍵長呈正相關。另外，從整體規律上來看，距離摻雜原子越遠，氧空位形成能越大，趨向于純CeO2（111）表面的2.7eV。由表1所示數據推斷，在催化氧化表面分子的反應過程中，與摻雜原子最近鄰的表層氧原子會優先被釋放，形成氧空位。

表1 Nd摻雜CeO2（111）表面不同位置的氧空位形成能

如上文所述，從NIST 查得的標準熱力學數據焓HY(T,Pθ)和熵SY(T,Pθ)是在某些特定溫度和標準壓強下測定的離散數值，不能滿足表面結構相圖所需的連續曲線關系。但根據理想氣體的熱力學關系，可用Origin 進行數據擬合使之成為關于溫度的連續函數H（T）=E+nRT和S（T）=CplnT-Rlnp+S0。擬合結果如圖3所示。

圖3 實驗反應溫度范圍內氣相H2分子的標準熱力學數據焓HH2(T,Pθ)和熵SH2(T,Pθ)的擬合曲線

因擴散速率通常比化學催化反應速率快幾個數量級，所以當固體氧化物燃料電池陽極催化反應經過一段時間達到平衡時，最終產物H2O 的生成速率也會迅速趨于穩定，另外，燃料電池運行過程中產生的燃料廢氣會及時向外界排出，所以氣相中H2O 的分壓也會迅速趨于一個定值，具體可參考實驗中的氣相組分測量值。而通入的反應氣體H2可以人為調控，所以H2的氣體分壓是氣相環境的決定因素。按照3.2節中所述方法代入數據處理，可得到以H2分壓PH2和溫度T 為橫縱坐標的熱力學穩定結構相圖，如圖4所示。在固體氧化物燃料電池陽極的實際反應條件下（黑色虛線標示，PH2=1atm，T=800K），一個氧原子會脫離晶格形成氧空位，即Nd 摻雜CeO2（111）表面可以提供一個晶格O 原子催化氧化H2分子。由此不難推斷，在上述反應條件下，體系的初始結構為H2吸附在無缺陷的Nd 摻雜CeO2（111）表面。經過一系列基元反應之后，最終H2被Nd 摻雜CeO2（111）表面氧化為H2O 分子并在表層留下一個氧空位。即確定了Nd 摻雜CeO2（111）表面催化氧化H2反應這一體系的初末態，為后期研究Nd 摻雜CeO2（111）表面的H2催化氧化反應機理奠定了較為準確的模型基礎。此外，也說明了Nd 摻雜CeO2（111）表面在此條件下是穩定的晶格結構，晶格中的Nd 摻雜原子仍然保持較為穩定的6個O 原子配位環境。只有當溫度過高，反應氣體中H2分壓較大時，才可能使表面形成3個以上的表層氧空位，金屬離子配位數急劇降低，從而發生原子聚集和Ostwald 效應，破壞表面的穩定性。

圖4 反應條件下Nd摻雜CeO2（111）表面相圖

不同的顏色區域表示在該溫度和壓強環境下的穩定結構。固體氧化物燃料電池陽極的反應條件用虛線標示（PH2=1atm，T=800K）。

4 結論

本實驗設計依托超算中心計算平臺，采用VASP 軟件，并結合從頭算熱力學方法研究了反應條件下Nd 摻雜CeO2（111）的表面相態，從微觀角度研究了摻雜CeO2表面結構隨溫度和壓強的變化，確定了反應條件下的Nd 摻雜CeO2表面穩定結構。結合從頭算熱力學繪制表面結構相圖的方法可以提供實驗條件下表面結構的熱穩定性，并推測不同溫度和壓力條件下的表面狀態及構型。