氧化釓摻雜氧化鈰的分子動(dòng)力學(xué)模擬

陳云俊，孫毅，劉一志，康高英

(哈爾濱工業(yè)大學(xué)航天科學(xué)與力學(xué)系，黑龍江哈爾濱150001)

目前，燃料電池作為一種清潔能源越來越受到人們的關(guān)注[1].氧化釓摻雜的氧化鈰(GDC)是目前應(yīng)用較為廣泛的固體氧化物燃料電池(SOFC)的電解質(zhì)之一，GDC可以在中低溫(500～700 K)下工作，且溫度較低時(shí)其電導(dǎo)率比傳統(tǒng)的YSZ(氧化釔穩(wěn)定的氧化鋯)電解質(zhì)高一個(gè)數(shù)量級(jí)以上[2].由于工作環(huán)境及材料尺度的限制，當(dāng)前很難用傳統(tǒng)的實(shí)驗(yàn)手段清晰地揭示GDC在各種工作環(huán)境下的體系變化.分子動(dòng)力學(xué)(MD)模擬作為一種在原子尺度下研究材料性能演化規(guī)律的技術(shù)，被廣泛用于GDC的性能研究.目前國際上關(guān)于MD模擬GDC的相關(guān)工作已經(jīng)取得很多成就[3-5].雖然在MD模擬過程中，積分算法、時(shí)間步、初始構(gòu)型等因素都能夠影響模擬結(jié)果，但作為描述模型中粒子間相互作用的函數(shù)，勢(shì)函數(shù)才是MD模擬成功的核心因素.當(dāng)前，采用不同勢(shì)參數(shù)對(duì)GDC進(jìn)行系統(tǒng)MD模擬的工作相對(duì)較少，此處總結(jié)了當(dāng)前使用較多的幾組勢(shì)參數(shù)，同時(shí)給出了一套新的勢(shì)函數(shù)，并對(duì)GDC模型進(jìn)行了系統(tǒng)的模擬，詳細(xì)解釋了模擬過程中體系發(fā)生的變化及其導(dǎo)電機(jī)理.

1 GDC勢(shì)函數(shù)

GDC是摻雜了氧化釓的氧化鈰固溶體，化學(xué)式為 CexGd1-xO2-0.5x(x 為 Gd3+的摩爾濃度)，其分子動(dòng)力學(xué)模擬是建立在Born固體模型的基礎(chǔ)上，陽離子之間只有庫倫作用.一般采用 Born-Meyer-Buckingham函數(shù)來描述GDC間粒子的相互作用，具體表達(dá)式為[6]

式中:qi、qj分別為離子 i與 j的電荷量;rij為 i，j間對(duì)應(yīng)的離子間距;ε0為介電常數(shù);A、ρ、C為對(duì)應(yīng)參數(shù)，表1列出了目前較常用的幾組參數(shù).其中，Ⅰ:Lewis提出的勢(shì)參數(shù)[6];Ⅱ:Gotte提出的勢(shì)參數(shù)[7];Ⅲ:Licia提出的勢(shì)參數(shù)[8].

表1 不同勢(shì)參數(shù)列表Table 1 Different potential parameters

表1中的勢(shì)函數(shù)及對(duì)應(yīng)參數(shù)主要基于經(jīng)驗(yàn)勢(shì)方法，其參數(shù)的獲得十分依賴實(shí)驗(yàn)值，應(yīng)用范圍相對(duì)狹窄.在此，依據(jù)量子物理和量子化學(xué)方法，提出了一種離子晶體勢(shì)函數(shù)的求取方法，應(yīng)用范圍較廣.勢(shì)函數(shù)形式與式(1)有所區(qū)別，具體如下[9]:

式中:Φ+-描述了陰陽離子之間的作用，Φ++描述了陽離子之間的作用，Φ--描述了陰離子之間的作用.a、b、c為對(duì)應(yīng)的勢(shì)參數(shù)，具體如表2.

表2 新勢(shì)函數(shù)及對(duì)應(yīng)參數(shù)Table 2 New potential and corresponding parameters

2 分子動(dòng)力學(xué)模擬

2.1 晶格常數(shù)

