分子动力学模拟铜薄膜热导率

第19卷第4期2002年7月

计　　算　　物　　理

　CHI NESE JOURNA L OF C OMPUT ATI ONA L PHY SICS 　

V ol. 19,N o. 4　

[文章编号]　10012246X (2002) 0420321204

分子动力学模拟铜薄膜的热导率

叶杰成, 　唐祯安

(大连理工大学电子工程系, 辽宁大连　116023)

[摘　要]　采用分子动力学(M D ) 方法模拟铜薄膜的热导率, ～100～600K 范围内

铜薄膜热导率对尺寸及温度的依赖关系.

[关键词]　薄膜; 热导率; 分子动力学[中图分类号]　TK 123; O41113

[0　引言

分子动力学计算机模拟是研究微尺度热传输特

[1,2]

性的有力工具. 分子动力学是一种决定论(deter 2ministic ) 的原子级数值实验方法, 它不仅可以得到系

度梯度, 以至于无法准确确定系统的温度, 从而不能

正确地得到热导率和温度的关系;2) 由于系统不是均质的, 就不便采用周期性边界条件, 这就造成所模拟的系统粒子数目太多因而计算量巨大. 而采用作为分子动力学模拟基础的E MD 方法, 只需计算出平衡态下的热流函数J , 然后利用G reen 2K ubo 公式和统计物理规律即可求得热导率

[6,7]

统中每个粒子的运动轨迹, 还可以提供系统的时间演化过程信息, 进而可以通过统计力学得知系统的静态(热力学) 和动态(传输) 特性. 对一些理论上难以说明或实验中难以观察的现象给出基于物质微观结构及分子动力学关系的解释, 特别适合研究薄膜材料的微尺度热传导特性.

在微尺度热传输领域中, 由于尺度效应, 材料的热物性参数与宏观体形材料相比有很大的不同, 并存在着很多的未知量表明

[1]

[1]

.

本文采用E MD 方法对厚度在100～400nm 、温度在100～600K 范围内的铜薄膜热导率作了分子动力学模拟.

1　分子动力学模拟的基本原理和方法

分子动力学模拟的基本原理是建立一个合适的粒子系统, 利用牛顿运动方程确定系统中粒子的运动, 通过求得粒子动力学方程组的数值解, 决定系统中各个粒子在相空间中的运动规律, 然后按照统计物理和热力学原理得出系统相应的宏观物理特性.

对于系统温度保持不变的正则系综(NVT ) , 描述该系统的运动方程如下

r i =p i Πm i ,

p i =f i -ξ(r i , p i ) ×p i ,

[8]

, 热导率就是其中之一. 研究

, 薄膜材料的热导率是薄膜厚度的函数. 但实

[1,2]

验测定薄膜材料热导率, 不仅对设备要求很高, 技术上也有相当大的难度

. 因此分子动力学模拟在确

定微尺度热物性参数方面发挥着重要的作用. 本文通过分子动力学方法, 模拟计算铜薄膜的热导率. 从统计物理的角度可以将分子动力学分为平衡

态(E ) 和非平衡态(NE ) 两种模拟方法. 在以往的一些有关分子动力学模拟方法文献中, 都是用非平衡

[3～5]

态分子动力学(NE MD ) 方法来计算热导率. NE MD 是通过对系统施加扰动, 根据扰动所产生的

(1) (2) (3)

r i 、f i 、p i 、U (r 1, r 2, …, r n ) 和m i 分别是第i 个粒子

响应来求解有关的特性. 该方法主要用于处理远离平衡态或非均质材料, 但存在如下缺点:1) 为了得到足够大的温度梯度, 往往需要对系统施加很大的温

[收稿日期]2000-08-07; [修回日期]2001-02-19[基金项目]国家自然科学基金(5999555025) 资助项目

的坐标、所受到的合力、动量、其它粒子对其的作用势和质量, 其中ξ是约束温度不变的参数. 由于温度不变, 系统动能E k 必须满足 E k =0, 所以参数ξ是

[作者简介]叶杰成(1975-) , 男, 浙江永嘉, 硕士生, 从事微尺度热传输数值模拟方面的研究.

Δ

f i =-

i

U (r 1, r 2, …, r n ) .

