劉賽賽 王立峰

（南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016）

氧化鋅納米線振動問題研究

劉賽賽 王立峰?

（南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016）

采用連續(xù)介質理論與分子動力學模擬相結合的方法，研究了氧化鋅納米線的振動問題.建立了氧化鋅納米線核殼模型，解釋其等效楊氏模量及壓電常數(shù)的尺寸效應.通過連續(xù)介質理論求得氧化鋅納米線振動固有頻率，并與分子動力學模擬得到的結果進行對比.研究表明，氧化鋅納米線在極化方向的等效拉伸楊氏模量隨著橫截面尺寸的增加而逐漸增大，且通過核殼模型分別求得核、殼拉伸楊氏模量.擬合得到的等效拉伸楊氏模量與分子動力學方法獲得的等效拉伸楊氏模量符合得很好.根據(jù)連續(xù)介質理論得到等效彎曲楊氏模量，發(fā)現(xiàn)等效彎曲楊氏模量也隨著橫截面尺寸的增加而增大.氧化鋅納米線極化方向的壓電耦合能力比一般壓電陶瓷好，壓電常數(shù)隨著橫截面尺寸的增加逐漸減小.氧化鋅納米線在不同溫度條件下的振動頻率沒有明顯變化，在不同外電場條件下的振動頻率有顯著變化.分子動力學模擬得到不同橫截面尺寸的氧化鋅納米線振動頻率不同.根據(jù)連續(xù)介質理論，求得懸臂Timoshenko梁模型相應尺寸的振動頻率，發(fā)現(xiàn)橫截面的尺寸越大，連續(xù)介質理論與分子動力學模擬得到的振動頻率越接近.

氧化鋅納米線， 分子動力學， 尺寸效應， 壓電效應， 振動

引言

氧化鋅納米線在光學、力學等方面具有獨特的壓電性、光電性，而使其備受關注.研究發(fā)現(xiàn)氧化鋅是一種寬帶隙半導體，在室溫下仍有很強的導電性，因此科學界認為氧化鋅有望取代氮化鎵成為半導體發(fā)光二極管和半導體激光二極管的新材料［1］.

國內外對于納米氧化鋅的力學、電學性能等做了大量研究.王中林等利用氧化鋅納米線發(fā)明了納米發(fā)電機［2，3］.Corso 等采用量子力學從頭計算方法，發(fā)現(xiàn)在四面共價半導體中，氧化鋅納米線的壓電性能最好［4］.Hill等通過局域密度近似的平面波偽勢密度泛函理論，發(fā)現(xiàn)可以通過改變氧化鋅的晶格常數(shù)實現(xiàn)其壓電性質的變化［5］.基于Berry相理論［6-8］，Noel等得到氧化鋅的壓電常數(shù)和 Wannier函數(shù)理論結果符合很好［9］，Catti等測量了纖鋅礦、閃鋅礦及硫化鋅的壓電張量［10］.Lucas等通過納米壓痕技術，利用原子力顯微鏡測量了氧化鋅納米帶的楊氏模量［11］，Asthana等實驗研究了氧化鋅納米線機械性能的尺寸效應［12］，Xu等用掃描電子顯微鏡測量了氧化鋅納米線在不同載荷作用下的機械性能［13］，Agrawal等利用實驗和計算相結合的方法得到氧化鋅納米線的體楊氏模量約為140GPa［14］.然而由于生產(chǎn)技術、納米氧化鋅各項異性、實驗設備精確度、納米氧化鋅材料的選取等原因，對于納米氧化鋅楊氏模量的實驗測量結果很不一致［15］.Ewald［16-18］求和方法的提出，使通過計算機模擬來研究氧化鋅納米線、納米管等的機械性質得到迅速發(fā)展［19-21］.Zaoui等研究了壓強對氧化鋅剪切模量的影響［22］，Wang等分析了溫度對氧化鋅納米線機械行為的影響［23］，Sun等研究了閃鋅礦壓強、溫度與體積之間的關系［24］.