首先計(jì)算了298 K時(shí)的晶格常數(shù).模擬過程中，采用NPT系統(tǒng)，時(shí)間步為1 fs.先根據(jù)晶體學(xué)知識(shí)建立尺寸為5l×5l×5l(其中l(wèi)為設(shè)定的初始晶格常數(shù)，值為5.41?)的模型，再對(duì)整個(gè)模型在298 K溫度下弛豫100萬步，模型完全穩(wěn)定后，取后50萬步的平均值作為計(jì)算結(jié)果.

圖1給出了不同摻雜濃度(文中摻雜濃度都的摩爾濃度)下GDC的晶格常數(shù)變化情況.由圖1可知，除了Gotte參數(shù)的計(jì)算結(jié)果，其他4組勢(shì)參數(shù)得到的結(jié)果與實(shí)驗(yàn)值[10]變化趨勢(shì)一致.實(shí)際上只給出了CeO2的勢(shì)參數(shù)，并未給出Gd2O3的參數(shù)[7]，為了對(duì)比不同參數(shù)的計(jì)算效果，采用Gotte的CeO2參數(shù)結(jié)合其他文獻(xiàn)給出的Gd2O3參數(shù)計(jì)算了GDC晶格常數(shù).

圖1 模擬晶格常數(shù)與實(shí)驗(yàn)值的對(duì)比Fig.1 Calculated lattice constants from different potentials and experimental data

對(duì)于沒有摻雜的CeO2材料，粒子間通過勢(shì)函數(shù)相互作用，整個(gè)體系處于穩(wěn)定狀態(tài).Gd2O3的摻雜引入了氧空位，體系中粒子受力平衡被破壞.為了達(dá)到新的平衡狀態(tài)，粒子的位置發(fā)生了改變，整個(gè)體系產(chǎn)生膨脹，從而導(dǎo)致晶格常數(shù)增加.圖1中4組勢(shì)參數(shù)得到的結(jié)果對(duì)比，可以看出Lewis的結(jié)果有較大區(qū)別，這是因?yàn)樵摻M參數(shù)提出于20世紀(jì)80年代，沒有充足的實(shí)驗(yàn)值對(duì)參數(shù)進(jìn)行驗(yàn)證，因此該套參數(shù)的計(jì)算結(jié)果與實(shí)驗(yàn)值差別較大，只能用于定性說明體系的變化情況;而Hideaki及Licia的勢(shì)參數(shù)是基于大量實(shí)驗(yàn)數(shù)據(jù)基礎(chǔ)上擬合得到的，因此這2套參數(shù)的計(jì)算值較為接近實(shí)驗(yàn)值，更能反映真實(shí)的情況.此外，對(duì)比不同參數(shù)得到的結(jié)果表明新參數(shù)的正確性.

2.2 構(gòu)型變化

徑向分布函數(shù)(RDF)描述了體系中區(qū)域密度與平均密度的比值，對(duì)于固體材料而言，通過徑向分布函數(shù)可以得到不同原子之間的最近臨距離，可以很好地反映不同摻雜濃度時(shí)晶格結(jié)構(gòu)的變化.圖2分別給出了不同勢(shì)函數(shù)得到的不同摻雜濃度下，不同原子之間的第一臨位距離的變化情況.

圖2 不同勢(shì)函數(shù)計(jì)算的各離子間最近鄰距離Fig.2 Calculated nearest neighbor between different ions from different potentials

由圖2可知，雖然勢(shì)參數(shù)不同，但計(jì)算結(jié)果的變化趨勢(shì)基本一致，驗(yàn)證了模擬結(jié)果的合理性.此外，圖2中不同勢(shì)參數(shù)得到的計(jì)算結(jié)果對(duì)比可以發(fā)現(xiàn)，隨摻雜濃度的不斷增加，陽離子與氧離子之間的最近鄰距離不斷變小，且Gd3+-O2-間距始終大于Ce4+-O2-間距，而O2--O2-最近鄰距離卻不斷增大.氧空位產(chǎn)生后，O2-受力情況發(fā)生改變，處于氧空位臨位上的O2-會(huì)產(chǎn)生趨向氧空位的運(yùn)動(dòng)，從而使整個(gè)體系中O2--O2-的最近鄰距離增大，而整個(gè)過程中陽離子沒有發(fā)生大的跳躍，始終在初始位置附近振蕩，即陽離子之間的距離并沒有發(fā)生大的變化，因此陽離子與氧離子之間的最近鄰距離整體呈現(xiàn)減小的趨勢(shì).另外，由于在GDC中庫倫力是主要的作用形式，Ce4+-O2-間的吸引作用大于Gd3+-O2-間吸引作用，因此Gd3+-O2-間距始終大于Ce4+-O2-間距.

