苗　竹　張　?！詈Ｈ?/p>

(清華大學(xué)熱能工程系，熱科學(xué)與動力工程教育部重點實驗室，北京100084)

?

升溫?zé)Y(jié)過程中TiO2納米顆粒的原子分類分析(I)：表面原子識別

苗竹張海*楊海瑞

(清華大學(xué)熱能工程系，熱科學(xué)與動力工程教育部重點實驗室，北京100084)

開發(fā)表面原子識別模型，對單個TiO2顆粒升溫?zé)Y(jié)過程中表面原子進行分類研究。模型對顆?？臻g立方體網(wǎng)格化，利用標(biāo)準(zhǔn)球形顆粒體積積分法確定最佳網(wǎng)格尺寸為0.3 nm。通過表面網(wǎng)格識別實現(xiàn)表面原子分類，采用近鄰網(wǎng)格中外部網(wǎng)格數(shù)量(Next)作為準(zhǔn)則數(shù)判斷目標(biāo)網(wǎng)格是否為表面網(wǎng)格，確定最佳Next=9?；贚AMMPS軟件模擬了半徑為0.75 nm顆粒的升溫過程，發(fā)現(xiàn)系統(tǒng)能量弛豫速度明顯高于結(jié)構(gòu)弛豫速度；利用表面識別模型分類分析原子特性，表面原子平均位移大于內(nèi)部原子，且表面O原子遷移活性高于Ti原子；表面原子配位數(shù)低于內(nèi)部原子，佐證表面結(jié)構(gòu)規(guī)律性較差。研究結(jié)果為深入分析納米材料活性位等結(jié)構(gòu)分布奠定基礎(chǔ)。

TiO2；納米顆粒；分子動力學(xué)；表面識別；表面原子

www.whxb.pku.edu.cn

1　引言

二氧化鈦(TiO2)是一種典型n型半導(dǎo)體金屬氧化物1,2，其銳鈦礦型晶相具有較高活性3,4，在氣體傳感器5－7、催化劑8,9和染料敏化太陽能電池10,11等方面得到廣泛應(yīng)用。近年來，人們發(fā)現(xiàn)火焰合成法制備TiO2納米顆粒不僅經(jīng)濟，而且具有粒度分布可控、產(chǎn)物純度高的特點12,13；同時，火焰合成法可以直接成膜制備氣體傳感器，具有明顯優(yōu)勢14,15。

TiO2納米顆粒的表面特性是影響傳感器氣敏感應(yīng)性能的主要因素。氣體分子與TiO2顆粒表面發(fā)生相互作用，顆粒表面相比內(nèi)部區(qū)域存在更多結(jié)構(gòu)缺陷，充當(dāng)活性位用于氣敏感應(yīng)，顆粒表面特性分析對氣體傳感器性能研究有重要作用。分子動力學(xué)是研究納米顆粒結(jié)構(gòu)的有效方法16,17，可以直觀確定顆粒結(jié)構(gòu)在不同情況下的變化。Lu等18基于分子動力學(xué)方法研究了TiO2納米顆粒之間的相互作用，直觀觀察納米顆粒團聚融合過程中的結(jié)構(gòu)變化；Mushnoori19和Tillotson20等分別通過分子動力學(xué)模擬TiO2表面與不同反應(yīng)物作用過程，確定了反應(yīng)物與表面活性位作用機制，Zhang等21利用分子動力學(xué)方法研究了納米顆粒之間偶極相互作用對顆粒碰撞-燒結(jié)過程的影響。然而，現(xiàn)有研究主要基于納米顆粒結(jié)構(gòu)變化的直觀觀察，缺少對納米顆粒表面結(jié)構(gòu)的針對性分析?；诒砻嬖臃诸?，可以進一步確定活性位等特征結(jié)構(gòu)分布，對表征和分析納米顆粒結(jié)構(gòu)特性有重要作用。