壓電梁振動問題一直是備受關注的問題，曹樹謙和郭抗抗研究了壓電發(fā)電懸臂梁非線性動力學響應問題，給出了不同結構參數(shù)及外激勵參數(shù)下系統(tǒng)響應特性的數(shù)值結果和實驗結果［25］.彭劍等研究了時滯反饋作用下壓電梁的參數(shù)共振問題，分析了軸力等因素對穩(wěn)定性區(qū)域的影響［26］.隨著微機電系統(tǒng)（MEMS）技術的發(fā)展［27，28］，為了在單位體積內實現(xiàn)更多的功能，人們開始研制納機電系統(tǒng)（NEMS）.納機電系統(tǒng)的研究有望促進醫(yī)療技術、信息傳播、軍事國防、生產(chǎn)制造等的發(fā)展.Liu和Wang通過半量子化分子動力學方法研究了單壁碳納米管振動問題，發(fā)現(xiàn)在低溫、高頻振動情況中，碳納米管具有明顯的量子效應［29］.Ke等基于非線性理論研究了壓電納米梁熱電力耦合振動（thermoelectricmechanical vibration）［30］.納米壓電梁結構是最簡單的納機電系統(tǒng)，壓電梁振動研究是基于振動的機械傳感器研制的基礎，對壓電懸臂梁振動的研究，有利于提高以新型壓電材料為基礎的能量采集器的效率.本文主要研究懸臂氧化鋅納米線振動問題，關注溫度、外電場強度對振動頻率的影響.

1 氧化鋅納米線力學性質的分子動力學模擬

1.1 模型建立

本文研究纖鋅礦結構氧化鋅納米線，初始結構如圖1（a），紅色為氧原子、灰色為鋅原子，圖1（b）為氧化鋅納米線模型，兩端黃色為固定層，橫截面如圖1（c），計算所用到的氧化鋅納米線尺寸、序號見表1．b、h、L分別為模型的寬、高和長，單位為納米，Δb、Δh 是殼不同邊的厚度［19，31，32］，取 Δb＝0．406 nm，Δh＝0．469nm，z軸垂直于延長方向．

圖1 氧化鋅結構圖Fig.1 Structures of ZnO

氧化鋅納米線原子間的相互作用采用Buckingham 對勢［33，34］:

式中rij是離子i和離子j之間距離，qi是離子i所帶電荷，A、ρ、C是該勢的參數(shù)，取值見表 2．等式右邊第一項是長程庫倫作用，第二項是短程相互排斥作用，第三項是短程相互吸引作用．氧離子之間既有庫倫作用又有短程作用，鋅離子之間只有庫倫作用，氧離子和鋅離子之間忽略短程相互吸引作用，只有庫倫作用和短程相互排斥作用，計算長程庫倫作用采用 Wolf算法［35，36］，對庫倫作用中 r-1項進行截斷，提高計算效率．本文阻尼系數(shù)取0．4，截斷半徑取1nm，可實現(xiàn)Wolf求和的最佳收斂［35］.在分子動力學中Wolf算法要求采用周期性邊界條件，我們將氧化鋅納米線模型放入較大的真空模擬盒中，避免不同模擬盒中的氧化鋅納米線之間存在相互作用，近似模擬單根納米線的情況．

表1 氧化鋅模型大小及序號Table 1 Size and serial numbers of ZnO nanowires

1.2 等效拉伸楊氏模量

壓電材料單位體積總勢能［37］:

式中Uv是單位體積總勢能，是外電場恒定時的彈性系數(shù)，εij是應變張量，是應變恒定時的介電常數(shù)，Ei是外電場矢量．當外電場為零時，z方向的楊氏模量簡化為:

式中U是總勢能，V是氧化鋅納米線的體積．

將模型完全弛豫后固定兩端，固定層向兩側各拉伸 ΔL＝0．002nm，應變 ε＝2ΔL／L，再次弛豫后得到總勢能U，如此重復50次，得到總勢能與應變關系曲線圖和應力與應變關系圖．圖2為表1中No.2納米線的總勢能、應力與應變關系圖．根據(jù)總勢能和應變關系，由式（3）可得到拉伸楊氏模量．