2.3 擴(kuò)散情況

一定溫度下，O2-能在GDC膜中擴(kuò)散是GDC膜導(dǎo)電的根本原因.為了研究氧離子的擴(kuò)散行為，模擬了溫度為1 273 K時(shí)不同摻雜濃度體系中氧離子的擴(kuò)散情況，系綜為NVT，時(shí)間步為1 fs，運(yùn)行100萬步.

均方位移(MSD)隨時(shí)間的變化情況體現(xiàn)了模擬過程中各類原子的擴(kuò)散行為，根據(jù)愛因斯坦的擴(kuò)散定律[11]，有

式中:D為粒子的擴(kuò)散系數(shù).當(dāng)時(shí)間較長(zhǎng)時(shí)，均方位移對(duì)時(shí)間變化曲線的斜率為6D.根據(jù)式(2)可以計(jì)算出不同勢(shì)參數(shù)下氧離子的擴(kuò)散系數(shù)，與實(shí)驗(yàn)值的對(duì)比結(jié)果如圖 3，圖 3中為了方便和實(shí)驗(yàn)值[9]對(duì)比，把CexGd1-xO2-0.5x中的x值換算成對(duì)應(yīng)的Gd2O3濃度.

圖3 1 273 K下不同勢(shì)參數(shù)計(jì)算的O2-的擴(kuò)散系數(shù)與實(shí)驗(yàn)值對(duì)比Fig.3 Calculated diffusion coefficients of oxygen from different potentials and experimental result

從圖3可以看出，隨著摻雜濃度的增大，O2-的擴(kuò)散系數(shù)開始不斷增大，但當(dāng)摻雜濃度達(dá)到一定程度后，擴(kuò)散系數(shù)卻開始逐漸減小，表明了O2-的擴(kuò)散行為由開始的激烈逐漸減緩，與實(shí)驗(yàn)現(xiàn)象一致.

事實(shí)上，O2-能夠在GDC膜中擴(kuò)散是由于氧空位的存在.初始階段，隨著摻雜濃度的增加，氧空位數(shù)量不斷增多，O2-的擴(kuò)散行為變得容易，即擴(kuò)散系數(shù)不斷增大;但同時(shí)要注意的是，氧空位顯示為正電荷，摻雜離子Gd＇Ce顯示為負(fù)電荷，兩者之間存在的庫倫吸引作用會(huì)阻礙O2-的擴(kuò)散，隨著氧空位不斷的增加，庫倫吸引不斷增強(qiáng)，對(duì)O2-擴(kuò)散的阻礙作用也不斷增強(qiáng);以上2種效果處于競(jìng)爭(zhēng)的狀態(tài)，當(dāng)摻雜濃度達(dá)到一定程度，后者的作用超過了前者，則O2-的擴(kuò)散行為開始變得困難，從而擴(kuò)散系數(shù)開始不斷減小，正如圖3所示.

此外，圖3中不同參數(shù)的計(jì)算結(jié)果雖然和實(shí)驗(yàn)值變化趨勢(shì)一致，但在具體數(shù)值上并不完全相同.實(shí)驗(yàn)中，當(dāng)Gd2O3摩爾濃度為10%時(shí)，擴(kuò)散系數(shù)最大;對(duì)于Lewis、Licia、Gotte這3組參數(shù)，當(dāng) Gd2O3摩爾濃度為8%～9%時(shí)得到的擴(kuò)散系數(shù)最大，且Lewis參數(shù)提出的時(shí)間最早，其得到的模擬結(jié)果絕對(duì)值與實(shí)驗(yàn)值差距最大，表明該參數(shù)適用范圍有限;而Hideaki的勢(shì)參數(shù)，最大擴(kuò)散系數(shù)對(duì)應(yīng)的Gd2O3為6%～8%;新提出的勢(shì)函數(shù)計(jì)算的結(jié)果與實(shí)驗(yàn)值變化情況較為接近，也驗(yàn)證了勢(shì)函數(shù)的合理性.同時(shí)，不同的計(jì)算結(jié)果說明了相對(duì)于晶體構(gòu)型而言，O2-的擴(kuò)散行為對(duì)勢(shì)參數(shù)的變化更為敏感.