顆粒表面結(jié)構(gòu)分析的基礎(chǔ)是表面原子識別，為此本文基于分子動力學(xué)原理建立了表面原子識別模型，對TiO2納米顆粒表面原子和內(nèi)部原子進行分類分析?？臻g網(wǎng)格化后以表面網(wǎng)格識別實現(xiàn)原子分類，通過近鄰網(wǎng)格種類分析判斷目標(biāo)網(wǎng)格是否為表面網(wǎng)格。利用開源軟件LAMMPS22,23模擬半徑0.75 nm顆粒升溫過程，通過表面原子識別模型對比分析升溫過程中表面原子和內(nèi)部原子遷移規(guī)律，確定顆粒表面特性。

2　模擬方法

本文計算模型為銳鈦礦型TiO2，晶胞參數(shù)a= 0.37845 nm、c=0.95143 nm、Z=4，屬于四方晶系24。構(gòu)建半徑0.75 nm的球形TiO2顆粒模型如圖1所示，隨機刪除若干外側(cè)原子保證顆粒整體電中性。圖中粉色原子為鈦(Ti)原子，綠色原子為氧(O)原子，顆粒中共有原子150個，其中Ti原子50個，O原子100個。

模擬通過開源軟件LAMMPS完成，將TiO2納米顆粒模型導(dǎo)入后模擬升溫過程。模擬在正則系綜(NVT)中進行，原子間相互作用勢選用Matsui-Akaogi(MA)勢24，該勢函數(shù)被認(rèn)為是刻畫金屬氧化物晶格結(jié)構(gòu)的良好選擇25。勢函數(shù)形式為：

其中Uij為i粒子與j粒子之間的相互作用勢，rij是i粒子與j粒子之間的距離，qi和qj分別是i粒子和j粒子的等效電荷，Ci和Cj分別是i粒子和j粒子之間的范德華力系數(shù)，f是一個標(biāo)準(zhǔn)力參數(shù)，其值大小為41.84 kJ·nm－1·mol－1，Ai和Aj分別為i粒子和j粒子的斥力半徑，Bi和Bj是與泡利斥力有關(guān)的參數(shù)。勢函數(shù)具體參數(shù)如表1所示。

系統(tǒng)邊界條件選擇非周期性(m,m,m)，計算時間步長為1 fs，每隔幾納秒輸出中間結(jié)果，對系統(tǒng)勢能、均方位移(MSD,mean-squared displacement)和配位數(shù)變化過程進行分析，其中MSD指顆粒中所有原子距離初始位置位移平方的平均值，用于描述系統(tǒng)結(jié)構(gòu)弛豫過程；配位數(shù)指該原子周圍晶格排布中配位原子個數(shù)，用于表征原子周圍晶格排布周期性。

單顆粒升溫模擬中，顆粒先在10 K下穩(wěn)定1 ns，然后經(jīng)過1 ns升溫至設(shè)定溫度，隨后穩(wěn)定弛豫約40 ns。在373－2373 K范圍內(nèi)選擇19個計算點進行穩(wěn)態(tài)計算，每隔幾個納秒輸出計算結(jié)果用于過程分析。模擬完成后，分析系統(tǒng)勢能和MSD變化規(guī)律，同時通過表面原子識別模型對表面原子進行識別和特性分析。

3　表面原子識別模型

表面原子識別模型旨在將表面原子和內(nèi)部原子進行分類分析，以便對比分析兩種原子的特性。由于顆粒弛豫之后原子坐標(biāo)分布無規(guī)則，且不連續(xù)，直接處理原子坐標(biāo)難度相對較大；為此，模型通過空間網(wǎng)格化后以網(wǎng)格為處理對象，根據(jù)網(wǎng)格內(nèi)是否有原子將網(wǎng)格分類為內(nèi)部網(wǎng)格和外部網(wǎng)格，再通過近鄰網(wǎng)格種類確認(rèn)目標(biāo)網(wǎng)格的所處位置，實現(xiàn)表面網(wǎng)格識別，表面網(wǎng)格中所包含的原子即為表面原子。

3.1空間網(wǎng)格化

采用立方體網(wǎng)格對納米顆粒所在空間進行網(wǎng)格化，理論上網(wǎng)格尺寸越小，原子類型識別的準(zhǔn)確性越高，然而由于顆粒中原子分布不連續(xù)，網(wǎng)格尺寸過小時顆粒內(nèi)部會出現(xiàn)無效空網(wǎng)格，干擾表面網(wǎng)格識別，如圖2所示，因此網(wǎng)格尺寸存在最佳優(yōu)化值。

