童 宇, 易法軍， 劉立勝， 張清杰

(武漢理工大學(xué) a.材料復(fù)合新技術(shù)國家重點實驗室； b.理學(xué)院， 湖北 武漢 430070)

熱電材料是一種將熱能和電能直接轉(zhuǎn)換的功能材料，在溫差發(fā)電和熱電致冷等領(lǐng)域具有極為重要的應(yīng)用前景。Bi2Te3基化合物是室溫附近性能最好的熱電材料。隨著納米研究技術(shù)的發(fā)展，近年來有關(guān)低維Bi2Te3材料取得高熱電優(yōu)值的報道不斷出現(xiàn)[1,2]。由于薄膜材料的熱電性能與塊體材料有所不同，Bi2Te3薄膜有較高的研究價值。人們在研究Bi2Te3納米薄膜材料的過程中探索了包括電化學(xué)沉積法在內(nèi)的諸多合成制備方法，并廣泛地關(guān)注其熱學(xué)和電學(xué)特性[3～5]。目前對熱電材料的研究仍主要集中在熱電性能，但是材料的力學(xué)性能對其服役行為與使用壽命具有重要的影響，然而對于Bi2Te3薄膜材料力學(xué)性能的研究幾乎處于空白狀態(tài)。另外由于Bi2Te3熱電材料在實際應(yīng)用中不可避免地要經(jīng)歷較大的溫度波動，研究材料的力學(xué)性能隨溫度的變化關(guān)系具有重要意義。因此對Bi2Te3納米薄膜材料力學(xué)性能的研究非常重要。但是由于納米尺度實驗本身存在一定困難，而分子動力學(xué)方法能從微觀研究材料原子運動過程和結(jié)構(gòu)力學(xué)變形細(xì)節(jié)，另外由于分子動力學(xué)方法能夠考慮溫度效應(yīng)，有助于人們對不同溫度下材料的力學(xué)行為進行預(yù)測。因此從原子、分子尺度來考察材料的結(jié)構(gòu)和特性的分子動力學(xué)模擬方法成為一種重要的研究手段，并取得了許多重要的進展。例如Yang等人利用MAEAM勢計算了Nb(110)薄膜的熔點、熔化機制以及相應(yīng)的動力學(xué)行為[6]。Cao等人利用分子動力學(xué)研究了不同厚度、取向和加載方向?qū)nO納米薄膜彈性性能的影響[7]。

本文用分子動力學(xué)方法模擬了單晶Bi2Te3薄膜沿其力學(xué)性能較弱的c軸方向的單軸拉伸過程。首先通過對系統(tǒng)的弛豫獲得了無初應(yīng)力以及能量最小化的穩(wěn)定狀態(tài)模型。在此基礎(chǔ)上考察了微觀尺度上薄膜材料整體的變形演化過程，得到了拉伸過程中能量、應(yīng)力的變化曲線。模擬不同溫度下的拉伸，分析溫度效應(yīng)對其力學(xué)性能的影響。當(dāng)前研究結(jié)果為研究Bi2Te3納米薄膜材料器件提供了有效的依據(jù)。

1 分子動力學(xué)模擬方法

分子動力學(xué)模擬是計算經(jīng)典多粒子體系的一種有效方法，計算整體系統(tǒng)的力學(xué)行為和性能，是納米計算力學(xué)的重要工具。模擬系統(tǒng)中的粒子運動符合經(jīng)典的牛頓方程，即

(1)

其中原子i所受的力可以直接由勢能函數(shù)U對坐標(biāo)ri的一階導(dǎo)數(shù)得到。因此原子間相互作用的描述方式是分子動力學(xué)模擬計算結(jié)果精確與否的關(guān)鍵。本文選用的Huang和Kaviany等人發(fā)表的Bi2Te3的多體勢。該勢函數(shù)最初用于計算Bi2Te3的晶格熱導(dǎo)率[8]，勢參數(shù)由第一性原理計算以及實驗參數(shù)擬合而來。利用該勢函數(shù)我們曾計算了Bi2Te3材料從0到600 K溫度范圍內(nèi)的力學(xué)性質(zhì)[9,10]。Bi2Te3多體勢函數(shù)的具體形式為：

