升溫燒結過程中TiO2納米顆粒的原子分類分析(I)：表面原子識別

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

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

?

升溫燒結過程中TiO2納米顆粒的原子分類分析(I)：表面原子識別

苗竹張海*楊海瑞

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

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

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

www.whxb.pku.edu.cn

1　引言

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

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

顆粒表面結構分析的基礎是表面原子識別，為此本文基于分子動力學原理建立了表面原子識別模型，對TiO2納米顆粒表面原子和內部原子進行分類分析。空間網格化后以表面網格識別實現原子分類，通過近鄰網格種類分析判斷目標網格是否為表面網格。利用開源軟件LAMMPS22,23模擬半徑0.75 nm顆粒升溫過程，通過表面原子識別模型對比分析升溫過程中表面原子和內部原子遷移規律，確定顆粒表面特性。

2　模擬方法

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

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

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

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

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

3　表面原子識別模型

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

3.1空間網格化

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

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

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

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

綜上，通過立方體網格實現納米顆粒空間網格化，并利用標準球形顆粒體積積分法確定網格最佳尺寸為0.3 nm，為表面原子識別模型確定最佳網格化方案。

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

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

3.2表面原子識別

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

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

綜上，通過表面網格識別實現表面原子和內部原子分類，開發表面原子識別模型。通過近鄰網格種類分析法，以近鄰網格中外部網格數量Next為準則數判斷目標網格是否位于顆粒表面，并確定最佳Next=9。

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

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

4　結果分析

4.1顆粒升溫過程

圖6　不同半徑(r)顆粒識別出表面原子與內部原子數之比與Next關系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顆粒升溫過程系統勢能隨時間變化曲線，由圖可知，系統勢能隨溫度升高迅速增大，溫度穩定后系統勢能也隨之穩定；溫度越高，穩態勢能越大。顆粒MSD隨時間變化曲線如圖9所示，對比勢能變化曲線可看出，MSD變化速度明顯慢于勢能變化速度；溫度穩定后，系統MSD仍緩慢增大而后趨于穩定。結果表明：系統勢能弛豫速度明顯快于系統結構弛豫速度，不能僅靠能量穩定判斷系統是否達到穩態。當溫度高于1173 K時，不同溫度穩定MSD趨于一致，表明顆粒在高于此溫度時會發生熔化現象，液相中原子運動性一致，對應一致的MSD結果。

4.2表面識別及特性分析

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

圖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　不同溫度下表面原子與內部原子配位數對比Fig.11　Coordination number comparison of interior atoms and surface atoms under different temperatures

5　結論

基于分子動力學開發納米顆粒表面原子識別模型，對顆粒中原子進行分類分析。采用立方體網格對顆粒空間網格化，利用標準球形顆粒體積積分法確定最佳網格尺寸為0.3 nm；通過近鄰網格種類分析識別表面網格，以近鄰外部網格數Next為準則數判斷表面網格，并確定最佳Next=9。模擬半徑0.75 nm納米顆粒升溫過程，分析系統勢能和MSD變化規律，發現系統能量弛豫速度明顯高于結構弛豫速度，不能通過能量穩定判斷系統穩態。通過表面原子識別模型對比分析表面原子與內部原子特性，其中表面原子平均位移大于內部原子，且表面O原子遷移活性高于表面Ti原子；內部原子配位數整體高于表面原子。綜上，表面原子識別模型效果穩定，能夠準確分辨不同溫度下表面原子和內部原子性質差異，可以用于進一步分析表面特征結構；該模型可以推廣應用與其他晶體顆粒結構，但最佳網格尺寸與該晶體的晶格參數有關。

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.[徐凌,唐超群,黃宗斌.物理化學學報,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.[趙萌,王金興,馮彩慧,鄒博,陳騁,王竹儀,吳鳳清,鄒樂輝.物理化學學報,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.[敖平,許響生,徐瀟瀟,李加衡,嚴新煥.物理化學學報,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.[高素雯,蘭章,吳晚霞,闕蘭芳,吳季懷,林建明,黃妙良.物理化學學報,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.[郭薇,王開,沈藝華,張賀,翁韜,馬廷麗.物理化學學報,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.[張易陽.基于滯止火焰合成的高溫納米顆粒動力學研究[D].北京:清華大學,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).

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

?Editorial office ofActa Physico-Chimica Sinica

［Article］