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

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

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

0　引言

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

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

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

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

(1)

(2)

圖1　非平衡模擬過(guò)程

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

(3)

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

2　模擬方法

2.1　系綜

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

2.2　作用勢(shì)

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

(4)

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

2.3　局域溫度的量子修正

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

(5)

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

(6)

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

ω=vq

(7)

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

(8)

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

(9)

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

3　結(jié)果與討論

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

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

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

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

表1 為溫度在300K時(shí)氧化鋅薄膜的熱導(dǎo)率的計(jì)算結(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之間時(shí)，熱導(dǎo)率隨薄膜厚度的增加而增加，與大體積的氧化鋅熱導(dǎo)率相比，具有明顯的尺寸效應(yīng).

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

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

(10)

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

(11)

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

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

(12)

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

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

4　結(jié)論

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

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

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

(4)用氣動(dòng)理論對(duì)研究結(jié)果進(jìn)行分析，觀察到氧化鋅熱導(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)率的分子動(dòng)力學(xué)模擬[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2006,38(5)：708-711.

[3]孔憲仁，吳國(guó)強(qiáng)，孫兆偉，等.單晶碳和鍺薄膜熱導(dǎo)率的分子動(dòng)力學(xué)模擬[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，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)率的測(cè)試與分子動(dòng)力學(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)率的分子動(dòng)力學(xué)模擬[J].物理學(xué)報(bào)，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.