322计　　算　　物　　理

第19卷　

ξ=

p ∑

∑|

i

i

×f i

p i |

2

其中V 为系统的体积:k B 为波耳兹曼常数, T 为系

.

(4)

统的温度, J 为热流函数, 可通过E MD 方法确定. 在

[5]

笛卡尔坐标下, 热流函数的微观表达式如下

J

α

在分子动力学模拟中, 关于材料系统内的粒子间相互作用势函数U (r 1, r 2, …, r n ) 的选取非常重要, 它决定了模拟结果正确与否. 本文采用的势函数是嵌入原子法(Embedded Atom Method , E AM ) 多体势函数

[9]

=

V

∑

j

E j v j -4

α

∑∑

j

r jk r jk (v j +v j ) .

r jk 9r k (≠j )

α

β

αβ

(11)

α、β=x , y , z , E j 为粒子, v j 为粒子j 的速度, 为粒子j , V 为系统的体积. E AM 是建立在H onhenberg 2K ohn 理论基础上

的, 该理论采用局域电子密度函数描述相互作用的原子中的一个原子对系统能量的贡献

U i =-f (ρi ) +

2

V ij

((6)

图1是厚度从100～400nm 的几种铜薄膜在不

同温度下的热导率. 从图中可以看出在某一温度下, 热导率随着膜厚的减小而减小. 热导率不是随着膜厚均匀地减小, 而是膜越薄, 热导率递减得越厉害. 由图1还可见, 在同一厚度, 热导率随温度增加而减少, 其原因是高温时晶格对声子的散射加强. 表2说明在温度为300K 时, 不同厚度相同厚度变化量时薄膜热导率之间的差异.

表2　相同厚度差不同厚度之间热导率之差

T able 2　Difference

of the therm al conductivity with various

thickness at the same thickness variation

不同厚度之间Πnm

热导率之差ΠJ ・m -1・K -1・s -1

400～30016

300～20029

200～10062

j

f (ρi ) , Φij , 且有

6

V (r ) =

k =1

∑a

2

k

(r k -r ) H (r k -r ) ,

3

3

(7) (8)

Φ(r ) =

k =1

∑

A k (R k -r ) H (r k -r ) ,

其中r 为距离量, H 为单位阶跃函数

H (x ) =

0, 　x

x ≥0.

(9)

r k 与R k 的选择满足r 1>r 2>, …, >r 6和R 1>R 2. r 1和R 1分别代表V 和Φ的截断半径; a 1, …, a 5, A 1, A 2由面心立方晶体平衡时的晶格长度、弹性模

量、内聚能(C ohesive Energy ) 等参量决定; a 6是中心势的短程势部分

[9]

, 其值见表1

[9]

.

表1　参数a k , r k , A k , R k 的值(nm)

(其中r k 和R k 是以晶格长度a =013615nm 为单位)

T able 1　V alues of p arameter a k , r k , A k , R k

参数

a 1

值

29. 059241

参数

R 1R 2A 1A 2

值

1. 22474491. 0000009. 80669416. 774638

参数

r r 2r 3r 4r 5r 6

值

1. 22474491. 15470541. 11800651. 0000000. 86602540. 7071068

a 2-140. 0568a 3a 4a 5a 6

130. 0731-17. 4813531. 8254671. 58749

图1　不同厚度下薄膜热导率随温度的变化

Fig 11　Variation of the thermal conductivity with tem perature for different thickness

　　在本文中, 粒子初始构形选择面心立方结构, 模拟采用的时间单位是110×10性边界条件.

对于本文所要计算的热导率, 根据G reen 2K ubo 法

[6～8]

-5

s , 模拟中引入周期

, 其计算公式为

λ=

〈J (t ) ×J (0) 〉d τ. 2

3k B T 0

∞

(10)

　　图2表示温度分别为100K 和600K 时, 热导率与薄膜厚度之间的关系. 可以看出, 在同一温度下热导率随膜厚非线形地变化(随着厚度的减小, 热导率急剧地下降) . 这是因为随着厚度的减小, 粒子受到薄膜表面边界散射的概率增大. 图2还表明, 在同一

　第4期

叶杰成等:分子动力学模拟铜薄膜的热导率323

电子参与导热. 实际上, 当膜很薄时, 晶界散射(grain

scattering ) 和缺陷对热导率的影响, 以及声子对热导率的贡献都不可忽视, 必须采用更加精细的物理模型. 对热流的计算既要考虑分子的动能, 也要考虑分子间的势能和相互作用, 并且要考虑分子团簇和分子聚集的影响, 这方面还有许多工作要做.

3. 100nm 、200nm 、300nm 和400

图2　Fig 12　Variation of with at 160～600K 范围内的热导率. 其

结果可用于分析微电子器件中铜薄膜的微尺度热传导问题.

[参　考　文　献]

[1]　T ien C L ,Chen G. Challenges in M icroscale C onductive and

Radiative Heat T rans fer [J].AS ME Journal of Heat T rans 2fer ,1994, 116::799-807.

[2]　Abrams on A R , T ien C L. Recent Developments in M i 2

croscale Therm ophysical Engineering [J].M icroscale Ther 2m ophysical Engineering ,1999, 3:229-244.

[3]　T enebaum A ,Ciccotti G, G allico R. S tationary N onequilibri 2

um S tates by M olecular Dynamics [J].Physical Review A , 1982, 25:2778-2787.

[4]　M ountain R D ,MacD onald R A. Thermal C onductivity of

Crystals:A M olecular 2Dynamics S tudy of Heat Flow in a T w o 2dimensional Crystal [J].Physical Review B ,1983, 28:3022-3025.

[5]　M ichalski J. Thermal C onductivity of Am orphous S olids

above the Plateau :Molecular 2Dynamics S tudy [J].Physical Review B ,1992, 45:7054-7065.

[6]　Castai G,F ord J ,Vivaldi F ,Visscher W M. One 2Dimensional

Classical Many 2Body System Having a N ormal Thermal C on 2ductivity [J ].Physics Review Letter , 1984, 52:1861-1864.

[7]　P oetzsch R H ,B ettger H. Interplay of Dis order and Anhar 2

m onicity in Heat C onduction :M olecular 2dynamics S tudy [J].Physical Review B ,1994, 50:15757-15763.

温度, , 但是

增加的速度随着薄膜厚度的增大而减慢, 逐渐趋于一稳定值, 即其宏观体形材料的热导率. 从图2中还可以看出600K 的曲线较100K 的平坦些, 这说明在温度较低时薄膜的热导率尺寸效应表现得更明显, 因此在低温的条件下, 对比较厚的样品容易观察到热导率的尺寸效应, 这是实验通常在低温下进行的原因.

图3给出了Nath 和Chopra 对铜薄膜的实验数[10]

据和本文结果的对比. 图中也给出了体形铜热导率值作为参考. 对于体形材料和厚度为300nm

的薄

图3　模拟结果与Nath 和Chopra

实验结果的比较

Fig 13　Results com pared with Nath &

Chopra experimental data

[8]　Allen M P , T ildesley D J. C om puter S imulation of Liquids

[M].Ox ford :ClarendonPress ,1989,230.

[9]　Ackland GJ ,T ichy G,Vitek V ,Finnis M W. S im ple N 2body

P otentials for N oble Metals and Nikel [J ].Philos ophical Magazine A ,1987, 56:735-756.

[10]　Nath P ,Chopra K L. Thermal C onductivity of C opper Films

[J].Thin S olid Films ,1974, 20:53-62.

膜材料, 本文的数据和Nath 和Chopra 的实验结果相

当吻合. 但对于厚度为10nm 的薄膜, 二者的数据有较大的出入. 这是由于在模拟中所采用的物理模型只考虑了几何边界散射对热导率的影响, 并且只有

324计　　算　　物　　理

第19卷　

MOL ECU LAR DY NAMICS SIMU LATION TO THE THERMAL

CON D UCTIVITY OF COPPER THIN 2FI LM

YE Jie 2cheng , 　T ANG Zhen 2an

(Department o f Electronic Engineering , Dalian Univer sity , , )

[Abstract ]　Thermal conductivity of copper films is simulated S ize and tem perature dependent effects of the films with 100～400nm thickness at 100～600K tem are to the simulation. [K ey w ords ]　thin 2film ; thermal [R eceived d ate ]2000-08-07; [R evised d ate ]2001-02-19

第19卷第4期2002年7月

计　　算　　物　　理

　CHI NESE JOURNA L OF C OMPUT ATI ONA L PHY SICS 　

V ol. 19,N o. 4　

[文章编号]　10012246X (2002) 0420321204

分子动力学模拟铜薄膜的热导率

叶杰成, 　唐祯安

(大连理工大学电子工程系, 辽宁大连　116023)

[摘　要]　采用分子动力学(M D ) 方法模拟铜薄膜的热导率, ～100～600K 范围内

铜薄膜热导率对尺寸及温度的依赖关系.

[关键词]　薄膜; 热导率; 分子动力学[中图分类号]　TK 123; O41113

[0　引言

分子动力学计算机模拟是研究微尺度热传输特

[1,2]

性的有力工具. 分子动力学是一种决定论(deter 2ministic ) 的原子级数值实验方法, 它不仅可以得到系

度梯度, 以至于无法准确确定系统的温度, 从而不能

正确地得到热导率和温度的关系;2) 由于系统不是均质的, 就不便采用周期性边界条件, 这就造成所模拟的系统粒子数目太多因而计算量巨大. 而采用作为分子动力学模拟基础的E MD 方法, 只需计算出平衡态下的热流函数J , 然后利用G reen 2K ubo 公式和统计物理规律即可求得热导率

[6,7]

统中每个粒子的运动轨迹, 还可以提供系统的时间演化过程信息, 进而可以通过统计力学得知系统的静态(热力学) 和动态(传输) 特性. 对一些理论上难以说明或实验中难以观察的现象给出基于物质微观结构及分子动力学关系的解释, 特别适合研究薄膜材料的微尺度热传导特性.

在微尺度热传输领域中, 由于尺度效应, 材料的热物性参数与宏观体形材料相比有很大的不同, 并存在着很多的未知量表明

[1]

[1]

.

本文采用E MD 方法对厚度在100～400nm 、温度在100～600K 范围内的铜薄膜热导率作了分子动力学模拟.

1　分子动力学模拟的基本原理和方法

分子动力学模拟的基本原理是建立一个合适的粒子系统, 利用牛顿运动方程确定系统中粒子的运动, 通过求得粒子动力学方程组的数值解, 决定系统中各个粒子在相空间中的运动规律, 然后按照统计物理和热力学原理得出系统相应的宏观物理特性.

对于系统温度保持不变的正则系综(NVT ) , 描述该系统的运动方程如下

r i =p i Πm i ,

p i =f i -ξ(r i , p i ) ×p i ,

[8]

, 热导率就是其中之一. 研究

, 薄膜材料的热导率是薄膜厚度的函数. 但实

[1,2]

验测定薄膜材料热导率, 不仅对设备要求很高, 技术上也有相当大的难度

. 因此分子动力学模拟在确

定微尺度热物性参数方面发挥着重要的作用. 本文通过分子动力学方法, 模拟计算铜薄膜的热导率. 从统计物理的角度可以将分子动力学分为平衡

态(E ) 和非平衡态(NE ) 两种模拟方法. 在以往的一些有关分子动力学模拟方法文献中, 都是用非平衡

[3～5]

态分子动力学(NE MD ) 方法来计算热导率. NE MD 是通过对系统施加扰动, 根据扰动所产生的

(1) (2) (3)

r i 、f i 、p i 、U (r 1, r 2, …, r n ) 和m i 分别是第i 个粒子

响应来求解有关的特性. 该方法主要用于处理远离平衡态或非均质材料, 但存在如下缺点:1) 为了得到足够大的温度梯度, 往往需要对系统施加很大的温