3 結(jié)論

1)4套勢(shì)參數(shù)及本文提出的勢(shì)函數(shù)計(jì)算的GDC構(gòu)型變化趨勢(shì)基本一致.但提出時(shí)間較早的Lewis參數(shù)計(jì)算結(jié)果相對(duì)實(shí)驗(yàn)值差距較大，更適合于定性的分析，新勢(shì)函數(shù)較好的符合了實(shí)驗(yàn)情況，驗(yàn)證了該勢(shì)函數(shù)的正確性.

2)隨著摻雜濃度的增加，GDC的晶格常數(shù)不斷增加，O2--O2-的最近鄰距離不斷增大，而陽離子與氧離子之間的最近鄰距離不斷減小，且Gd3+-O2-間距始終大于Ce4+-O2-間距.

3)工作過程中，陽離子始終在初始位置發(fā)生震蕩，沒有發(fā)生較大的位移;O2-在GDC膜中發(fā)生了擴(kuò)散.隨著摻雜濃度的增加，O2-的擴(kuò)散開始變得激烈，擴(kuò)散系數(shù)增大，當(dāng)Gd2O3摩爾濃度達(dá)到6%～8%時(shí)，摻雜離子與氧空位的相互作用會(huì)阻礙O2-擴(kuò)散，使得其擴(kuò)散系數(shù)逐漸降低;且不同勢(shì)參數(shù)計(jì)算的擴(kuò)散系數(shù)有較大差異，表明了O2-的擴(kuò)散對(duì)勢(shì)參數(shù)的改變更為敏感.

[1]MORALES M，ROA J J，CAPDEVILA X G，et al.Mechanical properties at the nanometer scale of GDC and YSZ used as electrolytes for solid oxide fuel cells[J].Acta Materialia，2010，58(7):2504-2509.

[2]WANG Y G，LINAN A，SARAF L V，et al.Microstructure and ionic conductivity of alternating multilayer structured Gd-doped ceria and zirconia thin films[J].Journal of Materials Science，2009，44(8):2021-2026.

[3]YE F，MORI T，OU D R，et al.Dopant type depency of domain development in rare-earth-doped ceria:an explanation by computer simulation of defect clusters[J].Solid State Ionics，2009，180(20):1127-1132.

[4]WEI X，PAN W，CHENG L F，et al.Atomistic calculation of association energy in doped ceria[J].Solid State Ionics，2009，180(1):13-17.

[5]HULL S，NORBERG S T，AHMED I，et al.Oxygen vacancy ordering within anion-deficient ceria[J].Journal of Solid State Chemistry，2009，182(10):2815-2821.

[6]SAYLE T X T，SAYLE D C.Elastic deformation in ceria nanorods via a fluorite-to-rutile phase transition[J].ASC Nano，2010，4(2):879-886.

[7]GOTTE A，SPANGBERG D，BAUDIN M.Molecular dynamics study of oxygen self-diffusion in reduced CeO2[J].Solid State Ionics，2007，178(25):1421-1427.

[8]MINERVINI L，ZACATE M O，GRIMES R W.Defect cluster formation in M2O3-doped CeO2[J].Solid State Ionics，1999，116(4):339-349.

[9]CUI Z W，SUN Y，CHEN Y J.Semi-ab initio interionic potential for gadolinia-doped ceria[J].Solid State Ionics，2011，187(1):8-18.

[10]INABA H，SAGAWA R，HAYASHI H.Molecular dynamics simulation of gadolinia-doped ceria[J].Solid State I-onics，1999，122(2):95-103.

[11]陳正隆，徐為人，湯立達(dá).分子模擬的理論與實(shí)踐[M].北京:化學(xué)工業(yè)出版社，2007:110-115.CHEN Zhenglong，XU Weiren，TANG Lida.Theory and practice of molecular dynamics simulation[M].Beijing:Chemical Industry Press，2007:110-115.