圖2 No.2氧化鋅納米線Fig.2 ZnO nanowire （No.2）

在軸向拉伸載荷作用下，氧化鋅納米線楊氏模量可看作是由核、殼的楊氏模量共同作用的結果如圖1（c）．設核的楊氏模量為Yc，殼的楊氏模量為Yc+Y1，則拉伸等效楊氏模量為［38］:

式中Ys為等效拉伸楊氏模量，A為模型橫截面面積，A1為殼的橫截面面積，A2為核的橫截面面積．

用截面內的原子個數(shù)比代替截面面積比，對表1中的氧化鋅納米線進行分子動力學模擬，得到的拉伸楊氏模量與截面內的原子個數(shù)比通過最小二乘法擬合，如圖3．得到核、殼楊氏模量分別為Yc＝184．02GPa，Yc+Y1＝81．10GPa，進一步得到連續(xù)介質等效楊氏模量Ys．從圖4可以看出，連續(xù)介質等效楊氏模量Ys與分子動力學模擬得到的楊氏模量Y符合得較好．

圖3 拉伸楊氏模量與原子數(shù)比擬合曲線Fig.3 Fitting curve of tension Young′s modulus with the ratio of atoms shell and bulk

圖4 拉伸楊氏模量與截面寬度的關系Fig.4 Relationship of Young′s modulus and width for ZnO nanowires

1.3 等效彎曲楊氏模量

為使模型簡化，假設梁在受載荷作用時，垂直于軸線的平面在彎曲過程中始終垂直于軸線，即梁受純彎曲載荷作用．

當梁受到純彎曲載荷作用時，截面上的彎矩為

式中ρ′為梁的曲率半徑，I1、I2分別是殼、核截面慣性矩．

由彎矩平衡方程得:

式中Yb為等效彎曲楊氏模量，I為模型整體的截面慣性矩．

根據(jù) Yc＝ 184．02GPa，Yc+Y1＝ 81．10GPa，得到等效彎曲楊氏模量與截面y方向尺寸關系如圖5所示，Yby、Ybx分別是氧化鋅納米線沿y、x方向的彎曲楊氏模量．

圖5 等效彎曲楊氏模量與截面高度的關系Fig.5 Relationship of equivalent bending Young′s modulus and height for ZnO nanowires

從圖5可以看出，氧化鋅納米線的等效彎曲楊氏模量都隨著截面尺寸的增加而增大．然而，由于氧化鋅的各項異性，如表1近似正方形橫截面的氧化鋅納米線，在上述兩個方向的等效彎曲楊氏模量不完全相同．

1.4 壓電常數(shù)的計算

極化強度由電子-原子核和鋅氧離子相對位移變化共同作用的結果［34］，當忽略前者時，氧化鋅極化方向的平均極化強度由下式?jīng)Q定:

本文僅研究氧化鋅納米線極化方向壓電性質，有:

由式（9）和式（10）得:

式中D是電位移，d是壓電常數(shù)，σ是應力．根據(jù)上式得到表1中壓電常數(shù)與其橫截面x方向尺寸關系如圖6所示．從圖6中可以看出隨著氧化鋅納米線尺寸的增加，極化方向的壓電常數(shù)逐漸減小，并趨于某一定值．對于序號10對應的橫截面最大的氧化鋅納米線，其壓電常數(shù)為 1．65×10-11C／N，與Momeni等的模擬結果符合較好［32］．

圖6 壓電常數(shù)與截面寬度的關系Fig.6 Relationship of piezoelectric constant and section width

2 固有振動

2.1 分子動力學模擬

本文分子動力學模擬過程中采用正則系綜（NVT），用 Nose-Hoover控溫方法［39，40］使系統(tǒng)的溫度保持在一定范圍內.

對序號2氧化鋅納米線分別進行300K、500K、1000K控溫模擬.從圖7幅頻曲線可以看出，溫度變化對氧化鋅納米線的振動頻率幾乎沒有影響.

