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

(遼寧師范大學)

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

0　引言

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

1　分子動力學方法

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

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

(1)

(2)

圖1　非平衡模擬過程

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

(3)

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

2　模擬方法

2.1　系綜

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

2.2　作用勢

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

(4)

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

2.3　局域溫度的量子修正

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

(5)

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

(6)

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

ω=vq

(7)

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

(8)

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

(9)

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

3　結果與討論

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

圖2　纖鋅礦結構氧化鋅晶體結構

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

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

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

表1　氧化鋅薄膜熱導率的計算結果

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

(10)

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

(11)

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

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

(12)

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

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

4　結論

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

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

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

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

[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]王亞輝，劉林華，孔憲仁.金剛石薄膜熱導率的分子動力學模擬[J].哈爾濱工業(yè)大學學報，2006,38(5)：708-711.

[3]孔憲仁，吳國強，孫兆偉，等.單晶碳和鍺薄膜熱導率的分子動力學模擬[J].哈爾濱工業(yè)大學學報，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].內蒙古大學.

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

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

[8]楊平，吳勇勝，許海鋒，等. 納米薄膜界面熱導率的分子動力學模擬[J].物理學報，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.