Nd摻雜CeO2（111）表面穩(wěn)定結(jié)構(gòu)的熱力學(xué)分析實(shí)驗(yàn)設(shè)計(jì)

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

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

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

1 實(shí)驗(yàn)設(shè)計(jì)依據(jù)

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

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

2 模型構(gòu)建及參數(shù)設(shè)置

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

2.1 模型構(gòu)建可分為如下幾個步驟

（1）從Material Studio 軟件自帶的結(jié)構(gòu)庫中導(dǎo)入CeO2的晶胞結(jié)構(gòu)；

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

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

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

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

2.2 計(jì)算參數(shù)設(shè)置如下

所有計(jì)算都是在以DFT 理論為基礎(chǔ)的VASP 軟件中進(jìn)行的。離子-電子間相互作用采用平面波增強(qiáng)方法（PAW）處理，電子間交換關(guān)聯(lián)勢能采用廣義梯度近似（GGA）處理。在靜態(tài)自洽計(jì)算中，Ce（5s，5p，6s，4f，5d），Nd（5s，5p，6s，4f，5d），和O（2s，2p）作為價(jià)電子處理，其余電子與原子核視為離子并用PBE 贗勢代替。在所有結(jié)構(gòu)優(yōu)化中設(shè)置截?cái)嗄転?00eV，離子位移收斂標(biāo)準(zhǔn)為0.03eV/?，電子能量收斂標(biāo)準(zhǔn)為10-6eV。所有計(jì)算均采用自旋極化。布里淵區(qū)內(nèi)的積分采用3×3×1的Monkhorst-Pack k 點(diǎn)設(shè)置。采用DFT+U 方法處理Ce-4f 軌道的強(qiáng)關(guān)聯(lián)相互作用，其中Hubbard U修正數(shù)值U-J設(shè)置為5eV。

氧空位形成能（Evac）計(jì)算公式是：

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

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

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

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

3 實(shí)驗(yàn)結(jié)果及討論

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

圖2 Nd摻雜CeO2（111）表面的優(yōu)化結(jié)構(gòu)及相關(guān)鍵長參數(shù)

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

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

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

圖3 實(shí)驗(yàn)反應(yīng)溫度范圍內(nèi)氣相H2分子的標(biāo)準(zhǔn)熱力學(xué)數(shù)據(jù)焓HH2(T,Pθ)和熵SH2(T,Pθ)的擬合曲線

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

圖4 反應(yīng)條件下Nd摻雜CeO2（111）表面相圖

不同的顏色區(qū)域表示在該溫度和壓強(qiáng)環(huán)境下的穩(wěn)定結(jié)構(gòu)。固體氧化物燃料電池陽極的反應(yīng)條件用虛線標(biāo)示（PH2=1atm，T=800K）。

4 結(jié)論

本實(shí)驗(yàn)設(shè)計(jì)依托超算中心計(jì)算平臺，采用VASP 軟件，并結(jié)合從頭算熱力學(xué)方法研究了反應(yīng)條件下Nd 摻雜CeO2（111）的表面相態(tài)，從微觀角度研究了摻雜CeO2表面結(jié)構(gòu)隨溫度和壓強(qiáng)的變化，確定了反應(yīng)條件下的Nd 摻雜CeO2表面穩(wěn)定結(jié)構(gòu)。結(jié)合從頭算熱力學(xué)繪制表面結(jié)構(gòu)相圖的方法可以提供實(shí)驗(yàn)條件下表面結(jié)構(gòu)的熱穩(wěn)定性，并推測不同溫度和壓力條件下的表面狀態(tài)及構(gòu)型。