白素媛，王立凡，支明蕾，王　斌，賈孟晗

(遼寧師范大學(xué))

*遼寧省自然科學(xué)基金(2015020079)

0　引言

半導(dǎo)體材料是電子信息技術(shù)發(fā)展的基礎(chǔ)，是人們一直研究、關(guān)注的焦點.伴隨著第一代、第二代、第三代半導(dǎo)體技術(shù)的快速發(fā)展，半導(dǎo)體器件的尺寸已經(jīng)進入到了微米、納米尺度.又由于物質(zhì)的局限性、量子效應(yīng)、巨大的表面及界面效應(yīng)，使微米、納米尺度下的器件的熱物性產(chǎn)生了明顯的尺寸效應(yīng)，因此器件的熱學(xué)性質(zhì)也開始受到了人們的廣泛關(guān)注.由于薄膜材料的熱導(dǎo)率會直接影響微納米器件的散熱性能,進而會間接影響器件的運行速度.因此,提高薄膜材料的熱導(dǎo)率[1],使微納米器件始終在最佳的溫度下運行則是半導(dǎo)體電子工業(yè)要解決的一個難題. 在微尺度傳熱的領(lǐng)域中，數(shù)值模擬發(fā)揮著很大的作用，它可以對實驗過程中遇到的不足之處進行彌補.數(shù)值模擬方法有很多種，其中應(yīng)用最廣泛的為分子動力學(xué)模擬方法(MD).王亞輝[2]、孫憲仁[3]、 Kato R[4]等均運用分子動力學(xué)方法分別對金剛石薄膜、單晶碳和鍺薄膜、二氧化硅薄膜的熱導(dǎo)率進行了研究.氧化鋅[5]作為一種過渡金屬氧化物，具有良好的光電磁等特性.與其它半導(dǎo)體相比，其制作工藝簡單、成本較低、無毒、化學(xué)性質(zhì)穩(wěn)定、熱穩(wěn)定性好.在溫度為500℃以下就可以制備出氧化鋅薄膜，這比制備其他半導(dǎo)體材料所需要的溫度都要要低，這也是氧化鋅成為研究熱點的原因.氧化鋅晶體結(jié)構(gòu)隨著環(huán)境的改變會形成不同結(jié)構(gòu)的晶體，分別為立方閃鋅礦結(jié)構(gòu)、立方巖鹽結(jié)構(gòu)、六邊纖鋅礦結(jié)構(gòu)，其中纖鋅礦結(jié)構(gòu)最為常見，并且其熱力學(xué)性質(zhì)最為穩(wěn)定性，因此纖鋅礦型氧化鋅成為研究熱點.對于研究薄膜材料的熱物性，熱導(dǎo)率的大小將直接影響器件的耐高溫性能以及器件的散熱性能，所以，熱導(dǎo)率的優(yōu)化和求解具有重要的實際意義和研究價值.黃正興[6]利用實驗方法研究了六角纖鋅結(jié)構(gòu)氧化鋅薄膜厚度對熱導(dǎo)率的影響；Kulkami、Zhou[7]等利用分子動力學(xué)方法模擬氧化鋅納米帶沿z方向上的熱導(dǎo)率，研究結(jié)果顯示氧化鋅納米帶具有明顯的溫度效應(yīng)和尺寸效應(yīng)；楊平[8]等利用非平衡分子動力學(xué)方法對閃鋅結(jié)構(gòu)氧化鋅薄膜界面的熱導(dǎo)率進行了模擬，結(jié)果顯示氧化鋅薄膜界面的熱導(dǎo)率隨厚度的增加而增加.該文采用的是非平衡分子動力學(xué)方法對纖鋅礦型氧化鋅薄膜的熱導(dǎo)率進行了模擬，并且分析了纖鋅礦型氧化鋅薄膜的厚度對熱導(dǎo)率的影響.

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