圖1　TiO2納米顆粒模型Fig.1　TiO2nanoparticle model

表1　TiO2模擬的MA勢函數(shù)參數(shù)表Table 1　Parameters of MApotential function for TiO2simulation

網(wǎng)格尺寸優(yōu)化使用標(biāo)準(zhǔn)球形顆粒體積積分法。以半徑2 nm顆粒為對象，在0.2－1.0 nm范圍內(nèi)選擇九個值進行網(wǎng)格化，積分內(nèi)部網(wǎng)格體積得到顆粒計算體積。如圖3所示為顆粒計算體積隨網(wǎng)格尺寸的變化關(guān)系。隨網(wǎng)格尺寸增大，顆粒計算體積逐漸增大；相比顆粒體積理論值33.5 nm3，網(wǎng)格尺寸在0.3 nm時顆粒計算體積最為合理，即為最佳值。當(dāng)網(wǎng)格尺寸小于0.3 nm時，內(nèi)部無效空網(wǎng)格導(dǎo)致顆粒計算體積小于理論值；隨網(wǎng)格增大，邊界網(wǎng)格增大會導(dǎo)致更多外部體積計入顆粒內(nèi)，導(dǎo)致顆粒計算體積隨網(wǎng)格尺寸增大而增大。

綜上，通過立方體網(wǎng)格實現(xiàn)納米顆?？臻g網(wǎng)格化，并利用標(biāo)準(zhǔn)球形顆粒體積積分法確定網(wǎng)格最佳尺寸為0.3 nm，為表面原子識別模型確定最佳網(wǎng)格化方案。

圖2　不同尺寸網(wǎng)格效果示意圖Fig.2　Schematic of the meshing results with different mesh sizes

圖3　半徑為2 nm的顆粒計算體積與網(wǎng)格尺寸關(guān)系Fig.3　Relation between calculated particle volume and mesh size with particle radius of 2 nm

3.2表面原子識別

空間網(wǎng)格化后，通過表面網(wǎng)格識別實現(xiàn)表面原子和內(nèi)部原子分類，原理如圖4所示，選定某目標(biāo)網(wǎng)格，統(tǒng)計其近鄰網(wǎng)格中外部網(wǎng)格數(shù)量(Next)，Next越大，該目標(biāo)網(wǎng)格位于顆粒表面的可能性越高。選擇不同若干Next值分別進行表面原子識別，對比識別效果，確定最佳Next值。

圖5為四種Next選值情況下表面原子和內(nèi)部原子識別效果圖，統(tǒng)計得到表面原子與內(nèi)部原子數(shù)之比與Next關(guān)系如圖6所示。隨Next增大，表面原子與內(nèi)部原子數(shù)比逐漸減小，變化速度先減小后增大；在Next=9時，曲線變化率最小，表明此時表面原子識別效果最穩(wěn)定，即最佳Next。對比半徑為1、2和3 nm顆粒識別結(jié)果基本一致，進一步確定表面原子識別最佳的Next=9。最佳Next也對應(yīng)與目標(biāo)網(wǎng)格一個側(cè)面的近鄰網(wǎng)格數(shù)吻合，具有一定的理論依據(jù)。

綜上，通過表面網(wǎng)格識別實現(xiàn)表面原子和內(nèi)部原子分類，開發(fā)表面原子識別模型。通過近鄰網(wǎng)格種類分析法，以近鄰網(wǎng)格中外部網(wǎng)格數(shù)量Next為準(zhǔn)則數(shù)判斷目標(biāo)網(wǎng)格是否位于顆粒表面，并確定最佳Next=9。

圖4　表面網(wǎng)格識別方法示意圖Fig.4　Schematic of surface mesh identification method

圖5　不同Next情況下顆粒表面原子識別效果圖Fig.5　Surface atoms identification result with different Nextvalues

4　結(jié)果分析

4.1顆粒升溫過程

圖6　不同半徑(r)顆粒識別出表面原子與內(nèi)部原子數(shù)之比與Next關(guān)系Fig.6　Relation between ratio of surface atom number to interior atom number and Nextfor particles with different radii(r)