對序號2模型，沿著z軸方向加外電場E，強度從-1.0V／nm 到 1.0V／nm，控制溫度為 300K.如圖8，橫坐標為電場強度，縱坐標為前四階振動頻率.可以看出，對于本文研究的氧化鋅納米線，當電場強度大于-0.20V／nm時，其振動頻率隨著電場強度的增大而增大；當電場強度小于-0.37V／nm，其振動頻率隨著電場強度的減小而增大；當電場強度在約-0.37V／nm到-0.20V／nm時，氧化鋅納米線運動幅度過大，無法確定振動頻率．

圖7 溫度對振動頻率的影響Fig.7 Effect of temperature on the frequencies

圖8 外電場強度對振動頻率的影響Fig.8 Effect of electric field intensity on the frequencies

2.2 懸臂Timoshenko梁模型

為了預測氧化鋅納米線振動行為，建立懸臂Timoshenko梁模型［41］，方向為沿軸向 x方向:

式中:E為楊氏模量，I為截面慣性矩，φ為彎曲轉角，y為撓度，k為截面因子，A為橫截面面積，G為剪切彈性模量，ρ為密度，t為時間變量，F(xiàn)＝k1nqE為軸向力，n為單層鋅（或氧）原子數(shù)，k1為比例系數(shù)，q為原子帶電量，E為外電場強度．

采用懸臂梁邊界條件［29］:

取截面因子 k＝5／6，密度 ρ＝5606kg／m3，泊松比υ＝0．3．將上文得到的等效彎曲楊氏模量Yby代入行列式，得到不同厚度的氧化鋅納米線沿y方向振動的前四階振動頻率和分子動力學結果（MD）對比如圖9．把等效軸力代入行列式，得到連續(xù)介質理論預測的振動頻率，圖10為氧化鋅納米線的振動頻率隨電場強度的變化．從圖9可以看出，氧化鋅納米線的橫截面尺寸越大，其連續(xù)介質力學結果和分子動力學模擬結果符合越好．從圖10可以看出，氧化鋅納米線連續(xù)介質力學模型得到的振動頻率和其分子動力學模擬得到的振動頻率符合得很好．

圖9 氧化鋅納米線振動頻率Fig.9 Frequencies of ZnO nanowires

圖10 軸力作用下氧化鋅納米線振動頻率Fig.10 Frequencies of ZnO nanowires under the axial force

3 結論

本文通過分子動力學模擬和連續(xù)介質理論相結合的方法，研究了氧化鋅納米線振動問題，得到如下結論:

（1）溫度變化對懸臂氧化鋅納米線的振動特性沒有明顯影響，氧化鋅納米線在極化方向的拉伸楊氏模量、垂直極化方向的彎曲楊氏模量隨著橫截面尺寸的增加逐漸增大.

（2）氧化鋅納米線極化方向的壓電常數(shù)隨著橫截面尺寸的增加逐漸減小，并趨向一定值；當極化方向電場強度大于-0.20V／nm（或小于-0.37V／nm）時，增加（或反向增加）極化方向外電場強度可以顯著提高氧化鋅納米線的振動頻率.

（3）懸臂Timoshenko梁理論求得的振動頻率和分子動力學模擬得到的振動頻率，在模型橫截面尺寸較小時差別明顯，在尺寸較大時符合得較好.

1 Ley L，Pollak R A，McFeely F R，et al.Total valenceband densities of states of III-V and II-VI compounds from x-ray photoemission spectroscopy.Physical Review B Condensed Matter， 1974，9（2）:600～621

2 Wang Z L，Song J H.Piezoelectricnanogenerators based on zinc oxide nanowire arrays.Science， 2006，312:242

3 Song J H，Zhou J， Wang Z L.Piezoelectric andsemiconducting coupled power generating process of a single ZnO belt／wire.A technology for harvesting electricity from the environment.Nano Letter， 2006， 6:1656

4 Corso A D，Posernak M，Resta R，et al.Ab initio study of piezoelectricity and spontaneous polarization in ZnO.Physical Review B， 1994，50（15）:10715～10721

5 Hill N A，Waghmare U.First-principles study of strain-electronic interplay in ZnO:Stress and temperature dependence of the piezoelectric constants.Physical Review B， 2000，62:8002