分子動力學(xué)方法是根據(jù)經(jīng)典牛頓運動學(xué)原理、應(yīng)用力場的變化所發(fā)展起來的一種研究凝聚態(tài)物質(zhì)的模擬方法. 非平衡分子動力學(xué)方法、平衡分子動力學(xué)方法是分子動力學(xué)的兩種方法.平衡分子動力學(xué)方法(EMD)是應(yīng)用線性相應(yīng)理論和漲落耗散理論，運用相關(guān)函數(shù)求得熱導(dǎo)率.所以，這種方法的優(yōu)點是尺寸效應(yīng)很小.非平衡分子動力學(xué)方法(NEMD)[9-10]是通過對系統(tǒng)施加相應(yīng)的溫度梯度或者施加相應(yīng)的熱流，根據(jù)傅里葉熱導(dǎo)率定律求得熱導(dǎo)率.該方法的主要的優(yōu)點是計算速度快，缺點是具有顯著地尺寸效應(yīng).該文研究的纖鋅礦氧化鋅薄膜的厚度在15.593～31.2096nm范圍內(nèi)，其模擬的結(jié)果與實驗值相比具有明顯的尺寸效應(yīng)，因此該文采用非平衡分子動力學(xué)模擬方法進行研究.

在模擬過程中，根據(jù)熱流形成的原因，采用非平衡分子動力學(xué)中的交換速度法進行模擬.交換速度法是通過交換粒子的速度來實現(xiàn)外加熱流，即先給定熱流，再統(tǒng)計溫度梯度，最后利用傅立葉定律計算求得熱導(dǎo)率，其計算公式為：

(1)

(2)

圖1　非平衡模擬過程

其中，nk為第k層所含的原子數(shù)，kB為玻爾茲曼常數(shù)，mi為原子質(zhì)量，vi表示原子速度.在圖1中，將第1層規(guī)定為冷層，第N/2+1層規(guī)定為熱層，并且符合周期性邊界條件.粒子速度不斷的交換將會引起熱層更熱、冷層更冷以及冷熱層的溫差會導(dǎo)致系統(tǒng)產(chǎn)生一個與熱流施加方向相反的熱量，從而產(chǎn)生溫度梯度.當(dāng)兩個方向相反的能量在傳遞過程中達到動態(tài)平衡時，就會得到穩(wěn)定的溫度梯度響應(yīng)曲線，同時通過原子之間的速度交換，可以計算外部施加的熱流，也就是在整個速度交換過程中所傳遞的能量，從而體系的熱導(dǎo)率可以結(jié)合(3)式得出：

(3)

其中，t表示模擬過程中粒子交換速率的總時間，A表示模型的橫截面積，vh表示冷層中能量最高粒子的速率，vc表示熱層中能量最低粒子的速率.由于周期性邊界條件對熱流方向有影響，體系內(nèi)的熱量將會從熱層沿著兩端流向冷層，所以，分母除以2.

2　模擬方法

2.1　系綜

在分子動力學(xué)模擬方法中，對溫度和壓力的控制是運用系綜的核心問題，該文選用微正則系綜(NVE).在模擬過程中，體系粒子的數(shù)目保持不變，粒子的體積保持不變，體系粒子間的能量保持守恒，在體系演化過程中，粒子是沿著相空間中的恒定能量軌道進行演化.并且，在整個模擬過程中整個體系內(nèi)的能量會一直不變，因此不需要重新控制體系的能量在模擬計算中首先將系統(tǒng)調(diào)節(jié)到特定的能量，同時給出初始溫度，然后再進行調(diào)整，最后直到體系的能量達到指定的能量值.微正則系綜(NVE)中，一般通過對溫度的標(biāo)定來實現(xiàn)對能量的調(diào)整.

2.2　作用勢

在分子動力學(xué)(MD)模擬中，用勢函數(shù)來描述分子之間、原子之間的相互作用力是動力學(xué)模擬的重要環(huán)節(jié)，勢函數(shù)的選取會直接影響模擬結(jié)果的精確性.常用的勢函數(shù)分為多體勢和兩體勢.多體勢主要為嵌入原子勢；兩體勢主要有Morse勢函數(shù)、Born-Mayer勢函數(shù)、Lennard-Jones勢函數(shù)、Buckingham勢函數(shù).該文采用Buckingham勢函數(shù)來描述氧化鋅薄膜內(nèi)部原子之間的作用力，Buckingham勢函數(shù)表達形式如下：

(4)

其中Eij表示原子i和原子j的相互作用能.

2.3　局域溫度的量子修正

在分子動力學(xué)(MD)模擬中，由Boltzmann能量均分定理可以計算出系統(tǒng)的溫度，其表達式如下：

(5)