通過LAMMPS模擬顆粒升溫過程，各升溫過程中溫度變化曲線如圖7所示。圖8為半徑0.75 nm顆粒升溫過程系統(tǒng)勢能隨時間變化曲線，由圖可知，系統(tǒng)勢能隨溫度升高迅速增大，溫度穩(wěn)定后系統(tǒng)勢能也隨之穩(wěn)定；溫度越高，穩(wěn)態(tài)勢能越大。顆粒MSD隨時間變化曲線如圖9所示，對比勢能變化曲線可看出，MSD變化速度明顯慢于勢能變化速度；溫度穩(wěn)定后，系統(tǒng)MSD仍緩慢增大而后趨于穩(wěn)定。結(jié)果表明：系統(tǒng)勢能弛豫速度明顯快于系統(tǒng)結(jié)構(gòu)弛豫速度，不能僅靠能量穩(wěn)定判斷系統(tǒng)是否達(dá)到穩(wěn)態(tài)。當(dāng)溫度高于1173 K時，不同溫度穩(wěn)定MSD趨于一致，表明顆粒在高于此溫度時會發(fā)生熔化現(xiàn)象，液相中原子運動性一致，對應(yīng)一致的MSD結(jié)果。

4.2表面識別及特性分析

通過表面原子識別模型將顆粒中表面原子與內(nèi)部原子分類，其中表面原子將Ti原子與O原子分類處理。統(tǒng)計得到內(nèi)部原子、表面原子、表面Ti原子和表面O原子分別在573、973、1173、1573和2073 K情況下的原子平均位移，如圖10所示。由圖可知，所有原子平均位移均隨溫度升高而增大；在同一溫度下，表面原子平均位移大于內(nèi)部原子，且表面中O原子平均位移大于Ti原子，當(dāng)溫度較高顆粒熔化時Ti原子與O原子的位移相差不大；表明TiO2納米顆粒表面原子運動強度高于內(nèi)部原子，且表面O原子遷移活性強于Ti原子。對比表面原子和內(nèi)部原子配位數(shù)如圖11所示，隨溫度升高兩類原子配位數(shù)均降低；各溫度情況下，表面原子配位數(shù)低于內(nèi)部原子配位數(shù)。

圖7　不同升溫過程溫度變化曲線Fig.7　Temperature curves of different heating processes

圖8　不同升溫過程中勢能隨時間變化曲線Fig.8　Variation of potential energy with time in different heating processes

圖9　不同升溫過程中顆粒均方位移(MSD)隨時間變化曲線Fig.9　Variation of particle mean-squared displacement(MSD)with time in different heating processes

圖10　不同溫度下各類原子平均位移對比Fig.10　Averaged displacement comparison of different kinds of atoms under different temperatures

圖11　不同溫度下表面原子與內(nèi)部原子配位數(shù)對比Fig.11　Coordination number comparison of interior atoms and surface atoms under different temperatures

5　結(jié)論

基于分子動力學(xué)開發(fā)納米顆粒表面原子識別模型，對顆粒中原子進行分類分析。采用立方體網(wǎng)格對顆?？臻g網(wǎng)格化，利用標(biāo)準(zhǔn)球形顆粒體積積分法確定最佳網(wǎng)格尺寸為0.3 nm；通過近鄰網(wǎng)格種類分析識別表面網(wǎng)格，以近鄰?fù)獠烤W(wǎng)格數(shù)Next為準(zhǔn)則數(shù)判斷表面網(wǎng)格，并確定最佳Next=9。模擬半徑0.75 nm納米顆粒升溫過程，分析系統(tǒng)勢能和MSD變化規(guī)律，發(fā)現(xiàn)系統(tǒng)能量弛豫速度明顯高于結(jié)構(gòu)弛豫速度，不能通過能量穩(wěn)定判斷系統(tǒng)穩(wěn)態(tài)。通過表面原子識別模型對比分析表面原子與內(nèi)部原子特性，其中表面原子平均位移大于內(nèi)部原子，且表面O原子遷移活性高于表面Ti原子；內(nèi)部原子配位數(shù)整體高于表面原子。綜上，表面原子識別模型效果穩(wěn)定，能夠準(zhǔn)確分辨不同溫度下表面原子和內(nèi)部原子性質(zhì)差異，可以用于進一步分析表面特征結(jié)構(gòu)；該模型可以推廣應(yīng)用與其他晶體顆粒結(jié)構(gòu)，但最佳網(wǎng)格尺寸與該晶體的晶格參數(shù)有關(guān)。