6 King-Smith R D，Vanderbilt D.Theory of polarization of crystalline solids.Physical Review B Condensed Matter，1993，47（3）:1651

7 Vanderbilt D，King-Smith R D.Electric polarization as a bulk quantity and its relation to surface charge.Physical Review B Condensed Matter， 1993，48（7）:4442

8 Resta R.Macroscopic polarization in crystalline dielectrics:the geometric phase approach.Reviews of Modern Physics， 1994，66:899

9 Novel Y，Zicovich-Wilson C M，Civalleri B，et al.Polarization properties of ZnO and BeO:An ab initio study through the berry phase and wannier functions approaches.Physical Review B， 2001，65（65）:014111

10 Catti M，Noel Y，Dovesi R.Full piezoelectric tensors of wurtzite and zinc blende ZnO and ZnS by first-principles calculations.Journal of Physics and Chemistry of Solids，2003，64（11）:2183～2190

11 Lucas M，Mai W J，Yang R，et al.Aspect ratio dependence of the elastic properties of ZnO nanobelts.Nano Letters， 2007，7（5）:1314

12 Asthana A，Momeni K，Prasad A，et al In situ observa-tion of size-scale effects on the mechanical properties of ZnO nanowires.Nanotechnology， 2011，22（26）:265712

13 Xu F，Qin Q Q，Mishra A，et al.Mechanicalproperties of ZnO nanowires under different loading modes.Nano Research， 2010，3（4）:271

14 Agrawal R，Peng B，Gdoutos E E，et al.Elasticity size effects in ZnO nanowires——A combined experimentalcomputational approach.Nano Letters， 2008，8（11）:3668

15 Dai L，Cheong C D，Sow C H，et al.Molecular dynamics simulation of ZnO nanowires:size effects， defects， and super ductility.Langmuir， 2010，26（2）:1165

16 Darden T，York D，Pedersen L.Particle mesh Ewald:An N?log（N） method for Ewald sums in large systems.Journal of Chemical Physics， 1993，98（12）:10089

17 Bródka A， Sliwiński P.Three-dimensional Ewald method with correction term for a system periodic in one direction.Journal of Chemical Physics， 2004，120（12）:5518～5523

18 Raymand D，van Duin A C T，Baudin M，et al.A reactive force field （ReaxFF） for zinc oxide.Surface Science，2008，602（5）:1020

19 Kulkarni A J，Zhou M，Ke F J.Orientation and size dependence of the elastic properties of zinc oxide nanobelts.Nanotechnology， 2005，16（12）:2749

20 Hu J，Liu X W，Pan B C.A study of the size-dependent elastic properties of ZnO nanowires and nanotubes.Nanotechnology， 2008，19（28）:285710

21 Momeni K.EnhancedmechanicalpropertiesofZnO nanowire-reinforced nanocomposites:A size-scale effect.Acta Mechanica， 2014，225（9）:2549

22 Zaoui A，Sekkal W.Pressure-induced softening of shear modes in wurtzite ZnO:A theoretical study.Physical Review B， 2002，66（17）:553～562

23 Wang J， Kulkarni A J， Ke F J， et al.Novel mechanical behavior of ZnO nanorods.Computer Methods in Applied Mechanics and Engineering， 2008， 197（41）:3182

24 Sun X W，Chua Y D，Song T，et al.Application of a shell model in molecular dynamics simulation to ZnO with zinc-blende cubic structure.Solid State Communications，2007，142（1）:15

25 曹樹謙，郭抗抗.壓電發(fā)電懸臂梁的非線性動力學建模及響應分析.動力學與控制學報，2014，12（1）:018（Cao S Q， Guo K K.Nonlinear modeling and analysis of piezoelectric cantilever energy harvester.Journal of Dynamics and Control， 2014，12（1）:018 （in Chinese））

26 彭劍，李祿欣，馬建軍.時滯反饋作用下壓電梁的參數(shù)共振分析.動力學與控制學報，2016，14（5）:412（Peng J， Li L X， Ma J J.Parametric resonanceof piezoelectric beam with time-delayed feedback.Journal of Dynamics and Control， 2016，14（5）:412 （in Chinese））