Etotal=Eelectr+Ebond+Eangle

(2)

其中Eelectr為靜電勢能，表達式為：

(3)

Ebond為Morse型鍵伸縮勢能，表達式為：

Ebond=D{[1-exp(-α(r-r0))]2-1}

(4)

Eangle為三體鍵角彎曲勢能，表達式為：

Eangle=K[cos(θ)-cos(θ0)]2

(5)

勢函數(shù)具體參數(shù)詳見參考文獻[8]。

圖1 Bi2Te3晶體結(jié)構(gòu)與坐標(biāo)系

圖2 薄膜初始模型

模擬之前，先將Bi2Te3的納米薄膜模型置于設(shè)定的溫度中，并按照Maxwell-Boltzmann分布對原子賦予初始速度。然后將系統(tǒng)弛豫100皮秒(ps)，在此過程中利用Berendsen恒壓算法使x和z方向的壓力保持為0 Pa，同時利用Nose-Hoover熱浴算法使系統(tǒng)的溫度保持在特定值，以使系統(tǒng)初始應(yīng)力為0 Pa。模擬時間步長設(shè)為1飛秒(fs)。弛豫過程完成后進行沿z軸單軸拉伸，在此過程中Berendsen壓力控制僅限制在x方向。為了模擬真實的拉伸試驗，在拉伸過程中模型一端固定而另一端原子沿z軸即[0 0 1]晶向運動。拉伸應(yīng)變率設(shè)定為1×107s-1。拉伸過程中每施加0.001的應(yīng)變系統(tǒng)弛豫2000步，以使系統(tǒng)原子回到平衡態(tài)，然后再施加應(yīng)變，如此重復(fù)施加應(yīng)變-弛豫的過程，直到模型發(fā)生破壞。本文中采用Verlet算法的速度形式對系統(tǒng)運動方程進行力學(xué)時間積分計算。模擬過程中每隔5000步輸出原子的位置、應(yīng)力分量、總能量、勢能、動能等參數(shù)。模擬過程中應(yīng)力的計算方法采用Virial算法。論文中所有模擬過程均利用LAMMPS分子動力學(xué)軟件完成[11]，原子構(gòu)型的顯示使用VMD軟件[12]。

2 模擬結(jié)果與討論

拉伸模擬之前，先對系統(tǒng)進行一個趨衡過程即弛豫過程。這一步驟可以使系統(tǒng)達到初始平衡狀態(tài)，同時消除初始速度帶來的隨機因素，以及自由表面的引入對納米薄膜的影響。圖3表示了Bi2Te3納米薄膜在300 K時系統(tǒng)總勢能、應(yīng)力以及溫度隨著馳豫時間步的變化曲線。從圖3(a)可以看出，在弛豫初期，系統(tǒng)的總能曲線經(jīng)過一段時間的馳豫之后，達到一個穩(wěn)定位置。圖3(b)顯示在初始狀態(tài)，系統(tǒng)的初應(yīng)力很大，約為-1.5 GPa。經(jīng)過一段時間的弛豫后，應(yīng)力值保持在零點附近，說明弛豫后薄膜整體是一個無外力的平衡體系。圖3(c)表示溫度經(jīng)過多次震蕩后，在4萬步左右進入一個相對穩(wěn)定的狀態(tài)，達到所設(shè)定的值。以上3幅圖可以看出，系統(tǒng)在經(jīng)歷馳豫過程之后達到平衡態(tài)。

圖3 弛豫過程中