References

(1)Liu,M.;Li,H.M.;Wang,W.J.Catal.Today 2016,264,236. doi:10.1016/j.cattod.2015.08.044

(2)Bet-moushoul,E.;Mansourpanah,Y.;Farhadi,K.;Tabatabaei, M.Chem.Eng.J.2016,283,29.doi:10.1016/j.cej.2015.06.124

(3)Zhang,J.;Wu,B.;Huang,L.H.;Liu,P.L.;Wang,X.Y.;Lu,Z. D.;Wu,G.L.;Zhang,E.P.;Wang,H.B.;Kong,Z.;Xi,J.H.;Ji, Z.G.J.Alloy.Compd.2016,661,441.doi:10.1016/j. jallcom.2015.11.225

(4)Xu,L.;Tang,C.Q.;Huang,Z.B.Acta Phys.-Chim.Sin.2010, 26(5),1401.[徐凌,唐超群,黃宗斌.物理化學(xué)學(xué)報,2010, 26(5),1401.]doi:10.3866/PKU.WHXB20100526

(5)Zhao,M.;Wang,J.X.;Feng,C.H.;Zou,B.;Chen,C.;Wang,Z. Y.;Wu,F.Q.;Zou,L.H.Acta Phys.-Chim.Sin.2007,23(7), 1003.[趙萌,王金興,馮彩慧,鄒博,陳騁,王竹儀,吳鳳清,鄒樂輝.物理化學(xué)學(xué)報,2007,23(7),1003.]doi:10.3866/ PKU.WHXB20070708

(6)Moon,H.G.;Shim,Y.S.;Jang,H.W.;Kim,J.S.;Choi,K.J.; Kang,C.Y.;Choi,J.W.;Park,H.H.;Yoon,S.J.Sens.Actuat. B-Chem.2010,149,116.doi:10.1016/j.snb.2010.06.014

(7)Chaudhari,G.N.;Bende,A.M.;Bodade,A.B.;Patil,S.S.; Sapkal,V.S.Sens.Actuat.B-Chem.2006,115,297. doi:10.1016/j.snb.2005.09.014

(8)Ao,P.;Xu,X.S.;Xu,X.X.;Li,J.H.;Yan,X.H.Acta Phys.-Chim.Sin.2014,30(5),950.[敖平,許響生,徐瀟瀟,李加衡,嚴(yán)新煥.物理化學(xué)學(xué)報,2014,30(5),950.] doi:10.3866/PKU.WHXB201403111

(9)Marin,R.P.;Ishikawa,S.;Bahruji,H.;Shaw,G.;Kondrat,S. A.;Miedziak,P.J.;Morgan,D.J.;Taylor,S.H.;Bartley,J.K.; Edwards,J.K.;Bowker,M.;Ueda,W.;Hutchings,G.J.Appl. Catal.A 2015,504,62.doi:10.1016/j.apcata.2015.02.023

(10)Gao,S.W.;Lan,Z.;Wu,W.X.;Que,L.F.;Wu,J.H.;Lin,J. M.;Huang,M.L.Acta Phys.-Chim.Sin.2014,30(3),446.[高素雯,蘭章,吳晚霞,闕蘭芳,吳季懷,林建明,黃妙良.物理化學(xué)學(xué)報,2014,30(3),446.]doi:10.3866/PKU. WHXB201401022

(11)Guo,W.;Wang,K.;Shen,Y.H.;Zhang,H.;Weng,T.;Ma,T.L. Acta Phys.-Chim.Sin.2013,29(1),82.[郭薇,王開,沈藝華,張賀,翁韜,馬廷麗.物理化學(xué)學(xué)報,2013,29(1),82.] doi:10.3866/PKU.WHXB201211071

(12)Zhang,Y.Y.;Li,S.Q.;Deng,S.L.;Yao,Q.;Tse,S.D.J.Aeros. Sci.2012,44,71.doi:10.1063/1.3555633

(13)Zhao,B.;Uchikawa,K.;Wang,H.Proc.Combust.Inst.2007, 31,851.doi:10.1016/j.proci.2006.08.064

