（固体力学专业论文）纳米晶铁和钽在单向拉伸载荷作用下的分子动力学模拟.pdf

暑北工越走雪龋士雪他佟立 摘要 本文利用分子动力学方法研究了纳米晶b c c 金属铁和钽在单向拉伸载荷作用 下的力学性能和微观结构演化。结果表明，在孔隙率非常低的情况下，纳米晶铁 和钽的弹性模量与试样的相对密度成线性关系。两者的拉伸强度随着晶粒尺寸的 增加而增大，显示出反常的h a l l p e t c h 关系。对纳米晶钽的模拟结果表明纳米 晶钽的流动应力随着温度的升高而增大，并随着应变率的增加而增大。在变形过 程中，两者都出现了晶粒旋转，晶界滑移和运动以及位错的运动。低温模拟结果 表明，纳米晶铁和钽在变形过程中都出现了应力诱导的相交，其中f c c 结构含量 的最大值和抗拉强度成线性关系，并由此可以推出相变发生的临界应力值。应变 率虽然不是影响相变结构含量最大值的根本因素，但却可以影响到相变发生的起 始点。纳米晶钽在变形过程中出现了明显的晶间断裂现象，而铁则没有，在本文 的研究范围之内，裂纹形成所对应的试样变形和晶粒尺寸无关而与应变率和模 拟温度有关。纳米晶钽在变形过程中发生局部孪晶化。 纳米晶b c c 金属的变形机制可解释为晶粒旋转，晶界滑移和运动以及位错的 运动。这与f c c 金属相似，而当应力足够大时试样局部通过发生相变来协调试样 的变形。对于纳米晶钽晶间断裂也是一种协调试样变形的机制。这不同于f c c 金属。 关键词：分子动力学，纳米晶钽，纳米晶铁，晶间断裂，孪晶，相变 蕾站工t 土事龋士簟位话盘 i i i i i 1 1 1 _ m o l e c u l a r d y n a m i c ss i m u l a t i o n so fn a n o c r y s t a l l i n e t h f eu n d e ru n i a x i a lt e n s i o n a b s t r a c t u s i n gm o l e c u l a rd y n a m i c ss i m u l a t i o n ，w eh a v ei n v e s t i g a t e dt h e m e c h a n i c a l p r o p e d i e sa n dt h em i c r o s t r u c t u r a le v o l u t i o no fn a n o c r y s t e d l i n et a & f e ( g r a i ns i z e v a r i e sf r o m3 2 5n l nt o1 2 9 9n mf o rt aa n dv a r i e sf r o m2 8 2 n mt o8 4 7 n r n lu n d e r u n i a x i a lt e n s i l el o a d i n g t h er e s u l t ss h o wt h a lt h ee l a s t i cm o d u l u so fn a n o c r y s t a l l i n e f e & t a , w i t hn e a r l yz e r op o r o s i t y , h a sal i n e a rr e l a t i o n s h i pw i t ht h er e l a t i v ed e n s i t yo f t h es a m p l e s t h et e n s i l es t r e n g t ho fb o t hf ea n dt ai n c l e a s e sa st h eg r a i ns i z ei s i n c r e a s e df o rt h eg r a i ns i z er e g i o ns t u d i e d r e v e a l i n gar e v e r s eh a l l p e t c he f f e c t t h e s i m u l a t i o nr e s u l t so fn a n o c r y s t a l l i n et ai n d i e a t et h a tt h ey i e l ds t r e s so ft ai n c r e a s e s w i t hi n c r e a s i n gs t r a i nr a t e ( 2x1 07 2 1 0 9 s 一) o rd e c r e a s i n gt e m p e r a t u r e ( 1 0 4 1 5 0 0 k ) g r a i nr o t a t i o n g r a i nb o u n d a r ys l i d i n go rm o v e m e n ta n dd i s l o c a t i o n m i g r a t i o na r eo b s e r v e dd u r i n gt h ed e f o r m a t i o np r o c e s so ff e & t as a m p l e s t h er e s u l t s o fs i m u l a t i o n sp e r f o r m e da tn e a r l yz e r ot e m p e r a t u r ed e m o n s t r a t et h a ts t r e s s - i n d u c e d p h a s et r a n s i t i o n sf r o mb o d y c e n t e r e dc u b i c ( b c c ) t of a c e - c e n t e r e dc u b i c ( f c c ) a n d h e x a g o n a lc l o s e - p a c k e d ( h c p ) s t r u c t u r e st a k ep l a c ed u r i n gt h ed e f o r m a t i o np r o c e s s o fb o mf ea n dt a , a n dt h em a x i m u mf r a c t i o no ff c ca t o m sv a r i e sl i n e a r l yw i t ht h e t e n s i l es t x e n g t h f r o mw h i c hw ec a nc o n c l u d et 1 1 a tac r i t i c a ls t r e s se x i s t sf o rt h ep h a s e t r a n s i t i o nt ot a k ep l a c e s t r a i nr a t em a yn o ti n f l u e n c et h em a x i m u mf r a c t i o no f f c c h c pa t o m h o w e v e r , i tw i l li n f l u e n c et h es t r a i na tw h i c hp h a s et r a n s i t i o nb e g i n s t oo c c r r i n t e r - g r a n u l a rf r a c t u r ew a sc l e a r l yo b s e r v e di nt h ed e f o n - n a t i o np r o c e s so f n a n o c r y s t a l l i n et a ，w h e r e a ss i m i l a rp h e n o m e n o nd i dn o to c c u rr e m a r k a b l yd u r i n gt h e d e f o r m a t i o no ff e i nt h er a n g eo fs i m u l a t i o n sp e r f o r m e dh e r e ，t h es t r a i na tw h i c ht h e c r a c kf o r m a t i o ni si n i t i a t e di ss e e nt ob en o ta f f e c t e db yt h ea v e r a g eg r a i ns i z e b u t i n f l u e n c e db yt h es t r a i nr a t ea n dt e r n p e r a t u r e ，t w i n i n gi sa l s oc l e a r l yo b s e r v e di nt h e s i r e u l a t i o n so f t a s i m i l a rt on a n o c r y s t a u i n ef c cm e t a l s t h ed e f o r m a t i o nm e c h a n i s mo f n a n o c r y s t a l l i n eb c cm e t a l sc a nb ee x p l a i n e d 嬲b e i n gd o m i n a t e db yg r a i nr o t a t i o n g r a i nb o u n d a r ys l i d i n go rm o v e m e n ta n dd i s l o c a t i o nm i g r a t i o n n o n e t h e l e s s ，w h e nt h e f l o ws t r e s si sg r e a te n o u 曲，p h a s et r a n s i t i o nw i l lo c c u rt oa c c o m m o d a t et h e d e f o r m a t i o no fb o t hf ea n dt a ：i na d d i t i o n 。d u r i n gt h ed e f o r m a t i o np r o c e s so ft a , i n t e r - g r a n u l a rf r a c t u r ei sa l s oa ne f f e c t e dm e c h a n i s mt h r o v g ht h ef o r m a t i o no f w h i c ha l a r g ec o m p o n e n to ft h es t r a i ni sa c c o m m o d a t e d ，n i ct w ol a t e rm e c h a n i s m sa r eq u i t e d i f f e r e n tf r o mw h a th a v eb e e no b s e r v e di nt h es i r e u l a t i o n so fn a n o c r y s t a l l i n ef c c m e t a l s k e y w o r d s ：m o l e c u l a rd y n a m i c s ；n a n o e r y s t a l l i n ef e ；n a n o e r y s t a i l i n et a ； i n t e r - g r a n u l a rf r a c t u r e ；d e f o r m a t i o nt w i n n i n g ；p h a s e t r a n s i t i o n ， n i 地z 叠太亏赣士簟任话盘 第一章引言 自从1 9 8 1 年德国学者g l e i t e r 等提出并首次应用惰性气体凝聚结合原位冷压 成型法在实验室制备出纳米晶体样品以来，世界各国对这种新兴材料给予广泛的 关注，并相继开展了相关的研究工作。近十几年来随着纳米材料制各工艺和材 料微观结构分析方法等高尖端技术的快速发展，人们对这种材料的研究亦越来越 深入。随着物质某一方向的特征尺寸降低到纳米量级，纳米材料表面电子结构和 晶体结构发生变化，产生了宏观物体所不具有的四大效应：小尺寸效应、量子效 应、表面效应和界面效应。这使得纳米材料具有传统材料所不具有的系列优异 的力、电、磁、光学和化学等宏观特性，从而作为一种新型材料在宇航、电子、 冶金、化工、生物和医学领域展现出广阔的应用前景。纳米材料的研究成为当今 世界材料科学、凝聚态_ 物理、化学等领域中的一个热门课题，并被誉为“2 1 世 纪的新材料”而得到了世界各国的承认和深入研究。 通常，人们将晶粒尺寸小于1 0 0 纳米的金属定义为纳米盒属，晶粒尺寸在 1 0 0 一1 0 0 0 纳米范围内的金属定义为超细晶金属，晶粒尺寸大于1 0 0 0 纳米的则 为一般的会属。目前关于纳米材料的力学性能研究，已经涉及到弹性模量，硬度， 断裂韧性，压缩和拉伸的应力、应变行为，应变率敏感性，疲劳，蠕变和超塑性 等方面。纳米金属与普通的粗晶金属相比，具有许多不同的性质和优良的性能。 具体表现在：最初的研究表明，纳米金属的弹性模量比粗晶多晶材料大幅降低， 最大降幅高达7 0 。g l e i t e r 的解释是纳米金属随着晶粒尺寸减小，其晶界与三 叉晶界所占体积分数较大，这些晶界区域的平均原子间距较大，若其原子间相互 作用势与晶粒内相同，则其弹性模量比晶粒内小，故纳米材料的平均弹性模量比 多晶粗晶材料下降。但是随后的深入研究表明，纳米金属的弹性模量与其孔隙率 大小密切相关。由于孔隙或缺陷减少了外加应力的有效受力面积，具有裂纹形状 的孔隙还会引起应力集中，因此纳米金属中的孔隙，缺陷或裂纹使其弹性模量降 低。 纳米金属材料与普通粗晶，多晶材料的另一个显著不同的地方涉及h a l l p e t c h 关系。h a l l - - p e t c h 关系，亦即屈服应力与平均晶粒尺寸的平方根的倒数成 线性关系，曾是长期以来提高金属材料强度的指导规律之一。根据h a l l - - p e r c h 关系，对屈服应力有 置瞎工 太掌硪t 掌位静支 盯2 吼+ 脚一1 7 2 其中民为移动单个位错时所需的点阵摩擦力。七是常数，随不同的材料而不同。 对硬度应有 h 。= 凰+ d 。” 其中甄，k 。是常数。它们都是半经验公式，可用位错塞积或位错网来解释。纳 米金属材料制备成功以后，当晶粒尺寸小于1 0 0 纳米以下时，h a l l - - p c t c h 关系 是否成立已经引起了许多学者的兴趣。目前关于这方面的研究甚多。在纳米尺度 范围内，这一经典公式遇到了挑战。对纳米纯金属c u 的研究表明，随着晶粒的 减小，既有正h a l l - - p e t c h 关系，也有反h a l l - - p e t c h 关系。但对纳米晶p d 的实 验表明，其显微硬度保持基本不变。随着对纳米材料的进一步了解，发现制备过 程中的许多因素，比如密度，压实的压力，都会对纳米材料的内在关系产生很大 的影响。y a n g 【1 l 等人考虑密度及固结压力对纳米晶体c u 的显微硬度的影响时发 现：随着压实压力增大，纳米晶体c u 的密度升高，其显微硬度也升高，最高达 到粗晶铜的6 倍，因此实际上是密度影响了c u 的硬度。但随着晶粒尺寸的减小， 呈现出混合的h a l l - - p e r c h 效应，即随着晶粒尺寸的减小，硬度有时增加，有时 减小。目前有一种倾向认为：由于改变制备过程中的原子簇大小而合成的纳米金 属材料，随着晶粒减小，硬度上升：而对于一个通过多次回火而形成的不同晶粒 尺寸的样品而言，随着晶粒减小而硬度下降。s c h i o t z 2 l 运用分子动力学模拟纳晶 铜后发现：随着晶粒尺寸的减小，屈服应力降低，显示出软化趋势。 大量的文献指出，纳米晶金属由于晶粒尺寸大大低于普通多晶材料，导致其 力学性能。塑性变形机制亦有所不同。在糨晶金属中，塑性交形主要是由于晶粒 内规则晶格的线缺陷位错的运动产生的。位错可以贯穿整个晶粒，且两两之间 可以相互作用。晶界经常阻碍位错的传播，导致位错在晶界上堆积，从而导致材 料的硬化。但是当晶粒尺寸在纳米尺度范围内时，晶界原子的百分比大大增加， 传统的基于位错塞积理论的塑性变形需要重新认识1 3 】。当晶粒尺寸减小到2 0 纳 米时，在载荷作用下，纳米晶体由仅可弹性畸变的晶粒和黏性、可物质扩散的晶 界组成，晶界滑动和晶粒转动必然发生n 通过对纳米晶铜进行分子动力学模拟 发现，当平均晶粒直径小于1 5 r i m 时，材料的塑性变形机制从以位错运动为主过 2 叠非工业太掌硪t 雪幢话立 一i i i 一 渡到以晶界滑移为主【5 l ，导致纳米晶金属表现出反常的h a l l p e t c h 关系。纳米 晶铜的平均晶粒尺寸在1 0 1 5 n m 范围内时，流动应力出现最大值。 纳米力学的研究受到三项技术发展的的推动，其一来自于材料科学与工程， 尤其是对更强固，更轻、更韧、更具有多功能性的材料需求。纳米材料的出现使 得这一要求成为可能。一般情况下一种材料可以具有很好的韧性，或者很高的强 度，但很难两者兼而有之，通常是材料的强度越高，韧性就越差，如陶瓷；韧性 越好，强度一般就越低。当晶粒细化为纳米尺度后，由于塑性变形机制的不同， 则可能实现高强度高韧性的结合，如韧性陶瓷。实现高韧高强度可以通过在纳米 晶材料中保留一部分微晶，这样在变形过程中，纳米尺度的晶粒保证材料具有很 高的强度，而微晶可以稳定金属的塑性变形。另外还可以在纳米晶材料中使高角 晶界占优而有利于晶粒的旋转，这样纳米尺度的晶粒尺寸保证材料具有很高的强 度，而晶粒旋转协调材料的塑性变形【6 】。因此理解纳米金属的微观结构与其表现 出的宏观的力学性能之间的关系，从而从微观上揭示变形机理对于设计满足定 应用要求的材料具有根本性的意义。 分子动力学计算机模拟是研究复杂凝聚态系统的有力工具。这一技术既能得 到原子的运动轨迹，还能像做实验一样做各种物理现象的观察。对于平衡系统， 可以在一个分子动力学观察时间内做时间平均来计算一个物理量的统计平均值。 对于一个非平衡系统过程，只要发生在个分子动力学观察时间内( 一般为l l o o p s ) 的物理现象也可以用分子动力学计算进行直接模拟。可见数值实验是对 理论和实验的有力补充，特别是许多与原子有关的细微细节，在实际实验中无法 获得，而在计算机模拟中可以方便的得到。这种优点使分子动力学在材料科学中 显得非常有吸引力。分子动力学假定原子的运动是由牛顿运动方程决定的，这意 味着原子的运动是与特定的轨道联系在一起的。当核运动的量子效应可以忽略， 以及绝热近似严格成立时，分子动力学的这一假定是可行的。绝热近似也就是要 求在分子动力学过程中，每一时刻电子均要处在相应原子结构的基态。大多数情 形下，这一条件都是满足的。要进行分子动力学模拟需要知道原子闯正确的相互 作用势，从而必须知道相应的电子基态。电子基态的计算是一个非常困难的量子 多体问题。密度泛函的弓l 入使这方面的计算有了很大的变化。在分子动力学模拟 中一般只能采用经验势来代替原子间实际作用势。一般就用l c n n a r d - - j o n c s 势之 i 地工盐太雪球t 零值凳支 类的二体作用势来替代。这种替代可能对惰性气体之类的系统进行模拟有效，但 是许多实际材料科学中感兴趣的问题都是多体效应的结果，不是只用二体势代替 能解决的。因此人们对金属核共价系统的模型傲了许多尝试，提出了各种多体势 的函数形式。并取得了一些显著的成果，在实际中得到了广泛的应用。 分子动力学研究不受样品制备和测试技术的限制，因此分子动力学模拟能够 找出反应纳米晶体微观结构和力学性质之间的内禀本质，有助于人们对纳米材料 的微观结构与性能之间关系的理解。分子动力学的研究是对实验研究做出理论解 释，补充和弥补实验的不足，为纳米材料的制各新工艺和纳米材料的工业应用奠 定理论基础。因此分子动力学模拟在纳米材料科学研究中得到了广泛的应用。分 子动力学计算机模拟已经成功地被用来模拟纳米纯金属( 比如c u ，n i ，f e ，p a ) 、 非金属( 如s i ) ，合金( 如n i a 1 ) ，陶瓷等。模拟的内容包括晶格畸变，晶粒 生长，弹性模量，拉一压应力一应变关系，蠕变行为，高温变形行为，扩散，沉 积，烧结，固结等等，均取得了比较满意的结果。分子动力学计算机模拟在纳米 材料研究中的作用是相当重要的。随着对纳米材料研究的深入，它的作用有望进 一步扩展。 到目前为止，人们对于纳米金属力学性能和变形机制的认识在计算机模拟方 面主要基于纳米晶纯f c c 金属的研究，对纳米晶b c c 、h c p 金属与合金在拉伸 或压缩载荷下的分子动力学模拟则相对较少。s + r e nl f r e d e r i k s e n l 7 1 等利用分子 动力学方法模拟了纳米品b c c 铝在高应变率下的塑性变形，结果发现，纳米晶 钼在塑性变形过程中不但存在晶界运动和位错的迁移，而且在晶界区有裂纹形 成，这显然不同于纳米f c c 金属的模拟结果。通过分子动力学模拟研究纳米晶 b c c 金属的微观结构的演化与变形机制对于全面认识纳米金属材料具有积极的 作用。 本文利用分子动力学方法模拟了纳米晶铁和钽在单向拉伸载荷作用下的变形行 为，研究了不同晶粒尺寸的纳米晶铁和钽在不同应变率，不同温度下的力学行为 及微观结构的演化，力求对纳米晶b c c 金属的变形机制给出一个合理而又明确 的解释。纳米晶铁的模拟采用康涅狄格大学的j o nr i f k i n 教授发布的x m d 分子 动力学计算软件，纳米晶钽的计算采用森地亚国家实验室发布的大规模分子动力 学软件l a m m p s 。本文主要内容包括以下几个部分。首先是研究方法的介绍，在 彳北工业太簟覃1 - 事幢怡盘 这一部分我们将介绍分子动力学的基本原理，本工作所用到的相关算法以及所采 用的势函数，试样应力的计算方法，纳米晶试样初始构型的构建及算法设计、程 序实现，用于分析原子构型的公共近邻分析方法。其次是纳米晶铁和钽的模拟过 程处理及结果分析，包括平均晶粒尺寸、湿度、应变率等对力学性能的影响，试 样在变形过程中的微观结构演化，如晶间断裂、相变等的出现，晶界结构的变化， 以及对纳米晶b c c 金属微观变形机制的解释。尤其是对纳米晶b c c 金属在变形 过程中出现的不同于纳米晶f c c 金属的晶间断裂和相变现象做了详细的分析。 在结论部分，我们将比较纳米晶f c c 金属与纳米晶b c c 金属的相同点、不同点， 从而力求对纳米晶b c c 金属的力学性能与变形机制给出一个合理的解释。 j 北王t 太事覃t 掌位话杰 一i _ 第二章研究方法及原理 为了模拟纳米晶体在单向拉伸载荷作用下的变形行为，首先要构造和实际晶 体结构相似的纳米晶试样，确定模拟系统的初始构型和组成该构型的各个原子坐 标。其次对该初始构型进行分子动力学驰豫和模拟单向拉伸载荷作用下的加载， 并计算出所需要的物理量。最后要对变形前后的构型进行拓扑分析，确定变形前 后构型原子组态的变化，从而在微观水平上解释所模拟试样的变形机制。根据上 述要求，本章将依次介绍纳米晶试样计算机合成方法，分子动力学原理，应力的 计算方法和公共近邻分析方法。 2 1 纳米晶试样的计算机合成算法 在用分子动力学计算纳米金属及合金的力学性能时，一般用以下三种方法产 生试样的初始构型。( 1 ) 将具有不同晶向的原子簇组装压缩成块状试样限9 1 。( 2 ) 熔化金属结晶法，即在熔化的金属液体中，预先布置若干个取向不同，位置固定 的晶核，随后对金属进行冷却，冷却过程中液态金属原予将会以固定晶核为中心 结晶，从而长大成一个具有不同取向的纳米晶试样1 啦“】。( 3 ) 基于v o r o n o i 图的 空间填充构建方法【2 ，3 ，1 2 , 1 3 , 1 4 , 15 , 4 0 , 4 ”。前两种方法产生的试样内部含有各种各样 的缺陷，从而对模拟纳米晶金属的本质力学性能产生不利的影响。相比之下，第 三种方法产生的试样晶粒内部不含任何缺陷，晶界以高角晶界为主：为了获得织 构试样，可以对晶粒的生长方向 进行控制：合成试样的密度最接 近实验密度( 对一个2 0 h m 的试 样可以达到9 9 2 6 ) 1 16 1 ，因此获 得了广泛的应用。 本文采用基于v o r o n o i 图的空 间填充方法来构建纳米晶粒，下 面给出具体的算法流程图及关键 部分的详细说明。为方便起见， 对原理的说明采用二维v o r o n o i 图，但可以很方便地扩展到三维 情况。 图2 1 程序流程图 j 她z 业太雪硪- j ：，幢话盘 2 1 1v o r o n o ；空间结构图的建立 v o r o n o i 图是计算几何学中的基本结构之一【1 7 】，最近点意义下的v o r o n o i 图是 指在空间给定的点集( p i ，p 2 ，p n ) 中，和其他点相比，距离点p i 最近的 点的区域记为q i ，则集合( q1 ，q 2 ，q n ) 在空间所形成的结构即为v o r o n o i 结构。用v o r o n o i 图思想在计算机模拟中构建纳米晶粒基于如下考虑：晶体在结 晶时会有预先的结晶成核点，结晶时晶体以每个成核点为球心向各个方向生长， 假设各个方向的生长速度相同，则所生长出的晶粒形状就是球形晶粒，当两个晶 粒的生长发生碰撞，在该方向的生长就停止了，直到所有方向的生长都停止下来。 根据v o r o n o i 图的知识，生长出的晶粒的集合形成的三维空闯结构就是v o r o n o i 图。 用计算机产生纳米晶粒构型需要下面几个参数：试样三个方向的尺寸l x ， l y ，l z ，晶格类型( b c c ，f c c ) ，随机种子( 相当于结晶成核点) 的最小允许 距离d ( 用以控制晶粒尺寸的分布) ，试样中包含的晶粒数目n 。图2 1 给出了 计算机合成纳米晶试样的总体流程图，下面就每一步分别给予详细说明( 第一步 省略) 。 2 1 2 随机种子和随机矩阵的产生 产生随机种子和随机矩阵的方法很简单，就是利用随机数发生器产生在( 0 ， 1 ) 区间上均匀分布的随机数。个随机种子具有三个坐标分量，因此要产生三 个随机数构成一个种子的坐标如果产生的在( o ，1 ) 上均匀分布的三个随机数 分别为r x ，r y ，r z 则随机种子的坐标为( l x r x ，l y g y ，l z r z ) ，这并 不难理解。因为试样三个方向的尺寸 不一定是l 。在种子产生之后，要检查 该种子是否满足要求，即该种子与其 他已经产生的种子的距离是否都大于 最小允许距离d ，如果满足要求，则继 续产生下一个种子，如果不满足，则 该种子无效，重新产生该种子的坐标， 直到其满足要求为止。注意，在产生 种子坐标时，如果d 和n 都比较大，程图2 2 髓机种子产生流程图 否北z t 太簟巧t 覃位话盘 序找不到满足这种要求的随机点集，就会陷入死循环。这种情况下就要终止程序， 减小晶粒数目n 或者最小允许距离d 到一个合适的范围后重颏运行该程序。 要产生方向随机分布的晶粒，需要对每个晶粒产生随机矩阵以控制晶粒的生 长方向。一个随机矩阵( 三行三列) 具有九个分量，因此需要产生九个随机数。 由于所器要的矩阵为坐标交换矩阵，在九个随机数产生之后，需要对产生的矩阵 进行正交化，正交化采用s e h m i d t 方法，详 2 1 3 晶粒生长 在用v o r o n o i 图产生纳米晶粒的时 候，要解决两个关键问题。一是纳米晶 粒如何生长，按什么顺序生长。二是如 何判断产生的点是否超出了边界。上述 y 两个问题考虑的着眼点是程序实现的 方便程度，计算机的执行效率和占用内 存的情况。 晶粒的生长采用正方体逐层生长 i 北工盐太簟壕 季值铪盘 外，该晶粒的生长即告完成。假设第i 个种子所产生基准点的局部坐标用l = ( h ，k 。1 ) 7 表示，原子全局坐标用 r = ( x ，v ，z ) 1 表示，第i 个种子的全 局坐标为r i = ( x i ，y i ，z i ) 1 ，q 为一 随机产生的正交矩阵，则对f c c 金属。 以此基准点所产生的原子坐标为l 0 = ( h ，k ，1 ) t ，l 1 = ( h k + 0 5 ，1 + 0 5 ) t ，l 2 = ( h + 0 5 ，k ，1 + 0 5 ) 7 ，l 3 = ( h + o 5 ， k + 0 5 ，1 ) t 对b c c 金属所产生的原子 局部坐标为l 0 = ( h ，k ，1 ) tl 1 = 图2 4 原子坐标判断示意图 ( h + 0 5 ，k + 0 5 ，1 + 0 5 ) to 然后利用下式进行坐标变换，式中l 为原子的局部坐 标。 r = l u + q l 对所产生的每一个原子坐标都要判断该坐标是否属于该晶粒，判断方法参考 图2 4 。图2 4 为一个v o r o n o i 图结构( 不包括以p 为圆心的圆) ，p 点所在的多 边形区域称之为p 区，同理，有p 1 ，p 2 区等等。将p 的最近邻种子按由近及远 的顺序排列，分别为p 1 ，p 2 ，p 3 ，p 4 ， p 5 ，p 6 。并记下各点到p 点距离的一半， 分别记为r l ，r 2 ，r 3 ，r 4 ，r 5 ，r 6 。图中 圆l ，2 ，3 ，4 ，5 的半径分别为r l ，r 2 ， r 3 ，“，r 5 ，可以采用如下的判别方法。 首先计算生长点到p 的距离d ，判断d 是否小于r l ，如果是的话，则该点在1 号圆内，显然落在p 区。若d 大于r l ， 则计算生长点到p 1 的距离d l 。若d 1 小 于d ，说明生长点已经不属于p 区，如 果d l 大于d ，则比较d 和r 2 ，若d 小于 r 2 ，则生长点位于2 号圆环域内，从圈 图2 5 判断算法 叠她工t 太掌硪t 事位论立 上可以看出，2 号圆环域跨越p 区和p 1 区。由d 1 大于d 可以判断生长点不在p 1 区，则只能落在p 区，若d 大于r 2 ，则计算生长点到p 2 的距离d 2 ，若d 2 小于d ， 则生长点不在p 区，若d 2 大于d ，则比较d 和r 3 ，若d 小于r 3 ，生长点位于3 号圆环域内，由前述分析知道生长点不在p l 区和p 2 区，则只能落在p 区，若d 大于r 3 ，则计算生长点到p 3 的距离，依次类推，直到判断出生长点是否属于p 区。该算法的计算量较小，不用计算生长点和每个种子的距离，只需做若干次比 较即可。程序的运行结果表明，产生一个含有2 7 0 个晶粒，约3 2 0 0 万原子的试 样，所需要的c p u ( 奔四2 8 g h z ) 时间大约为1 4 分钟。这对一个需要运行若干 天的m d 模拟说来，该时问消耗是微不足道的。流程图( 图2 5 ) 更条理地说明 了该判别方法。 2 1 4 晶粒的修剪 在晶粒生长的过程中，由于没有对生长出的原子坐标进行碰撞检测，边界附 近生长的原子坐标有可能距离太近甚至发生重叠。当两个原子距离太近，小于 m d 算法的最小允许距离时，就会有可能导致某个原子的电子密度超过程序所允 许的范围而使分子动力学程序意外终止。为了避免这种情况的发生，有必要对产 生出的纳米晶试样进行修剪消除原子距离太近甚至重叠的情况。具体采用的方 法是，如果两个原予距离太近甚至重合，则删除其中一个原子。这里可能要对所 有产生的晶格点进行搜索，并计算两两之间的距离，这种方法是最费时的方法。 一个比较高效的算法是先对所有的晶格点的坐标进行排序，按x y z 从小到大的顺 序排列，或者根据每一个格点的坐标就可以知道该点的位置，并能很快的找到该 点的最近邻原子点坐标，一般情况下这种最近邻原子点是不多的，然后计算这些 最近邻与该点的距离，如果超过了截断距离，则将该原子删除，否则，保留该原 子，继续进行下一个原子点的检查。 假设原子间的最小允许距离为d ，对一个给定的原子，在以该原子为球心，d 为半径的范围内是不应该存在其他原子的，否则要将该原子删除。不妨设d = 0 7 ， 为了程序上的便于实现，设置一个三维的数组x y z m l n l 阳，并假设盒子的尺寸 为l x ，l y ，l z ，则m = i m ( l x * 1 7 3 2 0 7 ) ，这样做的目的是每个数组里面只存储 一个坐标，如果有两个离的太近的坐标，则坐标取整以后有可那相同，这样先前 存储的坐标就会被后一个坐标覆盖，从而起到第一遍的过滤作用，这样就会自动 1 0 摹北工叠太簟覃士事传话立 将不符合条件的坐标删除，为了和数组的大小相吻合，每个坐标分量也要乘以一 个因子1 7 3 2 0 7 = 2 5 ，然后取整，这样根据每个索引就可以大致判断该数组中的 坐标大约为多少，如索引x y z i j 闯中的坐标大约为x = v 2 5 ，y = j 2 5 ，z = k 2 5 。 在读入坐标之前，先将这些坐标初始化为一1 0 0 ，因为所有原子的坐标都是大于 零的。而且可能有的数组根本没有用上。在将坐标读入以后，就要进行判断看是 否有的原子距离太近，这时候对每一个原予坐标只需要判断2 6 个近邻即可，即 该坐标上下左右前后的坐标，这些坐标可以根据索引很容易的找到。对于索引为 独的原子坐标，可能离其太近的原子坐标存放在索引i - l j k ，i j - l k ，i j k 一1 ，i + l j k ， i j + l k ，i j k + ，i - l j 1 k 等2 6 个数组元素里，因此只需要判断这几个坐标是否距离 给定的原子太近即可。如果存在距离太近的原子，则将该中心原子的坐标的三个 分量置为一1 0 0 ，相当于将该原子删除掉。最后一步就是把符合条件的坐标输出 到文件中，符合条件的坐标是指那些坐标分量不是一1 0 0 的坐标。这种搜索方法 是一种比较省力省时的方法，其运行效率应该是比较高的。但是如果原子数目比 较多的话，该算法占用的内存比较大，空间复杂度较高。 图2 6 b c c 纳米晶试样 图2 7 f c c 纳米晶试样 图2 6 是产生的纳米晶粒试样( b c c ，2 7 0 个晶粒，约1 6 0 0 万个原子) 效果 图，根据原子的中心对称参数【1 s l 赋予不同的颜色。灰色代表晶粒内原子，黑色代 表晶界原子。图2 7 为产生的f c c 纳米晶试样( 2 7 0 个晶粒，原子数目约3 2 0 0 万) 。 在该算法的基础上通过对种子的位置进行控制，即不使用随机数发生器产生 种子的位置，而是指定种子的位置，就可以产生各种形状的品粒，如六面体晶粒， 叠北工业太事摹士簟佳话t 八面体晶粒和其他形状的晶粒。通过控制随机矩阵以控制纳米晶粒的生长方向， 就可以产生织构试样。对该算法进行适当的修改还可以产生各种形状的试样，以 模拟纳米团簇的各种性能。因此该算法可以很方便地对晶粒的尺寸及分布，晶粒 形状及方向进行控制，从而产生出满足一定条件的模拟试样，是一种行之有效且 简单易行的纳米晶试样构建方法。 2 1 5 产生满足周期性边界条件的构型 在处理原子、分子的聚合体问题时，分子动力学能处理的原子数目要受到计 算机运行速度和储存能力的限制，为了模拟块体材料，人们引入周期性边界条件。 因此在产生纳米晶试样时也要考虑到周期性边界条件。图2 6 、2 7 所示试样并 没有考虑到周期性边界条件。为了产生满足边界条件的试样需要对上述算法做适 当的技术处理。 以二维情况为例，满足周期性边 界条件的纳米晶试样如图2 81 9 , 所 示，其周期性单元已经标出，其中 红色线框所示为基本单元，黑色线 框所示为复制单元。线框内原点为 随机种子。为了产生满足周期性边 界的试样，首先在寻找给定种予的 近邻种子并按距离排序时，要考虑 圈2 8 满足餍期牲边界条件的试样 到其周期性复制单元内的种予坐标。但在晶粒生长时只以基本单元内部的种子为 中心。按此方法生长出的晶粒形状如图2 8b 所示。根据给定的试样三个方向的 尺寸，当某些原子的坐标大于该方向的尺寸l ，令其对应坐标减去l ，当某个方 向的坐标小于零时，令其坐标加上l ，这样就能把所有原子的坐标按照周期性关 系平移到所规定的尺寸范围，如图2 8c 所示。另外在晶粒修剪时也考虑到周期 性边界的情况。不光要保证试样内部的原子两两之间的距离大于给定值，还要保 证基本单元和复制单元边界之间的原子也要满足这个条件。 2 2 分子动力学原理 目前已经有许多文章和著作详细介绍了分子动力学方法的原理及其实现方 法。因此仅仅出于完整性的需要，本文只是简单介绍一下和本文工作相关的分子 圈圈 罩她工口太事项t ，幢怡盘 动力学的知识。对于势函数的介绍也仅仅局限于内嵌原子方法，其他类型的势函 数以及基于第一原理的相关知识，这里不做介绍。具体详细的分子动力学的知识 可参考相关文献 1 9 - 2 3 】。 2 2 1 牛顿运动方程 分子动力学模拟是指对于由多个原子核和电子所组成的多体系统，用计算机 模拟每个原子核的运动轨迹，通过对原子的运动轨迹进行统计力学分析计算我们 感兴趣的物理量。假设一个含有n 个原子的体系，体系的状态由这n 个原子的位 置 r t 和动量 p ；) 或速度 v i 来标志。对于三维问题，每个原子有3 个坐标分量 ( x ，y ，z ) ，三个速度分量( v 。，v ，v ：) ，这样n 个原子的这六个分量组成6 n 维 相空间。体系的能量为位置和动量的函数，表示为h ( r i ，p i ) 。体系的运动方 程为： 旦口：一旦 瓦只一石爿 旦忙旦h 西。 印 在实际应用中，把上面的哈密顿方程化为下面的牛顿方程，并且用位置r 和速度 v 。做为描述体系的参量。 h = 去m ，v f 2 + 矿( ) ) m a 要：一昙 矿一丽”( ) v ( f r c 为原子间相互作用势能，通过求解上面的方程就可以得到体系在相空间 的轨迹，迸而求得我们所感兴趣的物理量的统计平均值。分子动力学模拟的关键 问题是原子间作用势的确定，核心问题是求解牛顿运动方程组，运动方程的求解 一般采用有限差分方法，如v e r l e t 算法、l e a p - f r o g 算法等，具体算法的差分 格式可参考相关书籍，这里不再赘述。 2 2 2 分子动力学的分类 根据对原子问作用势不同的简化处理方法，分子动力学可分为：多体分子动 力学( m a n y b o d ym o l e c u l a rd y n a m i c s ，m b m d ) 、可变电荷分子动力学 蕃玷z 业太fj 螽士簟住倍盘 ( v a r i a b l e c h a r g e m o l e c u l a rd y n a m i c s ，v c m d ) 、紧束分子动力学( t i g h t b a n d i n g m o l e c u l a rd y n a m i c s 。t b m d ) 、密度泛函分子动力学( d e n s i t y - f u n c t i o n a lm o l e c u l a r d y n a m i c s ，d f m d ) 、第一原理分子动力学( f i r s t - p r i n c i p a lm o l e c u l a rd y n a m i c s ， f p m d ) 。其中前两种属于经典的分子动力学方法。它通过实验结果或者经验模型 确定原子问的作用势，优点是计算量小，适用于解决较大规模的问题，缺点是可 移植性差，适用问题的范围小，不同的问题往往需要发展不同的势函数。后面三 种属于现代分子动力学方法，又称为从头计算分子动力学方法( a bi n i t i om o l e c u l a r d y n a m i c s ，a i m d ) 。这种方法直接从量子力学的轨道理论出发获取原子间的作用 势，具有精确度高，可移植性好的优点，但是由于该方法要考虑到电子的运动特 点，计算量很大，从而所能求解的问题的规模和经典分子动力学方法相比要小很 多。 2 2 3 势函数 在积分牛顿方程时需要每个原子所受到的作用力，该力通常通过势能函数的 梯度获得，这取决于粒子的位置。模拟的真实性因此取决于所选择的势在模拟运 行所选择的条件下能重新产生材料行为的能力。因此原子间作用势是影响分子动 力学模拟精度的主要因素，在所研究问题的范围之内，原子间作用势要能反映原 子间相互作用的真实情况。原子间势函数的发展经历了从对势到多体势的过程。 对势认为原子之间的相互作用是两两之间的作用，与第三个原子的位置无关，但 是实际上任何一个原子的存在都会影响其周围原子的电子云分布，从而影响原子 间的相互作用力，因此虽然对势具有简单明了，计算量小，所适用问题规模较大 的优点，但其无法考虑化学键的方位变化，因此对固体材料结构的剪切、弯曲、 扭转等变形很难描述。为此提出了多体势的概念。 在纳米金属的分子动力学模拟中用的最多的原子问作用势为内嵌原子势和 有效介质理论模型，这些势函数都可用下式来表达 e = 妒( 勺) + f ( 岛) j - ，i = 1 其中第一项为对势项，根据需要可以采取不同的形式。f 是嵌入能。n 是除第f 个原子以外的所有其他原予的核外电子在第i 个原子处产生的电子云密度之和， 可以表示为 誓她z 业太簟磺t 簟幢伦互 一i i i i _ - p ，= ，( ) j l 厂( ) 是第- ，个原