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

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

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

1 GDC勢函數(shù)

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

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

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

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

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

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

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

2.1 晶格常數(shù)

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

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

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

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

2.2 構(gòu)型變化

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

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

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

2.3 擴散情況

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

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

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

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

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

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

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

3 結(jié)論

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

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

3)工作過程中，陽離子始終在初始位置發(fā)生震蕩，沒有發(fā)生較大的位移;O2-在GDC膜中發(fā)生了擴散.隨著摻雜濃度的增加，O2-的擴散開始變得激烈，擴散系數(shù)增大，當(dāng)Gd2O3摩爾濃度達(dá)到6%～8%時，摻雜離子與氧空位的相互作用會阻礙O2-擴散，使得其擴散系數(shù)逐漸降低;且不同勢參數(shù)計算的擴散系數(shù)有較大差異，表明了O2-的擴散對勢參數(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á).分子模擬的理論與實踐[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.