（机械制造及其自动化专业论文）纳米切削过程的分子动力学及多尺度仿真研究.pdf

啥尔滨t 业大学1 二学硕士学位论文 摘要 随着机械加工精度要求的不断提高，超精密加工技术得到了不断发展， 切削厚度已达到纳米量级。由于尺寸效应，适用于宏观尺寸的一些切削机理 不再适用，势必需要一种基于微观力学的理论来对纳米切削机理做出解释。 分子动力学方法恰具备此特点，它从构成物质的基本粒子出发，具有沟通宏 观特性和微观结构的作用，对许多理论分析和实验观察上都难以了解的现象 做出一定的微观解释。因此，采用分子动力学来研究纳米切削机理，受到了 越来越多的关注。 本文简要介绍了分子动力学在纳米切削过程仿真中的应用，以及分子动 力学未来的发展趋势。对分子动力学的基本计算原理进行了讨论，主要包括 势函数、数值积分方法、温度控制及边界条件等。 实现了一套完整的分子动力学仿真软件。本软件不仅具有分子动力学仿 真的基本功能，而且可以同时应用多种势函数，适用于多种物质相互作用的 模拟；可以根据需要使用不同精度的数值算法求解运动方程；使用参数化的 模型建立方式，并提供可选的数据输出；采用模块化设计，便于将来扩展使 用。 应用该软件对单晶铜进行了切削仿真实验，重点对晶体晶向对切削过程 的影响进行了仿真研究。 从原理、特点、应用等方面对多尺度仿真方法进行了探讨。该方法结合 了分子动力学和有限元原理，不仅具有分子动力学仿真计算精确的特点，而 且可以跨越从宏观到微观的多个尺度，计算量也相对较小，扩展了分子动力 学的应用范围。 论文最后应用多尺度仿真方法，对单晶铜的纳米切削过程进行了仿真研 究，并与分子动力学的切削仿真结果进行了比较。 关键词 分子动力学：多尺度；仿真；纳米切削 a b s t r a c t a l o n gw i t ht h em a c h i n i n gp r e c i s i o nr e q u e s t su n c e a s i n ge n h a n c e m e n t ，t h e u l t r a p r e c i s ep r o c e s s i n gt e c h n o l o g y i s d e v e l o p i n gc o n s t a n t l y ，t h e t h i c k n e s so f c u t t i n gh a sa c h i e v e da n a n o m e t e rm a g n i t u d e a sar e s u l to ft h es c a l ee f f e c t ，s o m e c u t t i n g m e c h a n i s m s f o r t h e m a c r o s c o p i c s i z ea r en o l o n g e r s u i t a b l e a m i c r o s c o p i c m e c h a n i c s t h e o r y i sn e e d e dt o e x p l a i n t h en a n o m e t e r c u t t i n g m e c h a n i s m t h em o l e c u l a rd y n a m i c s ( m d ) w h i c hr e s e a r c h e st h ef u n d a m e n t a l p a r t i c l ep o s s e s s e st h e s ep r o p e r t i e s i tc a l ln o to n l yc o n t a c t t h em a c r of e a t u r e sa n d m i c r os t r u c t u r e s ，b u tg i v e su n i q u ee x p l a n a t i o n st os o m es p e c i a lp h e n o m e n at h a t a r eh a r d l yk n o w n b y m a c r o t h e o r y ，t h e r e f o r e ，t h ec u t t i n gm e c h a n i s m r e s e a r c hb y m o l e c u l a rd y n a m i c sh a sr e c e i v e dm o r ea n dm o r ea t t e n t i o n t h i st h e s i s b r i e f l y i n t r o d u c e dt h em o l e c u l a rd y n a m i c si nt h en a n o m c t r i c c u t t i n gp r o c e s ss i m u l a t i o na p p l i c a t i o n ，a sw e l la si t sf u t u r ed e v e l o p m e n tt e n d e n c y d i s c u s s i o nw a sc a r r i e do nt h em o l e c u l a rd y n a m i c sb a s i cc o m p u t a t i o np r i n c i p l e i t m a i n l yi n c l u d e dt h ep o t e n t i a lf u n c t i o n ，t h en u m e r i c a li n t e g r a t i o nm e t h o d ，t h e t e m p e r a t u r ec o n t r o la n d t h eb o u n d a r yc o n d i t i o n ，e t c a c o m p l e t e m o l e c u l a rd y n a m i c ss i m u l a t i o ns o f t w a r ew a sd e v e l o p e d b e s i d e s t h ee s s e n t i a la b i l i t i e s ，t h i sp r o g r a mi sa b l et o a p p l ym a n yk i n d so fp o t e n t i a l f u n c t i o n sa to n et i m e ，s u i t a b l ef o rt h ei n t e r a c t i v es i m u l a t i o no fm a n yk i n d so f m a t e r i a l s s e v e r a in u m e r i c a lc a l c u l a t i o nm e t h o d sw i t hd i f f e r e n tp r e c i s i o na r e c o v e r e da n do t h e rf e a t u r e ss u c ha sp a r a m e t e r i z e dm o d e lb u i l d i n g ，o p t i o n a ld a t a o u t p u ta n dm o d u l ed e s i g n f o re x p e d i t i o na r ea l s oc o n t a i n e d n a n o m e t r i cc u t t i n gs i m u l a t i o n sf o rs i n g l e c r y s t a lc o p p e rb ym d sw e r e i m p l e m e n t e d t h ei n f l u e n c eo fc r y s t a lo r i e n t a t i o no np r o c e s s i n gw a se m p h a s i z e d p a r t i c u l a r l yo n t h em u i t i s c a l es i r e u l a t i o nm e t h o dw a sd i s c u s s e di nt e r m so fp r i n c i p l e c h a r a c t e r i s t i c ，a n da p p l i c a t i o n t h i sm e t h o dh a sc o m b i n e dm o l e c u l e rd y n a m i c s a n df i n i t ee l e m e n tm e t h o d i tn o to n l yh a st h em d sp r e c i s i o n ，b u te x c e e d si t s r a n g eo fa p p l i c a t i o nf r o mm i c r o s c a l et om a r c r o s c a l ea n do w n ss m a l lc a l c u l a t i n g a m o u n tr e l a t i v e l y , e x p a n d i n gt h er a n g eo fa p p l i c a t i o no fm o l e c u l a rd y n a m i c s t h et h e s i sp e r f o r m e das i m i l a rc u t t i n gs i m u l a t i o nt om o l e c u l a r d y n a m i c so n i i 哈尔滨工业大学工学硕士学位论文 s i n g l ec r y s t a lc o p p e rb yt h em u l t i s c a l em e t h o df i n a l l ya n dc o m p a r e dw i t ht h e c u t t i n gr e s u l to f m o l e c u l a rd y n a m i c s k e y w o r d s m o l e c u l a r d y n a m i c s ，m u l t i s c a l e ，s i m u l a t i o n ，n a n o m e t r i cc u t t i n g 。，。丝鎏些塑羔塑篓鍪些。，一 1 1 课题的背景和意义 第1 章绪论 机械工业是国民经济发展的基础，因为它需要为其他生产部门提供技术装 备【”。随着科学技术的不断进步，近几十年机械制造技术在提高精度等方面发 展迅速，已从精密加工发展到超精密加工，其加工精度从微米级提高到亚微米级 乃至纳米级，对国防【2 】、航空航天、核能开发及微电子等高新技术领域有着重 要的影响。例如，大规模集成电路的硅基片、计算机的高密度硬盘、天体望远 镜系统的反射镜和多面棱镜、导弹的陀螺仪等都是超精密加工的产品。精密及 超精密加工技术的发展和推广提高了整个机械制造工业的加工精度和技术水 平，使机械产品的质量、性能及可靠性都得到了普遍的提高。超精密加工技术 已经成为目前先进制造领域的前沿课题，它不仅仅是门科学技术，更是一个 国家科学技术水平和综合国力的藿要标志，困此受到各工业发达国家的高度重 视f 3 l ，它已成为各国在国际竞争中取得胜利的关键技术之一。 目前，超精密加工所依赖的各种技术如机床制造技术、测量和误差补偿技 术、金刚石刀具的刃磨技术都达到了相当高的水平，超精密加工正逐渐向纳米 级、原子级加工精度逼近1 4 j 。1 9 8 7 年，美国劳伦斯实验室和日本大阪大学的学 者在超精密机床p e r l i i 上对电镀c u 进行的直角切削实验中证实，使用金刚石车 刀可以实现数个纳米的切削厚度1 5 j 。但是，超精密加工中的一些机理还没有被 明确认识，这极大地阻碍了超精密加工技术的进一步提高。 与传统加工过程相比，超精密切削加工发生在很小的区域，该区域只包含 数个原子层或数百个原子层。加工过程在本质上是原子的离散现象，工件材料 应看作是原子或分子的集合体，对加工过程采用建立在传统连续介质力学基础 上的切削理论来解释显然是不合适的。此外，当被加工工件的精度达到纳米量 级时，会出现一些特殊的物理现象如：微磨损、尺寸效应等。对这些现象必须 从微观的角度加以研究。 根据以往的经验，通过切削实验来获取切削规律是最好的办法。但通过实 验研究纳米切削规律，不仅需要进行大量的纳米级加工实验，成本很高，而且 当材料切削厚度达到几个原子层时，实验变得困难且耗时。因为实验对超精密机 床、检测仪器分辨率、切削条件、刀具及其几何形状的要求极其苛刻，目前的 哈尔滨工业大学工学硕士学位论文 实验条件难以或根本无法实现。而分子动力学( m o l e c u l a r d y n a m i c s ，m d ) 模拟方 法正好能克服这些困难。它的基本原理是：建立一个粒子系统来模拟所研究的 对象，系统中各粒子问的相互作用根据量子力学来确定。对于符合经典牛顿力 学规律的大量子系统，通过粒子动力学方程组的数值求解，决定各粒子在相空 间的运动规律和轨迹，然后按照统计物理原理得出该系统相应的宏观物理特性。 将分子动力学方法应用在超精密切削机理的研究中，可以弥补实验方法的 不足，并且模拟结果有相当的应用价值，因此，近些年来已有不少学者采用浚 方法对超精密切削中的一些现象进行了模拟研究，并且一致认为分子动力学模 拟方法是分析微切削过程十分行之有效的工具。 但是，目前尚无一种高效通用的分子动力学模拟软件供不同研究方向的专 业人员使用。如果想采用该方法就必须自己动手编制程序，抛开计算机本身的 因素，对于非计算机专业的人员需要耗费大量时间和精力来开发可供专业使用 的模拟程序，但结果经常是事倍功半，即使程序勉强可以运行，也是效率低下， 漏洞百出，并且实用性很差。 实现一种可以在机械相关领域通用的分子动力学仿真软件，正是本文的目 标之一，通过使用该软件，使需要进行仿真计算的人们，可以把工作重点从编 程转移到仿真本身上，扩展分子动力学技术的应用空间。 1 2 分子动力学仿真的国内外研究发展综述 1 2 1 分子动力学仿真方法发展简介 近年来，利用计算机模拟技术研究材料的力学行为，日益成为人们感兴趣 的课题。由于计算机处理速度的迅速提高，计算机模拟已经和实验观察、理论 分析并列成为本世纪科学研究的三种方法f 。计算机模拟方法可以用来比较、 验证各种近似理论；也可以用来对实验和模型进行比较，从而提供了评估一个 模型正确与否的手段。计算机模拟方法的优点就是可以沟通理论和实验。某些 量或行为可能是无法或难以在实验中测量的，而用计算机模拟方法，这些量可 以被精确的计算出来。目前已经发展了几种可行的计算机模拟方法，如分子动 力学法、蒙特卡罗方法( m o n t ec a r l o ，m c ) 、最小能量法( e n e r g ym i n i m i s a t i o n e m 、 和晶格动力学法( l a t t i cd y n a m i c s ) t t l 。 分子动力学模拟是指对于原子核和电子所构成的多体系统，用计算机模拟 原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在 啥尔滨工业大学工学硕士学位论文 全部其他原子核和电子所提供的经验势场( 势函数) 作用下按牛顿定律运动i 8 。分 子动力学方法以其建模简单、模拟结果准确的特征而倍受研究者们的关注。而 且可以说分子动力学方法是电子计算机的孪生物，5 0 年代初第一台电子计算机 诞生后不久，在著名的物理学家f e r m i l 9 】的建议下就开始了分子动力学研究，当 时模拟对象只是1 6 个粒子组成的一维系统。在早些时候一些计算机模拟基本 方法( 如周期性边界条件、截断势等) 也已经被提出。 6 0 年代初，a l d e r 署1 w a i n w r i g h t t t o l 应用分子动力学研究了硬球系统，这被公 认为是第一次分子动力学模拟，这也是第一次将分子动力学应用到物理和化学 物理领域。v i n e y a r d “j 等人研究了液体中粒子的连续作用势，这是分子动力学 模拟在材料科学上的第一次应用。v e r t e t t ”j 予1 9 6 7 年绘出了l e m a e s d - i o u e s 系统 相图，并提出了著名的v e r l e t 算法，这使得分子动力学的算法更有效，分子动力 学应用范围得到扩展。 进入7 0 年代以后，分子动力学得到了更加蓬勃的发展，涨落耗散理论 ( f l u c t u a t i o nd i s s i p a i o nt h e o r y ) 、密度泛函理论( d e n s i t y - f u n c t i o nt h e o r y ) 等新理 论纷纷被引入分子动力学，同时它又与统计物理等学科结合，使得经典分子动 力学得到充实，新的算法不断涌现。 到了8 0 年代，分子动力学的体系逐渐完善，大体可分为：平衡态分子动力 学( e q u i l i b r i u mm o l e c u l a rd y n a m i c ) 、非平衡态分子动力学( n o n e q u i l i b r i u m m o l e c u l a r d y n a m i c s ) 、从头计算分子动力学( a b i n i t i o m o l e c u l a r d y n a m i c s ) 和大规 模分子动力学( l a r g es c a l em o l e c u l a rd y n a m i c s ) t ”0 6 1 。 从9 0 年代开始，计算机技术的飞跃，极大的推动了分子动力掌方法的发展。 随着算法的不断改进，并行计算技术的应用，加之它具有沟通宏观特性与微观 结构的作用，对于许多在理论分析和实验观察上都难以理解的现象均可以做出 一定的解释，因此在物理、化学、生物学、材料科学、摩擦学等领域，都得到 了很好的应用【1 7 - 2 0 o 在超精密加工机理的研究中，分子动力学模拟作为一个行 之有效的工具也正在被广泛使用。 1 2 2 超精密切削机理的分子动力学仿真概况 由于分子动力学仿真可以直观的模拟超精密加工过程，研究加工机理，目 前世界上很多学者对这一领域的研究都做出了较大的贡献。其中美国和日本在 这方面的工作处于世界领先水平，下弱对他们的工作做简单介绍。 美国的劳伦斯实验室最先用分子动力学模拟研究了超精密加工机理。他们 眙尔演工业大学工学硕士学位论文 将分子动力学模拟技术用于微磨损、微压痕与微硬度的研究“，1 9 8 9 年h o o v e r 和s t o w e t s 等人发表了用分子动力学模拟研究单晶c u n z 表面的微摩擦问题的论 文【2 3 】。此后不久，b e l a k 和s t o w e r s 以单晶c u 为研究对象，进行了压痕和切肖4 过程 的分子动力学模拟，用传统的m o r s e 势函数计算分子间的作用力，建立了c u 的纳 米切削过程的模型，如图1 1 所裂2 4 j 。在这个模型中，金刚石刀具的切削刃钝圆 半径为0 5 n m ，被吃刀量为2 2 个原子层( o 5 3 n m ) ，被吃刀量与切削刃钝圆半径 的比值接近1 ：1 。可见，这是一个典型的纳米切削模型。这个模型被人们广泛采 用，成为第一代纳米切削过程的分子动力学模型。b e j 越恿过这个模型研究发现 在己加工表面和后刀面之间的第三剪切区在超精密加工过程中起着很重要的作 用。当被吃刀量与切削刃钝圆半径在一个数量级时，有效刀具前角和合成耕犁 力的变化对整个力系统有很重要的作用。 t u t t i r , d ：l 口t n d e p m o f e u t ：2 2a t o m i ci q 州0 j 3 n r a ) 图1 - 1 单晶铜纳米切削分子动力学模型拉4 j f i g ，1 一it h em d m o d e lo f u a n o m e t r i ee u t t i t , , go f s i u g l ec r y s t a leu f 2 4 进入9 0 年代后，劳伦斯实验室的b e l a k 和s t o w e r s 在分析了金冈0 石刀具压痕 实验和切削实验的表面微观形貌后，对金刚石一铜( 硅) 界面进行了压痕和切 削过程的分子动力学模拟随2 6 l ，用e a m ( t h ee m b e d d e d a t o mm e t h o d ) 方法计算 了铜原子的内聚力，用结合序列模型计算金刚石刀具中的碳原子与硅表面的共 价作用，并提出了新的边界条件，如图1 2 所示。 该边界条件采用周期性边界条件，将模拟的单元分为三个区域：( 1 ) 牛顿区， ( 2 ) 热影响区，( 3 ) 固定边界区。牛顿区的原子运动由牛顿方程加以描述；热影 响区的引入是为了使系统温度保持室温；边界区的原子始终保持不动，也没有 哈尔滨工业大学工学硕士学位论文 原子振动。这种边界条件为新一代分子动力学模拟模型的形成奠定了基础。 图1 2 分子动力学模拟边界条件 f i g 1 - 2b o u n d a r yc o n d i t i o no f m o l e c u l a rd y n a m i c ss i m u l a t i o n f 2 5 在这之后，美国的o k l a h o m a 州立大学的r k o m a a d u r i 教授指导的学生在此基 础上进行了单晶硅切削过程的分予动力学模拟1 2 ”，分子间作用力的计算采用的 是m o r s e 势函数的计算方法。他们的模拟在主频为3 3 3 h z 的口工作站上进行，研 究的重点是刀具几何形状对纳米加工的作用。通过研究，发现大切削钝圆半径 ( 相对于被吃刀量来说) 和负前角在硅等脆性材料的纳米加工中对塑性变形区的 产生有益，有助于硅以塑性方式去除。 日本的分子动力学的模拟研究几乎是与美国劳伦斯实验室同时开展的。 1 9 9 0 年，大阪大学的i k a w a 教授、s h i m a d a 教授和n a g o y a j 墼i 大学的i n a m u r a 教授 都发表了应用分子动力学模拟研究超精密加工机理的论文 2 8 , 2 9 i ，在日本的超精 密加工领域，开创了分子动力学模拟的先河。他们的模型基本上与b e l a k 等人提 出相同，属于第一代分子动力学模型。此后，他们对模型边界条件进行了深入 的研究，s h i m a d a 教授和i k a w a 教授于1 9 9 2 年推出了第二代分子动力学模型，如 图1 3 所示口吼，在这个模型中，他们将刀具原子和工件原子都分为牛顿层、恒温 层和边界层，各层原子的意义与b e l a k 等人提出的是一致的，模拟时对恒温层的 原子标度以保持温度的恒定。对原予间作用力的计算仍采用m o r s e 势函数的计算 方法。他们模拟的对象是6 0 7 6 个工件原子，当时计算了1 7 0 个小时。 从1 9 9 2 年以后，s h i m a d a 教授、i k a w a 教授和i n a m u r a 教授应用分子动力学 模拟对超精密加工机理的研究做了大量工作，推动了这项研究的发展，也取得 了显著的成果。 。：ii i ； 堕玺鋈三i 些奎耋i m i i i 三i i i i i 耋耋窑耋垒：鎏兰：。，：，：。： ! ! = ! ! 自= e = ；一= 女= = = = e = ! ! = ! = e ! ! = = = = = = = = 图1 3 纳米切削分子动力学模型【3 0 l f j g i - 3 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 m o d e l o f n a n o m e t r i cc u t t i n g 【3 0 】 国内的分子动力学模拟研究工作起步较晚。1 9 9 3 年，温诗铸院士领导的清 华大学摩擦学国家重点实验室，开始应用分子动力学对纳米摩擦进行了仿真研 究【2 0 i 。1 9 9 8 年天津大学于思远教授指导的林滨博士，开始应用分子动力学对单 晶硅的超精密磨削过程进行模拟研究c 3 1 , 3 2 1 。从2 0 0 0 年开始，哈尔滨工业大学的 梁迎春教授指导的罗熙淳博士、霍德宏博士、唐玉兰博士等，就硅、铜、铝等 材料的超精密切削机理进行了分子动力学仿真研究，并取得了一定的成果【3 3 3 6 1 。 这些工作为我国的分子动力学模拟研究的进一步开展，奠定了良好的基础。 1 2 3 分子动力学的未来发展方向 从以上介绍的分子动力学模拟研究的发展不难看出，从1 9 9 0 年到1 9 9 5 年 是分子动力学飞速发展的黄金时代，每年有大量的论文在c i r p 、j s p e 和b u l l e t i n 上发表，但是从1 9 9 6 年以后，其发展速度放慢了，这方面的文章数量很少，而 且基本上还是延续以前的方法，对第二代分子动力学模型没有大的改进。就目 前而言，分子动力学模拟的发展方向主要有以下几个方面： ( 1 ) 基础理论的研究开展分子动力学模拟，如果不借鉴量子力学、统计力 学、固体物理等学科的研究方法和最新研究成果，没有系统完备的现代物理知 识，很难将这项研究深入进行下去。 ( 2 ) 模型建立的细化对于不同的材料应该根据他们的不同属性( 共价键、金 属键还是离子键) ，来合理选择描述分子间作用的势函数，这样才能正确地描述 材料本身的特征，不能一味的使用简单的m o r s e 、l e n n a r d - j o n e s 等势函数，这将 哈尔滨工业大学工学硕上学位论文 使模拟结果的可信度受到质疑。 ( 3 ) 计算速度的提高提高分子动力学的计算效率有两条途径，一是改进算 法，这依赖于理论基础的研究；再就是利用先进的计算机技术进行并行计算， 这是目前提高计算量和计算速度的最有效途径。 ( 4 ) 发展多尺度模拟方法单纯的分子动力学方法能够模拟的原子数量毕竟 有限，即使采用了并行计算技术也不能达到模拟真实物理系统的目的。而采用 分子动力学与有限元方法结合的多尺度模拟方法正是解决这一问题的有效途 径，它不仅保持了分子动力学的计算特点，而且扩大了模拟尺寸，拓展了分子 动力学的应用范围，目前这种方法的研究还处于起步阶段。 ( 5 ) 研究范围的扩展目前分子动力学的应用范围还很狭窄，仅局限于结构 完整的单质、单晶体的模拟研究，而实际生产中材料以多晶体材料为主( 3 a 1 ，而 且多存在各种缺陷，模拟各种结构复杂的晶体及化合物是分子动力学有待解决 的问题之一。 l - 3 本课题研究内容 分子动力学模拟方法的可应用领域在不断扩大，但目前还没有像有限元分 析软件那样成熟的分子动力学模拟软件，供各领域中的学者们使用。如果要利 用分子动力学方法来进行模拟，首先必须编制复杂的模拟程序，这不仅要耗费 大量的时间，而且对于非专业编程人员来说，要想得到即功能完善，执行效率 又高的程序是相当困难的。 本文将详细介绍分子动力学模拟程序的组成，并实现运用该程序进行切削 仿真实验的目的，具体内容如下： 1 开发可应用于多种环境的分子动力学模拟软件。该软件可以使用各种不 同的势函数来描述不同的晶体结构；根据不同的精度要求及计算速度要求，使 用不同的有限差分方法进行计算；具有多种可供选择的结果输出：可以直观的 对模拟结果进行观察和分析。此外，该软件还应具有建模简单，方便操作，可 扩展性强等特点。 2 使用分子动力学模拟软件进行纳米切削的模拟实验，充分发挥模拟软件 自身的特点，通过分析模拟结果，对纳米切削机理进行探讨。 3 对分子动力学的一个发展方向多尺度模拟方法进行初步探讨。这是 一种目前还不完善的方法，但它的应用前景可见一斑，发展该方法对拓展分子 动力学方法的应用范围有重大意义。 哈尔滨工业大学工学硕士学位论文 第2 章分子动力学原理及软件实现 2 1 分子动力学仿真方法的基本原理 根据统计物理学原理可知，对于由大量粒子组成的物理系统，其宏观特性 是这些粒子运动状态的综合反映。因此，分子动力学的核心问题就是要计算所 有粒子的运动规律和轨迹特征。根据实际情况，计算中作如下假设： ( 1 ) 系统中粒子的运动都遵循经典牛顿力学规律； f 2 ) 粒子之间的相互作用满足叠加原理。 这两条假设忽略了量子效应和多体作用，使模型与真实物理系统存在一定 差异，但并不妨碍我们对一般性系统做出预测。 考虑具有个粒子的物理体系，其体系的h a m i l t o n 量可以用广义坐标的形 式表述如下口s 】： 日( p ，哼) = 足( p ) + u ( q ) ( 2 - 1 ) p = ( p l ，p ：，忽) q = ( g l ，吼，吼) 式中口粒子的广义坐标： p 粒子的广义动量； k f p ) 体系的总动能： ( ，f g ) 系统的总势能。 广义坐标叮可以有多种形式，视应用的方便而定。它既可以是笛卡尔坐标， 也可以是原子质心的笛卡尔坐标加原子的取向变量( 在需要把原子看成个刚 体的时候) 。 广义动量p 是一组共轭动量，通常动能部分k ( p ) 有如下形式 世( p ) 2 善 耐1 2 + p ：，+ 屯) ( 2 - 2 ) 势能部分v ( q ) 则包含了所有原予内和原子间相互作用的全部信息 矿( 口) 2 莩。) + 莩丢v 2 h 。j + z 。z 。州z ，v 3 ( 川+ ( 2 - 3 ) 当动能和势能都表达清楚后，就可以建立该体系的运动方程 一篓型鳖呈堡篓鐾垒一 在笛卡尔空间进行数值积分求解。体系的温度与各原子的平均速度有关 3 k 。r = 帆v j 行 ( 2 ，5 ) 悟l 式中k 。- b o l t 2 m a n n 常数； r 系统温度： m 第f 个原子的质量； v 第i 个原子的速度。 分子动力学模拟计算的基本过程如图2 - 1 所示。 图2 - 1 分子动力学基本原理流程图 f i g 2 ib a s i cm d s i m u i a t i o np r i n c i p l ef l o w ；h a r t 从原子的位置、链接方式、各种势能函数、原子的速度计算出体系的总能 量。然后计算各个原子在该力场中的势能梯度。每个原子在分子力场中所受的 力计算出来以后，按照牛顿第二定律，就可以计算原予运动行为。这是一个不 断迭代的过程。步长大约在i 飞秒( 1 0 。5 秒) 左右。足够次数的循环迭代将完成体 系运动方程的积分过程，而得到一个多体问题的解，和在相空间( 即原子的位置 及动量坐标的空间) 中的运动轨迹。这样，原予体系中每个原子的彳亍为都记录在 案，描述得清清楚楚。 嚆 珂1 1 一 一受一卜一 一量一一所一一一：一 一 硒面赢一 一 一系一个一瞅 一器量蓍 哈尔演工业大学工学硕士学位论文 2 2 分子动力学仿真的主要技术 2 2 1 势函数 由2 1 节的原理描述可知，要得到体系内原子的运动轨迹，就要对式( 2 4 ) 进 行求解。原子的动能可以根据其自身质量和速度求出，但原子的势能却不能由 简便的形式得到。原子间的相互作用十分复杂，目前还没有任何公式能够对其 进行准确的描述，只能根据一些表面现象和实验数据对其进行模拟，这就是目 前采用的势函数。根据势函数的计算才能得到原子的受力情况，从而计算出原 子的运动状态，可见势函数的准确性对计算结果的影响极大。 然而，势函数的研究一直发展得很缓慢，从一定程度上也制约了分子动力 学在实际研究中的应用。 原子间势函数的发展经历了从对势理论到多体势理论的过程。对势认为原 子间的相互作用是两两之间的作用，与其它原子的位置无关，而实际上，在多 原子体系中，一个原子的位置不同，将影响一定空间范围内的电子云分布，从 而影响其它原子之间的有效作用，故多体势理论能够更准确的描述多原子体系 的相互作用。 ( 1 ) 对势在分子动力学模拟初期，人们经常采用的是对势。应用对势的首 次模拟是a l d e r 和w a i n w r i g h t 在1 9 5 7 年的分子动力学模拟中采用的间断势【1 0 1 。 1 9 6 4 年r a h m a n 应用了非间断的对势对氩元素进行了研究0 3 9 】，他和s t i l l i n g e r 在 1 9 7 1 年也首次模拟了液体h 2 0 分子【4 们，并对分子动力学方法做出了许多重要贡 献。比较常用的对势主要有l e n n a r d j o n e s 势 4 l 】和m o r s e 势【4 孙。其势能及作用力如 图2 - 2 所示。 l e n n a r d j o n e s 势的具体表达为： ( _ ) = 4 占i ( o 盯) 一”一( o 盯) ”l ( 2 6 ) 式中盯原予与妒h ) 中零点势的距离； s 最小势能处的能量； m ，”调整系数，- 越t ( m ，n ) 取值为( 1 2 ，6 ) ，( 1 0 ，5 ) ，( 8 ，4 ) 。 l e r m a r d - j o n e s 势是为描述惰性气体分子之间相互作用而建立的，也有人用 它来描述铬、钼、钨等体心立方过渡金属。 啥尔滨工业大学工学硕士学位论文 i i i i l 一铖舭 i - 只d x a s i 厂7 j 一一 i 圈2 - 2l e r m a r d - j o n e s 势和m o r s e 势 f i g 2 - 2p a i r - w i s ep o t e n t i a la n di n t e r a t o m i cf o r c eo f l e n n a r d j o n e sa n dm o r s e m o r s e 势函数的形式为： ( 7 ；) = a p 一2 “( 一) 一2 e 一吐一i ( 2 7 ) lj 式中彳结合能； 口势能曲线梯度系数； 弗原予间作用力为零时的原子间距。 m o r s e 势函数经常用来模拟c u 等金属圆体，以及构造多体势的对势部分。 ( 2 ) 多体势多体势于2 0 世纪8 0 年代初期开始出现。多体势的描述方法主 要有：d a w 和b a s k e s 在1 9 8 3 年首次提出的嵌入原子法( e a m ) f 4 3 】；f i n n i s 和s i n c l a i r 在1 9 8 4 年提出的f s 势函数【4 4 l ；t e r s o f f ：p1 9 8 6 年提出的t e r s o 罐级势【4 5 l 等。 e a m 势函数的基本思想是把晶体的总势能分成两部分：一部分位于晶格点 阵上的原予核之间的相互作用对势，另一部分是原子核镶嵌在电子云背景中的 嵌入能，它代表多体式作用。构成e a m 势的对势与镶嵌势的函数形式都是根据 经验选取。 在嵌入原子法中，系统的总势能表示为： u = f ( 尼) + 丸( o ) ( 2 8 ) fi t f 式中，嵌入能； 一其它原子的核外电子在第i 个原子处产生的电子云密度之和。 哈尔滨工业大学工学硕士学位论文 式( 2 8 ) 中的第二项是对势项，根据需要可以取不同的形式。原理上f ( p ) 能 写成多项式形式其值可以通过拟合金属的各种物理性质的实验数值得到，也 可由规则金属晶格膨胀或压缩的状态方程求得。e a m 的最大特点便是根据嵌入 函数f 与电子云密度p 的关系，只要给定原子排列构型，就可以通过计算电子 云密度来确定系统的总能量，因此适于研究非对称体系。e a m 势很好的描述了 金属原子之间的作用，是描述金属体系最常用的一种势函数。本文将引入e a m 势函数，用于c u 年l l a l 等金属的模拟切削研究。 此外，基= y e a m 的势函数还有很多种 4 6 , 4 7 】，这些多体势也大都用于金属的 微观模拟中。 f s 势函数是根据金属能带的紧束缚理论发展出来的一种在数学上等同于 e a m 的势函数。它给出了多体相互作用势的函数形式，即将嵌入能函数设为平 方根形式。a c k l a n d 等在此基础上通过拟合金属的弹性常数、点阵常数、空位形 成能、聚合能及压强体积关系给出t c u 、a 1 、n i 、a g 的多体势函数h 引。其中( 2 8 ) 式中的多体项分别为 乃( 勺) = 4 ( 置o ) 3 h ( 墨一o ) ( 2 9 ) k = l 办( _ ) = q ( o ) 3 目( 咯一，；，) ( 2 1 0 ) k = l 当z 0 时，( 上) = 0 ：当j 足 r 2 r 6 。它们的值随金属物质不同而有所不同。 t e r s o 磷函数是一种用来描述共价键相互作用的原子间相互作用势，它主 要考虑了共价键之间的角度对原子能量的影响。具体形式如下： 力= 兀( o ) 厶( o ) + 厶( o ) ( 2 - 1 1 ) 式中 五( o ) = 以c x p ( 一毛o ) 厶( 芬) = 一b oe x p ( - , u q r e ) 一：筌霎鳖笙墼型塑鎏；。一 l ，；， r + c 。s 万( 勺一凡) ( 毛一毛) ，凡 。 毛 = 蜀( 1 + 缈舒) 。“1 乞= 厶( ) g ( ) l t ，j g ( 让，q 2 d ，2 一c 竹+ ( 曩一s ) 2 凡= ( 十一) 2 ，一= ( “+ 一) ，2 呜= ( 4 4 ) ，岛= ( e 哆j r , j = r r ，j ，岛= ( 置s ，) 其中：f ，k 为系统中每个原子的下标，l ，为原子i 与原子之间的键长。谚， 即表示原子扩间的能量关系，厶( _ ) ， ( o ) 分别表示排斥和吸引函数，厶( o ) 为 截断函数。为原子l 与原子i k 的键角，z ，为修正系数，是反映异极键强度的 独立参数，为考虑到原子、分子更加复杂几何形状分布的修正系数。式( 2 一1 1 ) 中各参数的具体取值见参考文献 4 5 1 。 t e r s o f f 势函数能够很好的描述碳族原予如c 、s i 、g e 等晶体的结构，在模拟 切削中常用来描述金刚石结构，本文将在模拟软件中实现t e r s o 瞄夸函数。 2 2 2 有限差分法 在得到势函数的具体形式及参数后，对其求导即可得到原子的受力情况如 式f 2 1 2 ) 所示 f ：= ( 一掣) ii j q 这样牛顿原子的运动可描述为 争= 去善砘，驴2 磊刍，嗨) ( 2 一1 2 ) ( 2 1 3 ) 哈尔滨工业大学工学硕士学位论文 原子的运动方程是一组微分方程，求解该运动方程一般采用有限差分法，以 一定的时间步长沿时间轴对方程逐步积分。但由于计算分子间作用力所占的工 作量极大，所以凡是要求在一个时间步中重复计算分子间作用力的算法是不合 适的，例如著名的r u n g - k u t t a 法在每一时间步中要四次计算分子间作用力，因此 在分子动力学模拟中不宣采用。 适用于分子动力学的模拟算法有e u l e r 法、g e a r 法、b e e m a n 法、v e r l e t 法、 v e l o c i t y v e r l e t 法和r a h m a n 法等。其中常用的有高阶g e a r 法和v e r l e t 法。 g e a r 于1 9 7 1 年提出了基于预测校正积分方法的g e a r 算法【4 9 1 。这种方法的 基本思想是用f 时刻位置和其初始的h 阶导数信息来预测在h f 时刻位置和其 初始的撑阶导数，然后计算在预测位置上的力( 加速度) 。然而，我们可以发现这 些加速度并不等于预测值。因此，调楚对加速度的预测以接近真实值。在预测 加速度和观测到的加速度不一致的基础上，再进一步尝试改进对位置及其”一1 阶导数的估计，这就是预测校正算法的校正部分。 考虑在f + f 时刻某原子坐标的t a y t o r ) 好式： ，( r + r ) = r ( ，) + 加百o r + 酉a t 2 萨0 2 r + 可a t 3 可0 3 r + ( 2 - 1 4 ) 使用记号 ( r ) = r ( r ) ，( r ) s f 象, x 2 ( r ) * 百a t 2 萨0 2 r ，弓s 等窘 可以通过屯( t + & ) 写出对x o ( t + f ) 的如下预测： x o ( t + a t ) = x o ( t ) + x ( t ) + x 2 ( t ) + x 3 ( r ) 搿篙三州。+ 2 x 2 ( ( t ) + 3 x 3 ( ( t 2 x 2t ) + 3 x 3 ，j ( 2 - 1 5 ) 屯( f + f ) = ( ，) v 1 叫 为( t + a t ) = 屯( f ) 这样，有了x o ( t + a t ) ，就可以计算在预测位置的力，计算z ：( t + a t ) 的校正值。 记“与矿“的差为蝇t 现在可以通过而写出对而的如下的预测： 。= 群砌“+ c a x 2 ( 2 1 6 ) 式中e 是常数，对某给定阶的算法其值固定，例如对于五阶预测校正算法f 即 至) ，e 的值为 哈尔滨工业大学工学硕士学位论文 c 0 = 嚣，q = 百3 ，c 2 = 1 ，g = 砭1 ，c 4 = 西l g e a r 算法具有较高的精度，但是每步要计算两次作用力，增加了存储量负 担，因此其效率并不是很高。现在应用最广的是v e r l e t 算法，它的优点在于存储 量小且稳定性好，易于程序实现。v e