（原子与分子物理专业论文）氦掺杂诱发金属镍脆化的原子水平研究.pdf

摘要 在核技术和核工业领域中，氦的产生和聚集对于聚变堆和裂变堆第一壁金属 材料造成的损伤一直是困扰核材料使用安全的重大问题。氦对材料的影响主要体 现在氦的聚合、氦泡形成、氦与材料内缺陷相互作用，而导致材料的肿胀和脆化 断裂等一些宏观现象的发生，恶化了材料的力学性能因此很有必要对氦在金属 体内的聚集、迁移、扩散以及与各类缺陷的相互作用进行深入的研究，进而探讨 这些微观结构演化与材料力学性能之间的密切关系。目前，人们利用理论和实验 的手段试图从各种层次进行相关的研究。其中计算机模拟方法( 包括第一原理计 算，分子动力学方法和蒙特卡罗方法) 可以洞察氦原子在材料体内的聚集迁移行 为。然而，氦与缺陷的相互作用影响到材料力学性能的改变，有关方面的研究仍 知之甚少。本文以f e e 镍金属为例，应用g r i f f i t h 理论计算裂纹扩展过程体系应 变能的改变，在原子水平研究金属材料的脆化。首先设计了一个有微裂纹的镍金 属衬底，并用分子动力学方法计算裂纹扩展时体系在弹性应变中释放的能量。在 同样的缺陷条件下，将有氦团簇时的应变能释放值( e r ) 与无氦时的相应值( e 帕) 进行比较，通过计算氦团簇注入后应变能释放值的改变，研究它对材料脆化的影 响。 我们发现，裂纹扩展时系统应变能释放值e r 显著高于无氦团簇掺杂时的应 变能释放，e 帕。由于氦掺杂致使e r 的增加将驱动裂纹更易扩展，促使了材料的 脆化。该应变能释放值e 随着氦团簇与裂纹尖端距离( d t c ) 的增大，而非线 性地减小，直至d t c 达到临界值6 5 a o ( a o 为6 ：c 镍的晶格常数) ，e r 接近e r o 这表明d t c 在该临界范围以外时氦对材料脆化的影响消失。进一步的研究发现 氦掺杂导致e r 的增大主要来自于h e - h i 间近距离的强相互作用。这种相互作用 与氦团簇中氦原子与空位数比值有密切关系，随着团簇中氦空位比的增大，致 脆效应更加明显。综上所述，当金属材料中氦浓度或缺陷浓度提高并超过某一临 界值时，金属材料的脆化效应将增强 关键词：镍，氮团簇，裂纹缺陷，应变能释放，分子动力学模拟 a b s t r a c t i nt h ef i e l do fn u c l e a rt e c h n o l o g ya n dn u c l e a ri n d u s t r y , t h ep r o d u c t i o na n d a c c u m u l a t i o no fh e l i u ma t o m sl e a dt oah u g ei n c n a c et ot h es c c u l eu s eo fn u c l e a r m a t e r i a l s s u c ha st h es t r u c t u r a lm a t e r i a l so ff i s s i o no rf u s i o nf i r s tw a l l t h e a g g r e g a t i o no fh e ，t h u sf o r m i n gh eb u b b l e ，u s u a l l yd e g r a d e st h em e c h a n i c a l p r o r , e l r t i e so fm e t a l sa n de l i t e sf a t i g u el i f e t i m eo fm a t e r i a l s s oi ti sv e r yn e c e s s a r y t os t u d yt h ed i f f u s i o na n da g g r e g a t i o no fh ea t o m s ，a n dt h ei n t e r a c t i o n sb e t w e e nt h e h ea t o m sa n dd e f e c t si nm a l r i x , f u t t h e m l r m o r e ，t ob ei n s i g h ti n t ot h er e l a t i o n s h i p b e t w e e nt h em i e r o s l r u e t u r ee v o l u t i o no fh ea t o m sa n dm e c h a n i c a lp r o p e r t i e so f m a t e r i a l s a tp r e s e n t , al o to ft h e o r e t i c a la n de x p e r i m e n t a lw o r k sh a v eb e e nd o n ei n t h i sa 溉a m o n gt h e m c o m p u t e rs i m u l a t i o n s ( f o re x a m p l ea bi n i t i oc a l c u l a t i o n , m o l e c u l a rd y n a l m i c $ a n dm o n t ec a r l 们a a v a i l a b l et od e s c r i b et h eb e h a v i o ro fh e a g g r e g a t i o nv e r yw e l la ta t o m i s t i cs c a l e h o w e v e r , t h e r et i t l e s t i l lm a n yq u e s t i o n s u n s o l v e d f o re x a m p l e ，t h ei n t e r a c t i o no fh ed u s t e rw i t hd e f e c t s ，w h i c hi sc l o s e l y a s s o c i a t e dw i t ht h em e c h a n i c a lp r o p e r t i e so fm a t e r i a l s i nt h i sw o r k , t h eb r i t t l e - d u c t i l e o ft h em e t a l ，e s p e c i a l l yt h ef e en i c k e l ，w a ss t u d i e da ta t o m i cl e v e lb a s e do nt h e c r r i f l f i t l at h e o r y f i r s t , ac r a c k e ds u b s t r a t ew a sb u i l tb yi n t r o d u c i n gab l u n tc r a c ki n t o t h e f e en is u b s t r a t e t h e n , t h es l r a i n - e n e r g yo f t h ec r a c k e ds y s t e mr e l e a s e dd u r i n gt h e c r a c kp r o p a g a t i o nw a sc a l c u l a t e du s i n gm o l e c u l a rd y n a m i c sm e t h o d t h ee f f e c to f h e - d o p i n go nt h eb r i t t l e d u c t i l eo ft h es u b s l x a t ei s e v a l u a t e db yc o m p a r i n gt h e r e l e a s e ds t r a i n e n e r g yw i t hh e d o p i n g ( 助w i t ht h a tw i t h o u th e - d o p i n g ( e 而 w ef o u n dt h a tt h ev a l u eo fe ri sh i g h e rt h a nt h a to ff r o i ti n d i c a t e st h a tt h e h e - d u s t e ri n d u c e st h es u b s t r a t et ob eb r i t t l es i n c et h ei n c r e a s e ds l r a i n - e n e r g yp r o v i d e s ad r i v i n gf o r c ef o rt h ec r a c kg r o w t h w i t hi n c r e a s i n gt h ed i s t a n c eb e t w e e nt h ec r a c k 却a n d t h eh ed u s t e r ( d t c ) ，e rd e c r c a sn o n l i n e a r l yu n t i ld t cr e a c h e sac r i t i c a l v a l u e ，6 5 a o ( a o i st h el a t t i c l ：c o n s t a n to fn i ) b e y o n dt h ec r i t i c a l v a l u e ，e r a p p r o x i m a t e l ye q u a l st oe t o i ti m p l i e st h a tt h ee f f e c to f t h eh e - d o p i n g0 1 1t h ec r a c k n p r o p a g a t i o nc a nb en e g l e c t e d d e t a i l e da n a l y s i sr e v e a l st h a tt h ei n c r e a s eo fe rm i r d y c o m e sf r o mt h ei n t e r a c t i o nb e t w nh ea t o m sa n dn ia t o m ss u r r o u n d i n gt h eh e c l u s t e r f o rt h eh ec l u s t e r 诵t hl l i g hh et ov a c a n c yr a t i o ，t h es t r o n gh e - n ii n t e r a c t i o n i si nf a v o ro ft h ec r a c kg r o w t h f u r t h e r m o r e ，o n rs t u d ys h o w st h a tt h eb r i t t l e m e n to f t h ef i a c t u 他si nn iw i l lb ee n h a n c e dw i t hi n c r e a s i n gt h ed e n s i t yo f e i t h e rh ec l u s t e ro r f r a c t u l e s k e yw o r d s ：n i c k e lh e l i u m - c l u s t e r , c r a c k , s t r a i n - e n e r g y , m o l e c u l a rd y n a m i c s i r e u l a t i o n 论文独创性声明 本论文是我个人在导师指导下进行的研究工作及取得的研究成果。论文中除 了特别加以标注和致谢的地方外，不包含其他人或其它机构已经发表或撰写过的 研究成果。其他同志对本研究的启发和所做的贡献均已在论文中作了明确的声明 并表示了谢意。 作者签名 论文使用授权声明 日期：) 丑：! ：鲨 本人完全了解复旦大学有关保留、使用学位论文的规定，即：学校有权保留 送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内 容，可以采用影印、缩印或其它复制手段保存论文。保密的论文在解密后遵守此 规定。 作者签名：导师签名：陪寸谈 日期：三! z ! 鲨 第一章引言 第一章引言 1 1 氦对材料的影响和危害 在核能技术许多领域使用的金属材料中普遍存在着氦的产生及其引起的脆 化问题【1 3 】，如聚变堆第一壁结构材料的内表面受氘、氚等离子体的聚变产物口 粒子和中子的辐照，口粒子的注入及中子引发的( n ，口) 核反应，会产生氦原子的 积累：储存氚的金属容器中由于氚的衰变造成容器壁中高浓度氦( 3 h e ) i 拘积累， 降低了储氚材料的力学性能和使用性能：高通量散裂中子源中高能质子与靶、包 壳、窗材料的( p 口) 反应会在这些材料中产生高浓度的氦；核燃料钚的口衰变在 材料内部会引起氦积累。氦聚集成泡造成材料肿胀，这一现象在快中子增殖堆中 尤为明显。结构材料中氦泡引起肿胀的同时还造成高温氦脆卜7 】。所谓氦脆 ( h e l i u me m b r i t t l e m c n t ) 是指引入材料中的氦达到一定浓度时通常会在材料的基 体、位错、晶界、相界处形成纳米尺寸的析出相氮泡。氦泡的形成会导致材 料的延展率、蠕变断裂时间、疲劳寿命显著变小等等。它有多种表现形式，如高 温晶界型脆化、低温脆化 氦脆的研究是国防安全和能源安全的必然要求。例如反应堆的外围防线之一 的压力容器，通常是一巨大的焊接钢结构，这种压力容器会出现特殊问题，其中 最显著的就是中子辐照对脆性到延性的转变温度的影响；在航天技术上，人们仍 在进行氚和3 h e 的聚变研究，探求安全、高效的聚变能源。另外氦脆在科学研究 中也是需要解决的问题，例如北京的散射中子源( b s n s ) 中很重要的一方面研 究就是研制抗辐射损伤的靶体材料。 对于金属材料中的氦的研究大致可以分为两个方面。一是氦在金属中的微观 机制，提供金属材料的微观图像，并为解释宏观现象提出理论依据；另方面是 在不同条件下( 如氦的产率、浓度、不同材料的结构等) 研究含氦金属材料宏观 力学性能的变化，从而研制用于抵抗氦脆的材料，下面分别来介绍这两个方面。 1 2 研究历史和现状及进展 早在2 0 世纪6 0 年代，b a r n e s 8 首先报告了受辐照的钢由于氦泡在晶界聚集 而引起的高温氦脆。7 0 年代中期可控热核反应研究的进展推动了聚变堆材料的 研究。也促进了这领域的发展。在过去的几十年里人们对金属中氦的析出相( 氮 泡) 的形成机制做了大量的研究。研究基本上是从理论和实验上分别进行的。实 第一章引言 验上主要是利用离子加速器在材料中注入氦，或者利用氚衰变法引入氦。用透射 电镜、正电子湮没谱仪、中子小角散射、热解析谱、电子能量损失谱等多种方法 做微结构分析。理论方面现在主要是基于计算机模拟的理论分析，通常利用 s m o l u c h o v s k i 速率方程或基于原子间相互作用势的分子动力学方法、从头计算的 第一性原理方法、以及蒙特卡罗( m o n t ec a r l o ) 方法进行不同层次的微观机制的 研究。通过第一性原理和分子动力学方法可以计算氦在金材料中的各种能量关 系 1 2 i 关于氦的引入 常温常压下，氦在金属中的溶解度非常低，不可能通过通常的扩散方法在固 体中得到较高浓度的氦。金属中产生氦主要有下列途径： 一是某些元素受高能辐照引发( n 口) 反应生成氦并形成氦泡，这一现象 在快堆和聚变堆中尤为明显例如快中子的( n ，口) 反应【9 】有： ；m + ：厅7 一z a - 2 3 m + ；上k ( 1 2 1 ) z a 朋十o i 一， z a 一- - 2 4 肘+ ：n + ：日 ( 1 2 2 ) 其中材料m 的原子质量为a ，核电荷为z 。这种反应存在于所有核素中。另外 某些特性元素和热中子也会发生反应产生氦。对镍来说，反应式如下： 姐5 1 1 v l 。十0 1 疗曲 罢 厅+ ， ( 1 2 3 ) 嚣5 9 v i 十o i 玎曲 盏 阿+ ：胁 ( 1 2 4 ) ( n ，口) 反应产生的氦原子能量比较高，一般为几个m e v ，因而能够造成晶格 缺陷且带有，放射性，产生的氦浓度分布比较均匀。 二是储氚材料中有氚衰变产生3 h e 。反应式为： ；日_ ；胁+ 一l p + v 。 ( 1 2 5 ) 其中，一、h 分别是电子和中微子。粒子的能量不高，平均只有几个k e v ， 因而3 h e 的反冲能比较低，大概几个e v ，不会引起材料的离位损伤，这种引入 法方法得到氦浓度分布比较均匀。 三是采用了离子注入等方法向材料中引入4 h e 或3 h e 。这是一种比较直接快 捷的方法。离子注入法的优点在于各种实验参数可精确控制，并且实验简捷，费 用小，短时间内就可以得到高浓度的氦。 四是高压充氦技术。样品在高压氦气中充气，实际上由于氦的不溶解性质， 该方法很少应用。 2 第一章引言 1 2 2 氦在金属中的形成构型和相关能量 晶体中氦的行为本质上取决于氦原子在晶体内不同位置处的能量状态。由于 惰性原子的特殊性质，氦在金属中很难和金属原子发生反应。研究发现氦在金属 中的溶解度极其低，例如在我们研究的金属镍中氦原子在间隙位置的形成能很 高；氦间隙原子的迁移能很低，因而具有高的迁移性；氦和空位的结合能很高， 且氦原子与位错、晶界以及氦原子之间具有正的结合能，说明氦原子易被捕获于 空位、位错和晶界等缺陷处，也容易聚集而形成氦团簇。w i l s o n 等人在7 0 年代 利用计算机模拟进行理论计算【1 0 】，得到氦原子在金属中的存在行为和能量状态， 发现了氦的自捕获机制。由于间隙h e 形成能很高，而且h e 原子之间有正的结合 能，当氦原子数足够多时，就会自发在基体点阵中挤出空位，并形成相应的自间 隙原子，促成氦的自馅。自馅的过程实际上也是氦泡形核的先兆。 k m o r i s h i t a 1 1 1 等人计算了氦原子、空位及自间隙与氦空位团簇的结合能 研究表明这些结合能取决于氦空位团簇中氦与空位的比值。比值越大，空位与 团簇的结合能越大，说明氦可以抑制空位的热激发使团簇更加稳定。此外还利用 m d 方法得到能量参数，采用动力学蒙特卡罗方法( k i n e t i cm o n t ec a r l o ) 模拟 团簇的长期行为，进一步证实了氦原子能够抑制空位激发和促进自间隙原子的产 生来使氦空位团簇稳定存在 1 2 】。 1 2 3 氦的扩散及氦泡的成核生长过程 氦原子的扩散是氦泡形核和生长的基本要求。在可忽略辐照损伤下的扩散 【1 3 1 4 】，氦原子可以通过基体晶格进行扩散和迁移。具体可以分为间隙迁移，空 位迁移，受阻的间隙迁移。如下图所示： 易乇d o 舭o oooo 【 i ooo o o 。司0 5 o ooo 1 9 吕o 旬o 一。、o o o 。 o o o o o 。 l c ) c a ) 间隙迁移( b ) 空位迁移( c ) 受阻迁移 间隙迁移( a ) 是指氦原子通过间隙扩散，这种迁移所需要能量很低，所以在 很低的温度下就可以进行：空位迁移b ) 是指替位氦原子与近邻的空位交换位置 3 第一章引言 而发生迁移；受阻的间隙迁移( c ) 也叫离解机制，指氦原子脱离替位位置到间 隙位，然后以间隙扩散方式迁移直到被另一个空位捕获。这种激活能一般比较高 【1 5 】。 在辐照损伤下【1 3 】氦原子的扩散方式有离位机制和替代机制。离位机制是指 由于氦原子直接离位而产生的扩散；替代机制是指辐照产生自间隙原子运动到氦 原子所占据的空位，与空位复合，并将氦原子发射到间隙位，促进氦的扩散。这 中情况一般发生在温度比较高的时候。另外如果存在缺陷的情况下，譬如位错， 则可形成氦原子扩散的快速通道，称之为“管道扩散”。 前文所述由于氦原子的极其难溶解性，而且容易聚集成团。实际中氦的聚集 往往是外力使然，譬如离子注入、辐照、衰变等无不伴随着能量的释放。所以， 氦在金属中的聚集往往伴随着某些程度的离位损伤，但这主要取决于产生氦的速 率与原子离位速率的比值。在实际情况下，一般影响形核的因素有温度、氦的产 生速率、原子离位率等。试验上通过控制这些参数，我们可以分析氦的形核特点。 一般形核的研究分为两种不同情形 1 6 - 1 8 1 。第一种情形是低温( t o 5 t m ) 下发生。具体请参见文献 2 2 ，2 3 1 。 虽然进行了很多研究，但由于该领域的敏感性及研究手段的局限性，国内乃 至国际上对于金属中氦行为的研究仍然很不完善，而且能够有效克服上述的各种 4 第一章引言 氦脆问题的方法也比较有限。随着计算机的日新月异和新的研究工具不断涌现， 该领域各种尺度的研究工作仍然在进行。例如，最近基于第一性原理的研究仍在 计算氦在金属中的扩散及相关能量【2 4 】 1 3 含氦金属材料力学性能的研究进展 氦对材料宏观性能，特别是力学性能的影响主要表现为两类：一类是与材料 辐照损伤的联合作用，即加剧或减轻材料的肿胀( s w e l l i n g ) ，导致尺寸及形状的不 稳定性等：另一类是氦脆问题。核工业中普遍存在这些问题，譬如聚变堆中堆芯 部分，包括包层、第一壁、真空室和屏蔽等都存在各种力学问题。其中包层中管 道的固体力学问题是，由于液态金属横穿磁场时的电磁阻力增大，造成沿管道方 向的压力降很大，使管壁受到的内压大，从而存在强度问题。 我们这里主要关注氦对金属材料的力学性能的影响。氦对于材料力学性能的 影响与氦在金属体内的微观机制以及环境条件密切相关。实验上通过对温度、氦 的注入浓度及注入时间等参数的控制来研究材料力学性能的变化。研究氦脆的材 料大多集中在不锈钢上面，这是因为很多核材料中都要用到不锈钢材料。从温度 上讲，氦脆分为高温氦脆和低温氦脆。 有关h e 对金属力学性能影响的理论模拟研究，采用断裂力学方法。 s h i r a i 蚵【2 5 】用二维连续介质、弹塑性有限元法，对h e 的晶界脆性进行模拟。该 方法基于 c o nm i s e s 屈服应力标准，j 2 流变理论以及n 阶指数型硬化定律，对晶 界中氦泡产生不同向的应力应变曲线进行数值模拟，得出晶界氦泡引起的局部 缩颈和塑性损失等氦脆现象以及存在一个临界泡尺寸的结论，即在临界尺寸之 下，是不会导致晶界氦脆的。 1 3 1 高温氦脆 高温氦脆是由于应力增强的氦泡在晶界上长大所引起的，在t 芝o 4 5t t n 仃指 材料熔点温度) 尤为明显，这是堆材料氦脆的主要形式。研究高温氦脆，通常预 先注入或预先辐照，然后再进行力学试验，研究材料的拉伸、疲劳和蠕变等力学 性能 2 6 - z g 。在离子注入情况下，还经常采用双束注入方法，一束较重的例子如 n i + ，以造成辐照损伤，另一束为h c + 。 在高温区的拉伸试验表明：经块堆辐照产生氦原子后( 浓度约2 5 x 1 0 勺。不 锈钢的延展率下降1 0 1 5 ；而且改变拉伸试验的应变率，则发现断裂应变的 最大下降发生在最低应变率，即蠕变条件下。预注入不同剂量的蠕变试验表明， 只有氦原子浓度增加到某一临界值时，断裂时间和断裂应变才会发生骤然下降： 5 第一章引言 断裂方式也由晶内断裂模式向晶间断裂模式变化；随着测试温度和应力的增加， 断裂时问和断裂应变减小，蠕变时间也急剧减少。有人提出了氦泡呼吸模型：疲 劳拉伸试验中，存在一个临界频率，当在低频率环载荷作用下，在拉伸阶段晶界 氦泡可因空位扩散而迅速长大成为空洞，而在压缩阶段不可能愈合。这些结论基 本上可以描述高温氦脆的特点。 1 3 2 低温氦脆 低温氦脆是在高温氦脆提出之后才提出来的，早期的研究者认为氦脆只在 高温试验或含氦材料经过高温退火后才表现出来。2 0 世纪8 0 年代以来许多报道 证实了低温氦脆的存在。美国有关氚处理和核武器研究单位，如s a v a n n a hr i v e r 实验室和s a n d i a 实验室等在这方面进行了系统的长期实验观察。相对于高温而 言，低温情形下通常取温度t 0 4t m 【2 9 ，3 0 ，但高于室温下注入。有的试验也 采用高压冲氚的方式，然后经过长期时效研究材料的力学性能。如不锈铜中的实 验，不锈钢的总延伸率、断裂应变以及断裂强度都变小，时问越长氦脆越严重； 随着离位剂量( d i s p l a c e m e n td o s e ) 增大，金属硬化的增量减小，在高于l d p a 下饱 和；在氦浓度很低的情况下，氦泡对金属脆化的贡献可忽略，只有高于l d p a 时 才显著；离位率( d p a ) 和氦浓度是影响低温氦脆的因素，所以低温氦脆取决于二 者或二者的结合。 这些方面结论可以总结如下：充氚的效应和充氢引起的氢脆相似，但长期时 效后随着愈来愈多的氚衰变成氦3 ，材料的屈服强度提高，延伸率和断裂韧性均 进一步下降，说明等量的氮与氢( 氚) 相比，其致脆作用更加严重1 3 1 1 。 根据氦在基体中的性质及上述的研究。人们提出了一些抗氦脆材料的方法： 在贮氚的应用中采用氚溶解度低，扩散系数小的抗氢不锈钢以降低氦的积累量； 在不锈钢中形成适度弥散的碳化物沉淀相( 例如啊c ) ，使氦泡形核密度高而不容 易长大；针对氦在金属中的不溶性，在基体中添加易于成为氮原子的强捕馅中心 的元素，以阻碍氦的晶内扩散：还有强化晶界。 总之，在与核工业及军事应用有关的材料中，氦的产生和积累会引发许多问 题。氦虽然不溶于基体材料，但很容易被材料中的缺陷捕获，形成氦空位复合 物，在较高浓度或较高温度下形成微观氦泡导致材料宏观性能退化。通过多年 的工作，人们对氦的行为及其影响虽然有一定的认识，但研究结果仍然比较零散， 许多问题并未完全解决。其中氦脆问题( 特别是低温氦脆的机制) 虽被普遍关心， 但其规律性还远未搞清楚。因此有必要继续深入开展材料中的氦的行为，包括从 电子、原子层次的微观结构到宏观效应的研究，以及改进和发展针对各种使用情 况的模拟方法。 6 第一章引言 1 4 本文主要研究内容 上文对国内外的研究历史和现状的综述表明，对于金属中氮行为研究仍处于 初始阶段，尚无有效的方法克服上述的各种氦脆问题。其中，氦对金属材料力学 性能影响的研究甚少，但是这又是一个亟待解决的现实问题，该领域还有很多方 面需要解决和完善。本研究工作基于发展的原子势的基础上，利用分子动力学方 法研究氦团簇与金属内缺陷相互作用的过程，观察了氦团簇对含有缺陷的金属的 脆化行为的影响，为实验上改善储氚材料和核材料的力学性能提供理论依据。 文章的内容具体安排如下：首先回顾本工作的研究背景：然后在第二章中详 细介绍了分子动力学模拟原理和原子间的相互作用势：第三章中从原子尺度上研 究了氦对金属的脆化效应，并提出了合理的物理解释，最后在第四章中对整个论 文进行了总结。 7 第二章理论基础和方法 2 1引言 第二章理论基础和方法 本世纪以来，人们逐渐认识到，物质的各种物理性质都是与它们的微观结构 密切相关的。由于真实体系中粒子间相互作用十分复杂，使传统的分析方法对大 多数问题都不能直接获得解析解。计算机的诞生为人们对复杂问题的研究提供了 新的方法和工具。人们针对具体问题设计一个合理的模型，由计算机完成大部分 的计算工作，从而获得复杂体系的微观信息，并从中分析得到体系的宏观物理性 质与行为，这就是计算机模拟方法的基本思路 新近发展的计算机模拟方法介于传统的实验和理论分析这两种方法之间，起 着双重角色的作用。图2 1 1 给出了这三者之间的关系。 一方面，计算机模拟和解析理论都依赖于某个基本的理论模型。解析理论通 常要在这个基本模型的基础上再做出更多的假设和近似，然后才能进行理论推导 和预言。而计算机模拟则是对这个基本的理论模型直接进行模拟和数值计算，从 而获得这个基本模型的准确解。所以，只有把计算机模拟的结果与解析理论预言 的结果进行比较，就可以判断解析解近似的正确与否。从这个意义上说，计算机 模拟可以起到检验理论的作用。而且在计算机模拟中，可以包含一些复杂却又比 较符合实际情况的因素，从这一意义上说，计算机模拟具有比较解析理论方法更 加优越的特点。另一方面。计算机模拟的结果可以和实验相比较，从而有助于解 释新的实验现象和为实验中所遇到的疑难问题提供有价值的启示。另外，当物理 实验尚不能或难以进行时，计算机模拟实验可以代替部分实验，或优化实验参数， 从而得到实验上所不能得到的结果。当然由计算机模拟得到的结果最终应该接受 实验的验证。 最早的计算机模拟实验【3 2 】奠定了现在称之为m o n t e c a r l o ( 简称m c ) 模拟的 基础，它是一种随机方法。在m o n l c c a r l o 模拟方法产生后不久，又诞生了另一 种确定性的模拟方法，这就是分子动力学( m o l e c u l a rd y n a m i c s ，简称m d ) 模拟方 法。 1 9 5 7 年和1 9 5 9 年，a i d e r 和w a i n w r i g h t 3 3 ，3 4 通过对原子经典运动方程( 也 就是牛顿运动方程) 求数值解，首次实现分子动力学方法对硬球相互作用体系的 模拟。1 9 6 4 年，r a h m a n 3 5 把m d 方法应用到l c n n a r d j o n e s 相互作用体系中， 并在模拟中采用在每一微小时间步长内计算每个粒子所受到的作用力的方法，称 8 第二章理论基础和方法 图2 1 1 计算机模拟与实验和解析理论方法之间的关系 为s t e pb ys t e p 的程序步骤，以研究体系的时间演化。从而使得分子动力学方法 对实际体系的应用成为可能。现在的分子动力学模拟基本上都采用r a h m a n 的 s t e pb ys t e p 模拟步骤 早期的计算机模拟所依据的模型都是非常简单的，模拟的内容主要集中在一 些很简单的单原子流体或固体整体性质研究。近三十年来，随着微电子技术的迅 速发展，计算机的功能越来越强，速度也越来越快，这就为计算机模拟的发展提 供了非常有利的条件。现在计算机模拟所针对的体系越来越复杂，模拟应用的范 围也越来越广。固体表面结构、液体结构、离子注入、辐照损伤、异质催化、酶 的行为、薄膜生长等，都可以通过计算机模拟进行研究。 分子动力学模拟方法中，势能函数的选取是材料的计算机模拟的核心问题。 势能函数选取的好坏直接关系到模拟结果的正确与否。势能函数总是表示成粒子 9 第二章理论基础和方法 坐标( r ) 的数学函数形式，它给出体系总能量与粒子坐标的关系。如果这个函数对 体系豹描述越准确，那么用它来进行模拟计算所得到的物理体系的性质就与真实 情况越接近。 严格地说，材料中原子和电子是遵循量子力学的运动规律的，因此要精确的 描述原子与原子，原子与电子以及电子与电子之间的相互作用，应该采用第一性 原理的计算方法( a bl a i t i o ) 。然而，我们知道，材料中它们之间的相互作用是典 型的多体问题，不可能从量子力学上严格求解，而且即便在进行大量的假设和近 似之后，依然需要巨大的运算量，这远非一般的计算机所能胜任。因此，基于原 子间的经验和半经验的相互作用势的计算机模拟，便成为一种简便易行而且行之 有效方法，并且起到越来越不可替代的重要作用。 早期采用的原子问作用势主要是半经验的对势，例如l e n n a r d - j o n e s 势 3 6 】、 m o r s e 势 3 7 1 和b o r n - m a y e r 对势【3 8 】等，模型较为简单，对于很多情况并不能有 效的拟合。近年来，发展出来的多体势，能够精确描述更为复杂的系统。如 f i n n i s s i n c l a i r 多体势 3 9 ，4 0 】以及嵌入原子法( e m b e d d e da t o mm e t h o d ，e a m ) 4 1 】 等方法已能比较成功地描述金属与合金中原子的相互作用后来，p e t t i f e r 4 2 ， 4 3 又推广了f i n n i s - s i n c l a i r 势中的紧束缚近似下的平方根项，把三体和更高阶相 互作用考虑在内，以预言周期表中结构行为的趋势。b a s k e s 4 4 提出的修正嵌入 原子法( e m b e d d e d d e f e c t m c t h o d ，e d m ) 包含了与角度有关的原子相互作用。 2 2 分子动力学方法( m d ) 模拟 分子动力学( m o l e c u l a rd y n a m i c s ，简称m d ) 模拟的开创性工作是1 9 6 0 年 j b g i b s o n 等人关于低能辐照损伤的研究。虽然四十多年来，m d 模拟随着一代 又一代新计算机的出现而有了长足的发展，但其基本原理并未改变。m d 模拟已 经用于解决固体中原子碰撞中的众多问题，如缺陷的形成和迁移、离位阈能的计 算、移位级联的空间构型、级联发展过程、溅射机理、高能量密度级联以及薄膜 生长、表面重构、界面混合等计算机模拟的新领域。随着原子间相互作用势的发 展，m d 模拟已被广泛应用到对c 、s i 、g - a 、a s 等重要的半导体材料的模拟上。 它已成为材料科学研究领域里对微观结构和微观过程进行研究的一个重要方法。 固体中原子碰撞现象的m d 模拟程序一般都是在下述近似基础上进行设计 的： 体系完全可以由经典力学处理：即粒子的运动服从经典运动方程，不考 虑相对论效应和量子效应。在低能情况下，经典近似是成立的； 粒子间的相互作用可用经验势函数近似描述； 所有的势函数都和粒子的动能无关； 1 0 第二章理论基础和方法 势能函数的形式不因粒子在碰撞过程中可能发生激发或电离而改变。 2 2 1 分子动力学方法原理 我们知道体系的微观状态可以由组成体系的n 个粒子的坐标和动量完全确 定。在原子层次的模拟中，原子为组成体系的基本粒子，体系的h a m i l t o n i a n 量 为原子坐标和动量的函数。如果认为体系中粒子的运动服从经典运动规律，那么 体系的h a m i l t o n i a n 量就可以进一步表示为所有原子的动能k 和原子间相互作用 的总势能。之和，电子对原子运动的平均效应在势能m 中反映。如果用 芦“= ( i ，云，昂) ，伊”= ( 暖，呒，露) ，卢”= ( 五，融，芦。) ( 2 2 1 ) 分别表示体系中n 个粒子的坐标、速度和动量，那么，体系的h a m i l t o n i a n 量为： 日g ”，多”) = k ”) + 中( ”) k g ”) = 善筹 ( 2 2 2 ) ( 2 2 3 ) 体系的运动方程可以有多种表达形式。最基本的一种表达方式是拉格朗日运动方 程： 旦罢一鲁：0 ( 2 2 4 ) d ta 只谚 其中工( ，) 是拉格朗e t 函数，它和体系的h a m i l t o n i a n 量的关系是： 把( 2 2 5 ) 式代入( 2 2 4 ) 式，就可以得到体系的h a m i l t o n i a n 正则方程： i 生：盟：立 i 2 一= ： i d t 8 p ,m ， 【亟：罢：吨蛳一)d 【t两 、 ( 2 2 5 ) ( 2 2 6 ) 、， 一n r 一r 娥 一 一b一叶 = 、， _ v r ，l l 第二章理论基础和方法 v ，= 苦r + 毒歹+ 毒石由于警= z ，所以- v s o ( 一是体系中所有其他粒子 对i 粒子的总作用力。 m d 方法的核心是数值解正则运动方程( 2 2 6 ) ，求出体系的相空间轨道 ( f ) ，( f ) 或芦”( f ) ，帚”( ，) 。显然方程( 2 2 6 ) 能否求解，求得的解是否正确， 关键在于能否给出正确的粒子间相互作用势函数( f ”) 。 2 2 2 数值积分方法 m d 模拟采用的势函数都是连续势。在连续势作用下，粒子各个时刻都受到 力的作用，运动轨迹一般不再是直线。系统运动方程组( 2 2 6 ) 的数值积分可以 有多种方法。算法选择的关键在于差分方程能否稳定收敛，计算是否准确可靠， 在此基础上再看算法是否简单。分子动力学常用的算法有中心差分法( c e n t r a l d i f f e r e n c es c h e m e ) 4 5 、平均力法( a v e r a g ef o r c em e t h o d ) 4 6 、预报修正法 ( p r e d i c t o r - c o r r e c t o rs c h e m e ) 4 7 、v e r l e t 算法( t 1 1 ev e r l e ts c h e m e ) 4 8 等。我们在 计算中使用的是v e r l e t 算法，下面对这一算法作一介绍。 v e r i e r 算法是分子动力学中较简单但使用却最为普遍的一种方法，把粒子随 时间的位移作t a y l o r 展开： 机咖枷埘( f ) + a t 2 哟+ 击学nd ( 掣 ( 2 2 7 ) 私咖删埘( f ) + 一 掣n d ( 掣 ( 2 z 8 ) 将( 2 2 7 ) 、( 2 2 8 ) 两式相加，忽略陋) 及以上的高阶项，得 弓( f + r ) = z ( f ) 一弓o 一，) + ，2 j ； ( ，) ( 2 2 9 ) 将( 2 2 7 ) 、( 2 2 8 ) 两式相减，忽略陋) 及以上的高阶项，得 艿( f ) = 【弓o + ，) 一亏( f 一r ) 1 ( 2 f ) ( 2 2 1 0 ) 在恒定温度的m d 模拟中经常采用蛙跳算法( t h el e a pf r o ga l g o r i t h m ) 4 9 ， 它是与v e r l e t 算法等价的一种算法： 第二章理论基础和方法 地+ 址) = 毒( f 一 血) + j ；| ( f ) f( 2 2 1 1 ) k ( f + f ) = 亏( f ) + 毒( f + f 必。 根据( 2 2 1 0 ) ，t 时刻速度由下式求得： i ( f ) = i o + 圭址) + i o 一 址) ( 2 2 1 2 ) 也可以用下式的计算来提高精度： 毒( ，) = 上【毒( f + 址) + 毒o 一 厶f 为+ 古晴o 一址) 一毒( f + f ) 1 f ( 2 2 1 3 ) 2 2 3 边界条件 材料的宏观性质是大量粒子的统计行为。模拟体系的粒子必须非常庞大，才 能准确再现真实体系的行为另一方面，受计算机内存及容量的限制，m d 模拟 能够处理的粒子数是有限的。m d 模拟通过在无限大的体系中取出一块小晶体作 为模拟体系，模拟体系的边界原子要受到外界原子的作用力，在平衡构型中，此 力正好抵消晶体内部对边界原子的作用。而且边界原子在受到扰动时，它要向外 界传递能量，甚至离位。这就要求对模拟体系的边界提出合适的计算模型。一个 最简单的方法是增加体系中的原子的个数，随着计算机的飞速发展，这是一个趋 势。目前，根据具体的研究问题不同，模拟体系的边晁通常有以下几种处理方法： 周期边界条件 通常情况下，为了解决真实体系粒子数的庞大和计算机处理能力有限之间的 矛盾，模拟中一般采用周期性边界条件。 一般取一立方体为模拟体系，体系中含有几百至几万个粒子。如果不模拟表 面，就使中心原胞在三维空间上重复排列；如要模拟表面，就使模拟立方体在二 维水平面上重复排列。于是体系中的像粒子将在三维和二维空间周期性的重复出 现。则l 粒子的像坐标为： ( 石+ z 厶， y t + m l , ， z + ，z 厶) l , m , n = or 1 ，土2 ， 其中厶，厶，厶为立方中心原胞三边边长。体系的势能包括粒子与像粒子 之间的相互作用在内。 应用周期边界条件后，整个体系的赝粒子数无限多。当中心原胞中的一个粒 子移动时，其像粒子也会发生相应地移动。这样，若一个粒子离开了中心原胞， 必定有它的个像粒子从另一面进入中心原胞。因而，中心原胞中的粒子数总是 保持不变的。 第二章理论基础和方法 应用周期性的边界条件隐含着一个要求，即相距为工( l = m i n 6 k ，句，1 0 ) 的两粒子的空间相关性是可以忽略的，也就是说，粒子问的相关长度应远小于模 拟计算的立方体的边长。这就要求模拟体系的粒子数要足够大。 固定边界条件 对于正则系综0 叮v t ) 的分子动力学模拟，有时采用固定边界条件，比如对于 粒子与固体表面相互作用体系，在模拟中将靶底层或者下面几层原子固定，以模 拟衬底对模型靶的相互作用。 自由边界条件 如果模拟原胞包含的粒子对于所研究的问题足够大，或者为了减小周期性边 界条件或固定边界条件对所研究问题的影响( 冲击波的反射等) ，m d 模拟中也常 采用自由边界条件，即对边界原子不加任何限制。 2 2 4 时间步长( a t ) 的选取 我们从前面介绍的体系运动方程组的数值积分可以看出，从计算精度考虑， 时间步长t 的越小越好。但计算效率决定t 不能太小t 的选择要协调考虑这 两方面的要求。我们在m d 模拟中采用变步长的方法。该方法先设置一个比较 小的时间步长，在计算过程中每隔一定的时间步自动调整步长。调整的原则是该 时刻体系中运动最快的粒子在一个时间步内所走的路程不超过晶格长度a o 的2 ， 即a t = o 0 2 a o v 。，v 。为该最快粒子的速度。上面这个百分比还可调整，尽可能 达到在满足计算精度的要求( 如体系总能量误差不超过