圖4為拉伸過程中的原子構(gòu)型?？梢娎熳冃纬跗冢珺i2Te3納米薄膜的原子排列規(guī)則(圖4(a)，(b))。隨著應(yīng)變的逐漸增加，沒有出現(xiàn)比較明顯的原子滑移或錯位現(xiàn)象。這種原子整齊排列的狀態(tài)一直持續(xù)在整個拉伸過程，直到應(yīng)變值達到0.048后，薄膜沿相鄰Te1-Te1層間(0 0 1)晶面突然發(fā)生斷裂并分離，斷裂截面平整(圖4(c))。這種斷裂形式與塊體相似。根據(jù)Drabble等人提出Bi2Te3的成鍵模型[13]。相鄰層Te1-Te1之間的相互作用比較弱，幾乎沒有電荷作用。Te1原子的價電子都只與近鄰的Bi原子成鍵，所以Te1-Te1之間表現(xiàn)為較弱的范德華相互作用，比離子鍵或者共價鍵弱(但比惰性氣體如Kr-Kr或Xe-Xe的范德華鍵能略高)[8]。因此Bi2Te3晶體極易沿此晶面發(fā)生破壞。

圖4 Bi2Te3薄膜c軸方向單軸拉伸結(jié)構(gòu)演化

圖5 原子能量-應(yīng)變曲線

圖6 應(yīng)力-應(yīng)變曲線

圖5和圖6分別為Bi2Te3納米薄膜300 K時沿c軸方向拉伸的系統(tǒng)能量-應(yīng)變和應(yīng)力-應(yīng)變曲線。由于在分子動力學(xué)模擬過程中，系統(tǒng)能量的變化趨勢可以體現(xiàn)材料變形的特征，因此結(jié)合系統(tǒng)能量的變化曲線進行說明。由圖6可以發(fā)現(xiàn)，初始階段的應(yīng)力應(yīng)變曲線為線性變化，應(yīng)力隨著應(yīng)變的增加而增加，模型中的原子逐漸偏離原來的理想位置，對應(yīng)系統(tǒng)總能量也逐漸上升。隨著拉伸應(yīng)變的增大，原子間距逐漸增大，原子的運動也逐漸加劇，原子能量上升的速度也隨之增大，從而導(dǎo)致體系能量急劇升高(圖5)。當(dāng)應(yīng)變值達到0.03后，應(yīng)力應(yīng)變曲線變得略微平坦，表現(xiàn)出非線性特點。當(dāng)應(yīng)變達0.048后，隨著應(yīng)變的進一步增加，應(yīng)力急劇下降到0，說明此時Bi2Te3納米薄膜已經(jīng)破壞。由能量-應(yīng)變曲線可以看出，由于原子間相互作用此時急劇減小，系統(tǒng)勢能得到釋放，導(dǎo)致總能量降低。由應(yīng)力-應(yīng)變曲線以及拉伸過程中的原子構(gòu)型變化可以看出Bi2Te3納米薄膜沿c軸方向拉伸的破壞形式表現(xiàn)為明顯的脆性斷裂的特征。

根據(jù)以前對Bi2Te3薄膜沿a軸方向的研究[10]，比較兩個主軸方向的Bi2Te3薄膜拉伸過程，發(fā)現(xiàn)兩方向上共同的特點是拉伸初始階段應(yīng)力應(yīng)變幾乎呈線性變化，表現(xiàn)為脆性斷裂。但是沿c軸的斷裂由于直接發(fā)生在作用較弱的相鄰層Te1原子之間，因而從強度值上兩者相差較大，表明Bi2Te3納米薄膜具有明顯的各向異性的特點。