27 Ekinci K L， Roukes M L， Nanoelectromechanical systems.Review of Scientific Instruments， 2005，76（6）:061101

28 Craighead H G.Nanoelectromechanical systems.Science，2000，290:1532

29 Liu R M，Wang L F.Thermal vibration of a single-walled carbon nanotube predicted by semiquantum molecular dynamics.Physical Chemistry Chemical Physics， 2015，17:5194

30 Ke L L，Wang Y S.Thermoelectric-mechanical vibration of piezoelectric namobeams based on the nonlocal theory.Smart Materials and Structures， 2012，21（2）:025018

31 Chen C Q，Shi Y，Zhang Y S，et al.Size dependence of Young′s modulus in ZnO nanowires.Physical Review Letters， 2006，96:075505

32 Momeni K，Odegard G M，Yassar R S.Finite size effect on the piezoelectric properties of ZnO nanobelts:A molecular dynamics approach.Acta Materialia， 2012，60（13-14）:5117-5124

33 Binks D J.Computational modelling of zinc oxide and related oxide ceramics.The Department of Chemistry of the University of Surrey，1994

34 Dai S X，Dunn M L，Park H S.Piezoelectric constants for ZnO calculated using classical polarizable core-shell potentials.Nanotechnology， 2010，21（44）:445707

35 Wolf D，Keblinski P，Phillpot S R，et al.Exact method for the simulation of Coulombic systems by spherically truncated，pairwise r-1summation.Journal of Chemical Physics， 1999，110（17）:8254

36 Gdoutos E E，Agrawal R，Espinosa H D.Comparison of the Ewald and Wolf methods for modeling electrostatic interactions in nanowires.International Journal for Numerical Methods in Engineering， 2010，84（13）:1541

37 ANSI／IEEE Standard 176-1987

38 Wang L F， Hu H Y.Size effectson the effective Young′s modulus of nano crystal copper wires.International Journal of Computational Method， 2005，2（3）:315～326

39 Nose S.A molecular dynamics method for simulations in the canonical ensemble.Molecular Physics， 1984，52:255

40 Hoover W.Canonical dynamics:Equilibrium phase-space distributions.Physical Review A， 1985，31:1695

41 Huang T C.The effect of rotatory inertia and of shear deformation on the frequency and normal mode equations of uniform beams with simple end conditions.Journal of Applied Mechanics， 1961，28（4）:579～584

VIBRATION OF ZINC OXIDE NANOWIRES

Liu Saisai Wang Lifeng?
（State Key Laboratory of Mechanics and Control of Mechanical Structures，Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China）

The vibration of Zinc Oxide（ZnO） nanowires is studied via molecular dynamics（MD） simulation and continuum theory.The size effect of equivalent Young′s modulus and piezoelectric constant for the ZnO nanowires are described by core-shell model.The equivalent tensile Young′s modulus of ZnO nanowires in polarization direction increases gradually with the increase of the cross section size.The equivalent tensile Young′s modulus predicted by continuum theory is in a good agreement with the MD result.The equivalent bending Young′s modulus also increases with the increasing cross section size.Meanwhile， the piezoelectric constants of ZnO nanowires are larger than that of piezoelectric ceramics.The piezoelectric constants of ZnO nanowires decrease with the rising of the cross section size.In addition，the vibration of the cantilevered nanobeam made of ZnO is simulated by MD.The vibration frequencies of a ZnO nanowire keep constant at different temperatures.When the cross section size becomes larger，the vibration frequencies predicted by continuum theory get closer to those obtained by MD well.

ZnO nanowires， molecular dynamics， scale effect， piezoelectric effect， vibration

1 October 2016，revised 22 November 2016.

10.6052／1672-6553-2017-063

2016-10-01收到第1稿，2016-11-22收到修改稿.

?通訊作者 E-mail:walfe＠nuaa.edu.cn

? Corresponding author E-mail:walfe＠nuaa.edu.cn