其中N是模擬的粒子總數(shù)， kB是玻爾茲曼常數(shù),TMD是系統(tǒng)模擬得到的溫度，mi是粒子的質(zhì)量，vi是粒子的速度.當(dāng)體系的溫度遠低于材料的Debye溫度時，由于體系的比熱容與其溫度有關(guān)，所以(5)式不再成立，也就是說系統(tǒng)實際的溫度T≠TMD.為了得到精確的結(jié)果，務(wù)必對(5)式中的模擬溫度TMD、模擬結(jié)果kMD進行量子化修正.采用Debye模型進行系統(tǒng)溫度統(tǒng)計的量子修正，修正積分方程如下:

(6)

通過求解(6)式，可以得到和局域溫度TJ,MD相對應(yīng)的晶格溫度Tj.該文選擇的Debye模型為:

ω=vq

(7)

Debye模型所給出的晶格熱容為:

(8)

其中Cv是質(zhì)量定容熱容，ΘD是Debye溫度.在Debye近似下，能量方程(6)轉(zhuǎn)化為:

(9)

該文通過求解方程(9)，實現(xiàn)了局域溫度的量子修正.

3　結(jié)果與討論

在半導(dǎo)體材料中，氧化鋅是一種過渡金屬氧化物、Ⅱ-Ⅵ 族寬禁帶半導(dǎo)體材料，具有較高的激子結(jié)合能，其化學(xué)性質(zhì)和熱物性非常穩(wěn)定，既可以工作在室溫下，又可以工作在高溫下，其制備溫度比其它寬禁帶半導(dǎo)體都要低. 由于在低維條件下，氧化鋅的熱流密度和電場強度都很大，所以其熱物性變得尤其重要.氧化鋅晶體結(jié)構(gòu)[11]如圖2所示，其參數(shù)分別為：禁帶寬度大約為3.37eV、熔點為1975 ℃，晶胞參數(shù)a=3.289 ?，c=5.307 ?.該文采用的是非平衡分子動力學(xué)模擬方法，模擬系統(tǒng)在xy截面上采用4×4個晶胞，在z方向上取20、24、28、32、36、40個基本單元，系統(tǒng)溫度為300K，系統(tǒng)的時間步為1.1fs，所模擬的計算步數(shù)為150萬步.

圖2　纖鋅礦結(jié)構(gòu)氧化鋅晶體結(jié)構(gòu)

以4×4×96系統(tǒng)，薄膜厚度為24.963nm為例，在經(jīng)過150萬步計算之后得到如圖3所示的溫度分布曲線，根據(jù)圖像可以看出，各層的溫度分布曲線比較光滑，黃色原點為高低控溫層，由于熱流將從控溫層通過，所以控溫區(qū)內(nèi)存在溫差.

圖3　厚度為24.963 nm的氧化鋅薄膜中各層溫度分布情況

表1 為溫度在300K時氧化鋅薄膜的熱導(dǎo)率的計算結(jié)果，由表可知，薄膜厚度在15.593～31.2096nm范圍內(nèi)的氧化鋅薄膜，其熱導(dǎo)率在0.87～1.52165W/(m·K)之間.圖4給出了氧化鋅薄膜熱導(dǎo)率與厚度之間的關(guān)系.由圖可知，當(dāng)模擬的薄膜厚度在15.593～31.2096nm之間時，熱導(dǎo)率隨薄膜厚度的增加而增加，與大體積的氧化鋅熱導(dǎo)率相比，具有明顯的尺寸效應(yīng).

表1　氧化鋅薄膜熱導(dǎo)率的計算結(jié)果

在半導(dǎo)體材料中，熱傳導(dǎo)是通過聲子之間的相互作用完成的，并且聲子的運動決定了材料的傳熱特性.根據(jù)聲子氣的動力論[12]，熱導(dǎo)率可以表示為：

(10)

其中，v表示聲速，l表示聲子的平均自由程.根據(jù)固體物理理論，當(dāng)發(fā)生尺寸效應(yīng)時薄膜的厚度必將小于其聲子的平均自由程.根據(jù)該文所模擬的溫度，可用下式近似表示其聲子平均自由程：

(11)

其中， γ是Gruneisen常數(shù)，通常情況下γ≈2，Tm是材料的熔點溫度，a是晶格常數(shù).該文模擬的氧化鋅材料，可以根據(jù)(11)式估算其平均聲子自由程，并且其估算值要大于該文所模擬的薄膜厚度.所以，整個薄膜界面顯示出尺寸效應(yīng).

另外，由于薄膜邊界對聲子的散射也會影響熱導(dǎo)率的大小，因此，可以根據(jù)Matthiessn法則來確定薄膜中有效的聲子平均弛豫時間，即：