[收稿日期]2000-08-07; [修回日期]2001-02-19[基金项目]国家自然科学基金(5999555025) 资助项目

的坐标、所受到的合力、动量、其它粒子对其的作用势和质量, 其中ξ是约束温度不变的参数. 由于温度不变, 系统动能E k 必须满足 E k =0, 所以参数ξ是

[作者简介]叶杰成(1975-) , 男, 浙江永嘉, 硕士生, 从事微尺度热传输数值模拟方面的研究.

Δ

f i =-

i

U (r 1, r 2, …, r n ) .

322计　　算　　物　　理

第19卷　

ξ=

p ∑

∑|

i

i

×f i

p i |

2

其中V 为系统的体积:k B 为波耳兹曼常数, T 为系

.

(4)

统的温度, J 为热流函数, 可通过E MD 方法确定. 在

[5]

笛卡尔坐标下, 热流函数的微观表达式如下

J

α

在分子动力学模拟中, 关于材料系统内的粒子间相互作用势函数U (r 1, r 2, …, r n ) 的选取非常重要, 它决定了模拟结果正确与否. 本文采用的势函数是嵌入原子法(Embedded Atom Method , E AM ) 多体势函数

[9]

=

V

∑

j

E j v j -4

α

∑∑

j

r jk r jk (v j +v j ) .