(14)Tolmachoff,E.;Memarzadeh,S.;Wang,H.J.Phys.Chem.C 2011,115,21620.doi:10.1021/jp206061h

(15)Zhao,B.;Uchikawa,K.;McCormick,J.R.;Ni,C.Y.;Chen,J. G.;Wang,H.Proc.Combust.Inst.2005,30,2569.doi:10.1016/ j.proci.2004.08.146

(16)Fereidoon,A.;Aleaghee,S.;Taraghi,I.Comp.Mater.Sci.2015, 102,220.doi:10.1016/j.commatsci.2015.02.044

(17)Gautam,S.;Cole,D.Chem.Phys.2015,458,68.doi:10.1016/j. chemphys.2015.07.012

(18)Lu,J.;Liu,D.M.;Yang,X.N.;Zhao,Y.;Liu,H.X.;Tang,H.; Cui,F.Y.Appl.Surf.Sci.2015,357,1114.doi:10.1016/j. apsusc.2015.09.142)

(19)Mushnoori,S.;Chong,L.;Dutt,M.Materials Today: Proceedings 2016,3,513.doi:1016/j.matpr.2016.01.082

(20)Tillotson,M.J.;Brett,P.M.;Bennett,R.A.;Ricardo G.C.Surf. Sci.2015,632,142.doi:10.1016/j.susc.2014.09.017

(21)Zhang,Y.Y.;Li,S.Q.;Yan,W.;Yao,Q.;Tse,S.J.Chem.Phys. 2011,124,084501.doi:10.1063/1.3555633

(22)Freitas,R.;Asta,M.;Koning,M.D.Comp.Mater.Sci.2016, 112,333.doi:10.1016/j.commatsci.2015.10.050

(23)Martin,M.G.;Thompson,A.P.Fluid Phase Equilib.2004,217, 105.doi:10.1016/j.fluid.2003.06.007

(24)Zhang,Y.Y.Dynamics of Nanoparticles in Stagnation Flames. Ph.D.Dissertation.Tsinghua University,Beijing,2013.[張易陽.基于滯止火焰合成的高溫納米顆粒動力學(xué)研究[D].北京:清華大學(xué),2013.]

(25)Zhang,H.Z.;Chen,B.;Banfield,J.F.Phys.Chem.Chem.Phys. 2009,11(14),2553.doi:10.1039/b819623k

Atom Identification and Analysis of TiO2Nanoparticles in the Heating and Sintering Process(I):Surface Atom Identification

MIAO ZhuZHANG Hai*YANG Hai-Rui
(Key Laboratory for Thermal Science and Power Engineering of Ministry of Education,Department of Thermal Engineering, Tsinghua University,Beijing 100084,P.R.China)

In this paper,we develop a surface atom identification model to perform atom identification and analysis of the surface atoms of a single TiO2nanoparticle during the heating and sintering process.Cubic mesh was used to obtain the particle structure mesh,and the optimal mesh size of 0.3 nm was determined by volumetric integration of a spherical particle.Surface atoms were classified according to identification of surface meshes,and the number of external meshes in all of the neighbor meshes(Next)was set as a criterion to determine whether the target mesh was a surface mesh.The optimal value was Next=9.LAMMPS was used to simulate the heating process of a particle with a radius of 0.75 nm.The results show that energy relaxation is significantly faster than structure relaxation.The atom classification analysis with the developed surface atom identification model shows that displacement of the surface atoms is larger than the interior atoms,and surface O atoms are more active in migration than surface Ti atoms.The coordination number of surface atoms is lower than that of interior atoms.The present study provides fundamental information for analyzing the active structure distribution of nanoparticles.

TiO2;Nanoparticle;Molecular dynamics;Surface identification;Surface atom

March 14,2016;Revised:April 26,2016;Published on Web:April 26,2016.

O647

10.3866/PKU.WHXB201604262

*Corresponding author.Email:haizhang@tsinghua.edu.cn;Tel:+86-10-62773153.

The project was supported by the National Natural Science Foundation of China(51176095,51476088).

國家自然科學(xué)基金(51176095,51476088)資助項目

?Editorial office ofActa Physico-Chimica Sinica

［Article］