(12)

其中， τb是弛豫時間， p是薄膜邊界對聲子的鏡反射率，d是薄膜的厚度.根據(jù)(12)式以及晶格導(dǎo)熱系數(shù)可以得出，當(dāng)薄膜的厚度和薄膜邊界對聲子的鏡發(fā)射率越大時，其有效弛豫時間也越大，從而薄膜的熱導(dǎo)率也會越大，反之則越小.這就是薄膜的熱導(dǎo)率隨著薄膜的厚度增加而增大的原因.

文獻[5]中運用實驗方法，即通過組建瞬態(tài)熱反射測試系統(tǒng)對采用溶膠凝膠法的氧化鋅薄膜的熱導(dǎo)率進行了測試.研究溫度為500 K、薄膜厚度在80～276 nm之間，最后得出熱導(dǎo)率的范圍在1.4～6.51 W/(m·K)之間，研究結(jié)果表明熱導(dǎo)率隨氧化鋅薄膜厚度的增加而增加.與該文所得到的結(jié)果相類似，可見該文的研究是可信的和合理的，至于還存在差異，是由于采用的方法不同而導(dǎo)致的.

4　結(jié)論

( 1)采用微正則系綜，利用Buckingham勢函數(shù)，運用非平衡分子動力學(xué)方法計算氧化鋅薄膜的熱導(dǎo)率.

(2)模擬溫度為300K，模擬的氧化鋅薄膜厚度在15.593～31.2096 nm之間，其熱導(dǎo)率計算結(jié)果范圍在0.870575～1.52165 W/(m·K).

(3)當(dāng)系統(tǒng)溫度一定時，氧化鋅薄膜的熱導(dǎo)率隨著薄膜的厚度的增加而增加.

(4)用氣動理論對研究結(jié)果進行分析，觀察到氧化鋅熱導(dǎo)率的模擬結(jié)果具有明顯的尺寸效應(yīng).

[1]Afshin Ahmadi Nadooshan. An experimental correlation approach for predicting thermal conductivity of water-EG based nanofluids of zinc oxide[J]. Physica E,2017,87:15-19.

[2]王亞輝，劉林華，孔憲仁.金剛石薄膜熱導(dǎo)率的分子動力學(xué)模擬[J].哈爾濱工業(yè)大學(xué)學(xué)報，2006,38(5)：708-711.

[3]孔憲仁，吳國強，孫兆偉，等.單晶碳和鍺薄膜熱導(dǎo)率的分子動力學(xué)模擬[J].哈爾濱工業(yè)大學(xué)學(xué)報，2006,38(4)：517-519.

[4]Kato R. Conductivity Measurement of Thermally-Oxidized Si Films on Silicon Wafer Using a Thermo-Reflectance Technique[J].International Journal of Thermophysics,2005,26(1):179-190.

[5]宋淑芳.納米ZnO薄膜的特性研究[D].內(nèi)蒙古大學(xué).

[6]黃正興.薄膜熱導(dǎo)率的測試與分子動力學(xué)模擬研究[D].大連理工大學(xué).

[7]Ambarish J Kulkami, Zhou Min.Size-dependent thermal conductivity of zinc oxide nanobelts[J].APPLIED PHYSICS LETTERS,2006,88(14):19211-19213.

[8]楊平，吳勇勝，許海鋒，等. 納米薄膜界面熱導(dǎo)率的分子動力學(xué)模擬[J].物理學(xué)報，2011,60(6)：066601(1-7).

[9]Wang Zhiguo, Gao Fei, Crocombette Jean-Paul. Thermal conductivity of GaN nanotubes simulated by nanequilibrium molecular dynamics[J].PHYSICAL REVEW B,2007,75(15):33031-33034.

[10] Yasushi Ichikawa,Yasuaki Hiwatari.Thermal conductivity of uranium dioxide by nonequilibrium molecular dynamics simulation[J].PHYSICAL REVIEW B,1991,60(1): 292-298.

[11] Chen Mingyuan, Hong Zhenghan, Lin Shiangjiun,et al. Buckling instability of zinc oxide nanobelts under uniaxial compression investigated using molecular dynamics[J]. Computational Materials Science,2014,85:217-222.

[12] Toru Yoshikawa,Takashi Yagi,Nobuto Oka,et al.Thermal Conductivity of Amorphous Indium-Gallium-Zinc Oxide Thin Films[J].Applied Physics Express,2013,6(2):11011-11013.