r jk 9r k (≠j )

α

β

αβ

(11)

α、β=x , y , z , E j 为粒子, v j 为粒子j 的速度, 为粒子j , V 为系统的体积. E AM 是建立在H onhenberg 2K ohn 理论基础上

的, 该理论采用局域电子密度函数描述相互作用的原子中的一个原子对系统能量的贡献

U i =-f (ρi ) +

2

V ij

((6)

图1是厚度从100～400nm 的几种铜薄膜在不

同温度下的热导率. 从图中可以看出在某一温度下, 热导率随着膜厚的减小而减小. 热导率不是随着膜厚均匀地减小, 而是膜越薄, 热导率递减得越厉害. 由图1还可见, 在同一厚度, 热导率随温度增加而减少, 其原因是高温时晶格对声子的散射加强. 表2说明在温度为300K 时, 不同厚度相同厚度变化量时薄膜热导率之间的差异.

表2　相同厚度差不同厚度之间热导率之差

T able 2　Difference

of the therm al conductivity with various

thickness at the same thickness variation

不同厚度之间Πnm

热导率之差ΠJ ・m -1・K -1・s -1

400～30016

300～20029

200～10062

j

f (ρi ) , Φij , 且有

6

V (r ) =

k =1

∑a

2

k

(r k -r ) H (r k -r ) ,

3

3

(7) (8)

Φ(r ) =

k =1

∑

A k (R k -r ) H (r k -r ) ,

其中r 为距离量, H 为单位阶跃函数

H (x ) =

0, 　x

x ≥0.

(9)

r k 与R k 的选择满足r 1>r 2>, …, >r 6和R 1>R 2. r 1和R 1分别代表V 和Φ的截断半径; a 1, …, a 5, A 1, A 2由面心立方晶体平衡时的晶格长度、弹性模

量、内聚能(C ohesive Energy ) 等参量决定; a 6是中心势的短程势部分

[9]

, 其值见表1

[9]

.

表1　参数a k , r k , A k , R k 的值(nm)

(其中r k 和R k 是以晶格长度a =013615nm 为单位)

T able 1　V alues of p arameter a k , r k , A k , R k

参数

a 1

值

29. 059241

参数

R 1R 2A 1A 2

值

1. 22474491. 0000009. 80669416. 774638

参数

r r 2r 3r 4r 5r 6

值

1. 22474491. 15470541. 11800651. 0000000. 86602540. 7071068

a 2-140. 0568a 3a 4a 5a 6

130. 0731-17. 4813531. 8254671. 58749

图1　不同厚度下薄膜热导率随温度的变化

Fig 11　Variation of the thermal conductivity with tem perature for different thickness

　　在本文中, 粒子初始构形选择面心立方结构, 模拟采用的时间单位是110×10性边界条件.

对于本文所要计算的热导率, 根据G reen 2K ubo 法

[6～8]

-5

s , 模拟中引入周期

, 其计算公式为

λ=

〈J (t ) ×J (0) 〉d τ. 2

3k B T 0

∞

(10)

　　图2表示温度分别为100K 和600K 时, 热导率与薄膜厚度之间的关系. 可以看出, 在同一温度下热导率随膜厚非线形地变化(随着厚度的减小, 热导率急剧地下降) . 这是因为随着厚度的减小, 粒子受到薄膜表面边界散射的概率增大. 图2还表明, 在同一

　第4期

叶杰成等:分子动力学模拟铜薄膜的热导率323

电子参与导热. 实际上, 当膜很薄时, 晶界散射(grain

scattering ) 和缺陷对热导率的影响, 以及声子对热导率的贡献都不可忽视, 必须采用更加精细的物理模型. 对热流的计算既要考虑分子的动能, 也要考虑分子间的势能和相互作用, 并且要考虑分子团簇和分子聚集的影响, 这方面还有许多工作要做.

3. 100nm 、200nm 、300nm 和400

图2　Fig 12　Variation of with at 160～600K 范围内的热导率. 其

结果可用于分析微电子器件中铜薄膜的微尺度热传导问题.

[参　考　文　献]

[1]　T ien C L ,Chen G. Challenges in M icroscale C onductive and

Radiative Heat T rans fer [J].AS ME Journal of Heat T rans 2fer ,1994, 116::799-807.

[2]　Abrams on A R , T ien C L. Recent Developments in M i 2

croscale Therm ophysical Engineering [J].M icroscale Ther 2m ophysical Engineering ,1999, 3:229-244.

[3]　T enebaum A ,Ciccotti G, G allico R. S tationary N onequilibri 2

um S tates by M olecular Dynamics [J].Physical Review A , 1982, 25:2778-2787.

[4]　M ountain R D ,MacD onald R A. Thermal C onductivity of

Crystals:A M olecular 2Dynamics S tudy of Heat Flow in a T w o 2dimensional Crystal [J].Physical Review B ,1983, 28:3022-3025.

[5]　M ichalski J. Thermal C onductivity of Am orphous S olids

above the Plateau :Molecular 2Dynamics S tudy [J].Physical Review B ,1992, 45:7054-7065.

[6]　Castai G,F ord J ,Vivaldi F ,Visscher W M. One 2Dimensional

Classical Many 2Body System Having a N ormal Thermal C on 2ductivity [J ].Physics Review Letter , 1984, 52:1861-1864.

[7]　P oetzsch R H ,B ettger H. Interplay of Dis order and Anhar 2

m onicity in Heat C onduction :M olecular 2dynamics S tudy [J].Physical Review B ,1994, 50:15757-15763.

温度, , 但是

增加的速度随着薄膜厚度的增大而减慢, 逐渐趋于一稳定值, 即其宏观体形材料的热导率. 从图2中还可以看出600K 的曲线较100K 的平坦些, 这说明在温度较低时薄膜的热导率尺寸效应表现得更明显, 因此在低温的条件下, 对比较厚的样品容易观察到热导率的尺寸效应, 这是实验通常在低温下进行的原因.

图3给出了Nath 和Chopra 对铜薄膜的实验数[10]

据和本文结果的对比. 图中也给出了体形铜热导率值作为参考. 对于体形材料和厚度为300nm

的薄

图3　模拟结果与Nath 和Chopra

实验结果的比较

Fig 13　Results com pared with Nath &

Chopra experimental data

[8]　Allen M P , T ildesley D J. C om puter S imulation of Liquids

[M].Ox ford :ClarendonPress ,1989,230.

[9]　Ackland GJ ,T ichy G,Vitek V ,Finnis M W. S im ple N 2body

P otentials for N oble Metals and Nikel [J ].Philos ophical Magazine A ,1987, 56:735-756.

[10]　Nath P ,Chopra K L. Thermal C onductivity of C opper Films

[J].Thin S olid Films ,1974, 20:53-62.

膜材料, 本文的数据和Nath 和Chopra 的实验结果相

当吻合. 但对于厚度为10nm 的薄膜, 二者的数据有较大的出入. 这是由于在模拟中所采用的物理模型只考虑了几何边界散射对热导率的影响, 并且只有

324计　　算　　物　　理

第19卷　

MOL ECU LAR DY NAMICS SIMU LATION TO THE THERMAL

CON D UCTIVITY OF COPPER THIN 2FI LM

YE Jie 2cheng , 　T ANG Zhen 2an

(Department o f Electronic Engineering , Dalian Univer sity , , )

[Abstract ]　Thermal conductivity of copper films is simulated S ize and tem perature dependent effects of the films with 100～400nm thickness at 100～600K tem are to the simulation. [K ey w ords ]　thin 2film ; thermal [R eceived d ate ]2000-08-07; [R evised d ate ]2001-02-19