圖7為0～600 K溫度下Bi2Te3納米薄膜沿c軸方向拉伸的應(yīng)力應(yīng)變曲線?？梢婋S著溫度的升高應(yīng)力應(yīng)變曲線的斜率逐漸降低。極限強度和破壞應(yīng)變也隨之降低。其中擬合的彈性模量(應(yīng)變?yōu)?.01)從0K時的42 GPa下降到 600 K時的31 GPa，下降幅度為26%。極限強度從0K時的1.5 GPa下降到 600 K時的1.04 GPa，下降幅度為30%。破壞應(yīng)變從0 K時的0.051下降到 600 K時的0.043，下降幅度為15%。

圖7 不同溫度下的應(yīng)力-應(yīng)變曲線

3 結(jié) 語

本文根據(jù)Huang和Kaviany提出的Bi2Te3的多體勢函數(shù)利用分子動力學(xué)方法對Bi2Te3納米薄膜沿著力學(xué)性能較差的c軸方向進行單軸拉伸模擬。拉伸之前對系統(tǒng)進行弛豫，考察了弛豫過程中系統(tǒng)能量、溫度和初始應(yīng)力的變化。根據(jù)拉伸過程中原子構(gòu)型發(fā)現(xiàn)此方向的斷裂是沿著晶體中結(jié)合力較弱的Te1-Te1層斷開。由拉伸過程中的應(yīng)力-應(yīng)變和能量-應(yīng)變變化曲線發(fā)現(xiàn)Bi2Te3薄膜沿c軸方向拉伸的破壞形式表現(xiàn)為明顯的脆性斷裂的特征。受溫度的影響納米薄膜的彈性常數(shù)、極限強度與破壞應(yīng)變隨溫度升高而減小。

[1]Tang X F, Xie W J, Li H, et al. Preparation and thermoelectric transport properties of high-performance P-type Bi2Te3with layered nanostructure[J]. Applied Physics Letters, 2007, 90: 012102-012104.

[2]Xie W J, Tang X F, Yan Y G, et al. High thermoelectric performance BiSbTe alloy with unique low-dimensional structure[J]. Journal of Applied Physics, 2009, 105: 113713-113715.

[3]Walachova J, Zeipl R, Zelinka J, High room-temperature figure of merit of thin layers prepared by laser ablation from Bi2Te3target[J]. Applied Physics Letters, 2005, 87: 081902-081907.

[4]Miyazaki Y, Kajitani T. Preparation of Bi2Te3films by electrodeposition[J]. Journal of Crystal Growth, 2001, 229: 542-546.

[5]Kwon S D, Ju B K, Yoon S J, et al. Fabrication of bismuth telluride-based alloy thin film thermoelectric devices grown by metal organic chemical vapor deposition[J]. Journal of Electronic Materials, 2009, 38: 920-924.

[6]Yang X Y, Wu D. The melting behaviors of the Nb(1 1 0) nanofilm: a molecular dynamics study[J]. Applied Surface Science, 2010, 256: 3197-3203.

[7]Cao G X, Chen X. Size dependence and orientation dependence of elastic properties of Zno nanofilms[J]. International Journal of Solids and Structures, 2008, 45: 1730-1753.

[8]Huang B L, Kaviany M. Ab initio and molecular dynamics predictions for electron and phonon transport in bismuth telluride[J]. Physical Review B, 2008, 77: 125209-125223.

[9]Tong Y, Yi F J, Liu L S, et al. Molecular dynamics study on thermo-mechanical properties of bismuth telluride bulk[J]. Computational Materials Science, 2010, 48: 343-348.

[10]Tong Y, Yi F J, Liu L S, et al. Molecular dynamics study of mechanical properties of bismuth telluride nanofilm[J]. Physica B, 2010, 405: 3190-3194.

[11]Plimpton S J. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics, 1995, 117: 1-19.

[12]Humphrey W, Dalke A, Schulten K. VMD-visual molecular dynamics[J]. Journal of Molecular Graphics and Modelling, 1996, 14: 33-38.

[13]Drabble J R, Goodman C H. X-ray spectroscopic investigation of bismuth telluride[J]. Journal of Physics and Chemistry of Solids, 1958, 5: 142-147.