分子模拟（精选十篇）
分子模拟 篇1
为解决这个难题,让学生可以轻松地理解高分子物理中的现象,比如高分子结晶,我们创造性地将计算机分子模拟方法引入到高分子物理的课堂教学过程中,通过分子模拟方法直接从分子原子尺度展示高分子的微观聚集态结构的形成过程,从而使学生更容易理解高分子物理中的相关现象和理论。下面我们以讲授高分子结晶为例,具体地介绍分子模拟方法在教学过程中的应用。
高分子结晶过程分为两个阶段: 结晶成核和晶体生长阶段。其中成核阶段所形成的晶核尺寸较小,因此成核过程难以在实验中被直接观察到。学生在学习过程中也难以对成核的微观过程产生一个直观的印象。为了帮助学生理解成核过程,我们利用分子模拟的方法对聚合物的成核过程进行了模拟,直接向学生们展示微观成核过程,获得了不错的效果。下面我们列举两个列子,具体地说明分子模拟在教学过程中的应用。
1 温度对结晶成核的影响
根据传统的结晶成核理论［2］,温度越低,聚合物体系的过冷度越大,成核自由能位垒就越小,因此其成核速度就越快;反之,温度越高,成核速度就越慢,体系就越难以成核。在教学过程中,我们利用动力学蒙特卡洛分子模拟方法从分子尺度向学生展示了不同温度下的高分子成核过程。图1 显示了较低温度下高分子的成核过程。从图中可以看出,随着结晶时间的增加,晶粒逐渐变大。此外,晶粒在增大的过程中,也在不断地发生波动,即某个高分子键参与形成晶核之后,会融化掉,然后又再次参与结晶。而这种成核过程中的波动在实验中难以被观察到。通过分子模拟方法可以直接观察到这个成核的波动过程,从而便于学生理解。图2 显示了较高温度下晶核的演变过程。随着时间的增加,在局部空间先形成一个很小的晶核,然后晶核变大,但是当它的尺寸达到某个值时,晶核开始发生熔融,逐渐变小,最后消失了。也就是说,在较高的温度下,晶核无法不断增大,而是在热涨落的作用下先形成而后又消失。这样通过分子模拟,我们成功地向学生展示了实验中难以观察到的高分子的成核过程。
2 拉伸诱导结晶
拉伸过程也会显著地改变高分子的结晶行为。拉伸会使高分子链发生取向,减低构象熵,导致熔点上升,过冷度变大,从而提高高分子的结晶能力,使其发生应变诱导结晶［3－5］。但是如何在课堂上直接向学生展示拉伸诱导结晶的过程,是教学的难点。为了解决这个问题,我们同样使用分子模拟方法直接模拟拉伸诱导高分子结晶的过程,然后在课堂上向学生展示,使他们更容易理解拉伸对高分子结晶过程的影响。图3 显示了拉伸过程中高分子晶体的演变过程。从图中可以看出,在应变较小时( 75% ) ,体系中几乎没有形成晶体。随着应变越来越大,晶体开始出现且越来越大,而且形成的晶体倾向于沿着拉伸方法发生取向。这样我们就成功地在课堂上向学生演示了拉伸诱导高分子结晶的微观过程。
另外,通过分子模拟方法,我们还可以发现一些课本上没有展示的信息。比如,图4 显示了两根高分子链在不同应变下的构象情况。在应变较小时,可以看出高分子链中存在一些折叠链构象的晶体,而当应变变大时,含有折叠链构象的晶体几乎全部转化为伸直链构象的晶体。这说明,不同应变下的高分子晶体会展现出不同的形貌: 低应变时以折叠链构象晶体为主,而高应变时以伸直链构象晶体为主。
3 结论
计算机分子模拟可以从分子尺度直接展示一些高分子物理现象的微观演变过程。我们创造性地将分子模拟的结果引入教学过程中,与理论及实验现象密切结合,给学生直接展示相关高分子物理现象的微观过程,使神秘复杂的物理现象变得简单而直观,从而易于学生理解。这是高分子物理教学过程中的新的探索。
摘要：高分子物理是研究高分子材料微观结构、运动模式、物理性能及它们相互之间关系的一门重要基础性科学。高分子材料微观结构和运动模式难以在实验上直接被观察到,较抽象,这使得学生难以理解。我们利用分子模拟方法直接从分子尺度展示高分子材料的微观聚集态结构的演变过程,使学生更容易理解高分子物理中的相关现象和理论。
分子模拟 篇2
铜配合物的分子几何结构模拟方法研究
对比和评价了用分子力学MM+法和半经验PM3法模拟23个铜配合物分子几何结构的结果.MM+法对四配位配合物的几何模拟结构与晶体结构差别大而不适于模拟这类分子,其可以模拟五配位三元配合物,预计这种配合物Cu-N键和Cu-O键的.偏差将分别高达(9.45±0.062)%和(18.2±0.110)%.PM3方法可以较好地模拟所考察的各种类型铜配合物,预计四配位配合物的Cu-N键和Cu-O键的偏差分别为(3.26±0.023)%和(2.77±0.023)%.而3+1+1型配合物的该偏差分别为(5.18±0.032)%和(4.31±0.037)%.
作 者：刘琼 李海平徐小燕 黄凤仙 陈希 王明召 Liu Qiong Li Haiping Xu Xiaoyan Huang Fengxian Chen Xi Wang Mingzhao 作者单位：北京师范大学放射性药物教育部重点实验室,北京师范大学化学学院,100875,北京刊 名：北京师范大学学报(自然科学版) ISTIC PKU英文刊名：JOURNAL OF BEIJING NORMAL UNIVERSITY(NATURAL SCIENCE)年，卷(期)：200743(5)分类号：O6关键词：铜配合物 几何结构模拟 MM+ PM3通信作者
分子模拟 篇3
关键词：GaussView软件；中学化学；分子结构；晶体结构
文章编号：1005–6629（2015）6–0090–04 中图分类号：G633.8 文献标识码：B
在中学化学教材中关于粒子微观结构的教学内容占据了较大的比例，如初中化学内容关于分子、原子、离子的概念，元素的概念等，再如高中化学内容关于物质结构、有机化合物以及晶体结构与性质等[1]。这类知识在中学阶段大致可以归类为分子结构和晶体结构两部分内容。在化学学科中，物质的性质由物质的几何结构和电子结构决定[2]，因此，学生只有知道分子和晶体的几何结构以及原子和一些简单分子的轨道图形才能更好地理解和掌握化学知识。
1 图形软件在中学化学教育教学中的应用现状
随着信息技术的全面发展，将图形软件作为一种信息资源辅助教学已经成为一种新的教育理念。目前，已经有许多图形软件应用于中学化学教学中，如GaussView、Chemoffice、Shockwave3D、Authorware等软件，其中Shockwave3D软件对教师的计算机水平要求较高[3]，在中学阶段难以推广；Chemoffice软件在平面分子和3D分子的设计以及化学实验装置的绘制方面具有较大的优势[4]，但是该软件不能绘制晶体结构；Authorware软件则在多媒体课件的制作方面有着强大的功能[5]，但对于分子结构及晶体结构的模拟不够精确。
由于分子结构和晶体结构作为中学化学的难点，同时也是学生今后学习有机化学和结构化学的基础。对于分子结构的教学，目前CAI技术、GaussView软件、Chemoffice软件均可以实现其空间结构的模拟，后两者对分子结构的模拟更加精确。比较而言，Chemoffice软件更适合绘制复杂的有机分子的结构以及化学实验装置。对于晶体结构的教学，尚没有成熟的软件使用。根据我们使用GaussView软件的经验，发现该软件在分子和晶体结构的设计方面具有很多优势，且该软件界面清晰、操作简单，较好地凸显了其辅助特性。
2 Gaussview软件介绍
GaussView软件是Gaussion98计算程序的辅助软件，该软件最大的特点是软件自身小且专业性强，其模版库可视化效果好，生成文件小，操作简单、快速（如转速快），便于学生观察感知，特别有利于学生空间思维的培养，从而彻底改变学生对抽象概念认知模糊的状况，使教学过程变得简单、清晰，是目前理论化学首选的分子设计软件[6]。
2.1 GaussView软件的界面介绍
GaussView软件主界面包括一大一小两个窗口，前面紫色的窗口为绘图窗口，见图1；模型库界面包括元素周期表、环状官能团、链状官能团三个界面，其可视化效果非常好，见图2。
2.2 GaussView软件的功能介绍
GaussView软件最大的特点是模型库清晰直观，可以直接通过添加模型库中的基团构型来绘制分子结构，操作起来简单快速。具体来说主要功能有几点，一是快速而准确地画出各种分子及其过渡态的空间结构，包括有效地呈现分子的同分异构体，可以通过键长、键角、二面角设置功能来直观地呈现分子的键长、键角和二面角等参数；二是设计晶体结构；三是显示分子轨道，即显示Gaussion98软件对分子结构的计算结果[7]。在绘制微观粒子空间构型的时候可以利用其剪贴功能任意添加、删除原子或者原子基团，可以任意角度的旋转粒子结构，也可以在X、Y、Z轴方向直观地查看粒子空间结构。比较其他化学图形软件，GaussView软件在分子结构和晶体结构的绘制方面有较大的优势。
3 GaussView软件的实际应用
3.1 GaussView软件对分子结构的绘制
对于分子结构的绘制，目前中学化学教学中使用较多的是Chemoffice软件。该软件虽然在平面分子和3D分子的设计以及化学实验装置的绘制方面具有一定的优势[8]，但对于分子结构的绘制，由于其组件软件Chem3D界面不清晰，所需基团不能从模型库直观找到，而是需要根据官能团名称来查找，加大绘制难度；再者，绘制的分子结构不能通过调整二面角来查看分子空间结构，分子结构的可控性较差[9]。
再如，使用GaussView软件可设置分子的二面角，通过旋转基团可查看不同角度的分子结构，图4中“5”“6”“7”分别为任意二面角的乙烷分子结构以及甲基基团的交叉和重叠（二面角纵向结构）。该软件还可以任意替换原子或基团绘制其他分子结构。下图5中“9”溴苯分子结构是通过将“8”中苯分子中的任意一个氢原子替换为溴得来。该软件还有自动调整分子功能，使其变得更加规范。同样用改变二面角和替换官能团方法可绘制丁烷和环己烷的两种不同结构[10]，如图6中“10～13”。
综上所述，GaussView软件在分子结构的绘制方面有如下优势：一是模型库界面清晰直观，可直接用来绘制分子结构；二是绘制方法多样，除了直接绘制还可以通过替换官能团来绘制分子结构；三是可以通过调整分子的二面角来呈现不同角度的空间结构。以上几点是Chemoffice软件所不具备的优势。
3.2 GaussView软件对晶体结构的绘制
目前，在晶体结构的绘制方面，虽然有许多软件在使用，比如在结构化学中通常使用CarIne Crystqallography、Materials Studio、diamond等软件来绘制晶体结构，但以上几款软件对于中学化学中简单晶体结构的绘制来说参数设置复杂，安装复杂，并且大部分软件不支持中文，对于初学者来说使用困难。所以，在中学化学晶体结构的教学中，GaussView软件更适合中学教师的使用。以下使用GaussView软件来绘制晶体结构，具体操作如下：
建立晶体结构之前，首先查找晶体的空间群信息，空间群的标记有Hermann–Mauguin记法或空间群序号（Space group number），如面心立方的H-M标记为Fm-3m，空间群序号为225，打开链接http：//www.webelements.com/gold/crystal_structure. html可查看晶体晶胞参数，本文以面心立方中的铜金属晶体和氯化钠离子晶体为例，建立二者的晶体结构模型，查得
“Cu”晶胞参数
Spacegroup：Fm-3m（Spacegroupnumber：225）
Structure：ccp（cubicclose-packed）
a：361.49pm b：361.49pm c：361.49pm
α：90.000° β：90.000° γ：90.000°
“NaCl”晶胞参数：
Spacegroup：Fm-3m（Spacegroupnumber：225）
Structure：ccp（cubicclose-packed）
原子分数坐标
A：（0，0，0）、（1，1/2，1/2）、（1/2，0，1/2）、（1/2，1/2，0）
B：（1/2，0，0）、（0，1/2，0）、（0，0，1/2）、（1/2，1/2，1/2）
第一步：新建一个文档：File→New→Create molecule group；
Symmetry设置，lattice dimension（晶格大小）：three；最小晶格系统：cubic；空间群：225[Fm-3m]；Cell设置，celllengths、cellangles按照铜的晶胞参数设置；View设置（显示所有边界原子）；重复内容显示：normal；Content设置：添加铜原子（系统默认为5个）进行分数坐标设置；最后点击Combine生成晶体结构。下图7中“14”为铜金属晶体，“15～16”为不同角度氯化钠晶体，其中紫色为钠离子，绿色为氯离子。
从氯化钠晶体结构中可以直观地看到一个NaCl晶胞含有的Na+和Cl-的个数均为4个。以上是GaussView对晶体结构的绘制，而目前为止Chem3D对晶体结构的绘制只限于软件本身模型库中的少量子模型，对于种类繁多的晶体结构模型的绘制则没法实现。
4 结论
通过对比可知，GaussView软件在中学化学教学中有其独特优势，特别适合绘制微观粒子的空间结构模型，再加上该软件操作简便快捷、安装简单，对计算机水平要求较低，适合于中学阶段的教学使用；虽然GaussView软件在绘制晶体结构方面对结构化学专业知识要求较高，但由于该软件本身设计良好，绘制晶体结构时只需知道晶体的空间群、晶胞参数、分数坐标等信息，这些信息可通过专业网站查找，而其他软件对计算机水平要求比较高，难以在中学阶段使用。鉴于中学阶段涉及到的分子结构的教学比较简单，不涉及复杂的有机分子，且考虑到中学教师的计算机技能强弱的问题，就分子结构和晶体结构这两块教学来说，GaussView软件更适合在中学教学中使用，既便于教师快速操作使用，更便于学生感知。但该软件本身也有一些缺陷，比如目前是非免费使用软件，普及使用有待于进一步改进。
参考文献：
[1]中华人民共和国教育部制定.义务教育化学课程标准[S].北京：北京师范大学出版社，2011.
[2]中华人民共和国教育部制定.普通高中信息技术课程标准（实验）[S].北京：人民教育出版社，2003.
[3]王清华. Shockwave3D三维动画技术在化学教学中的尝试[J].化学教学，2004，（11）：38～40.
[4]王灵，胡乔生. ChemBioOffice2010在化学教学中的应用[J].软件导刊（教育技术），2010，（8）：87～89.
[5]孙鸿. Authorware在中学化学教学中的应用[J].株洲师范高等专科学校学报，2005，10（5）：89～90.
[6][7]张翼. Gaussian98和Gaussview在中学化学教学中的简单应用[J].化学教学，2005，（5）：47～49.
[8]刘芬芬.三种工具软件在化学教学中的应用探究[D].西安：陕西师范大学硕士学位论文，2013.
[9]陈兰美. Chem3D在有机化学教学中的应用[J].广东医学院学报，2009，27（6）：712～714.
一维纳米摩擦分子动力学模拟 篇4
世界上使用的能源大约有1/3～1/2消耗于摩擦, 如果能够尽力减少无用的摩擦消耗, 便可大量节省能源。另外, 机械产品的易损零件大部分是由于磨损超过限度而报废和更换的, 如果能控制和减少磨损, 则既减少设备维修次数和费用, 又能节省制造零件及其所需材料的费用[1,2]。作为研究相对运动的作用表面间的摩擦、润滑和磨损, 以及三者间相互关系的理论与应用的一门边缘学科的摩擦学自古至今都得到了相应的重视。摩擦学的研究得到了相应的摩擦学库仑定律和Amoton规律, 这些经验规律能够合理地解释生活生产中的摩擦现象, 对于指导人类日常生活和社会生产的进行起到了重要的推动作用;但是本世纪新兴的纳米技术和纳米材料研究的不断深入, 使得人们有能力制造出尺度越来越小, 而性能越来越完善的微型装置, 尤其在备受关注的生物、环境控制、医疗器械、航空航天、数字通讯、传感技术等前沿领域快速发展, 纳米机械的发展使得摩擦学问题再次成为科学研究热点。因为对于纳米尺度的物质系统, 基于连续介质假设的传统理论可能不再有效, 纳米尺度下的摩擦学系统由于纳米材料自身的特性将表现出不同的摩擦学规律, 而“纳米摩擦学”即是为解决微纳米系统中摩擦的规律和控制而逐步形成的。
近年来纳米摩擦学的发展较为迅速, 主要基于两方面原因:一方面是纳米摩擦测试技术的进步使得对纳米摩擦现象的实验测试成为可能。在实验室中可以通过原子力显微镜[3], 摩擦力显微镜[4], 石英晶体微天平[5,6]等对纳米摩擦过程进行观测, 从而能够更好地理解摩擦的起源问题;另一方面是计算机计算能力的快速增长和计算软件的不断更新。计算能力的快速增长使得复杂的摩擦现象可以通过分子动力学进行相对真实的模拟再现, 从而为深层次地理解摩擦现象提供了新的研究手段[7,8,9]。
本文将经典FK模型和能量耗散机制相结合开展类Xe原子在类Ag表面滑动摩擦的分子动力学模拟[10]。为了得到在一维纳米摩擦的一般规律, 而不强调求解某种材料的摩擦特性, 所以在模型中设定相应的参数调节来得到不同的摩擦环境。将Xe与Ag的相互作用势能改变为一种组合周期势能进行描述。模拟过程中原子在受到不同大小附加拉力的作用下, 在不同温度和不同匹配度的衬底表面进行滑移。模拟过程除了设置了界面势能大小, 匹配度等参数 (图1) , 还设置复合周期势能中正弦势能振幅比例UA, 势能周期比例UP两个重要的参数。振幅比例UA, 势能周期比例UP的调节可以很好地体现出一般衬底晶格中含有的两种原子之间的相对距离和两原子与滑移原子作用相对强弱对纳米摩擦过程的调制作用。
本文的第二部分描述分子动力学的模型和模拟过程中的参数设置。模拟的相关结果和讨论在文章的第三部分给出。第四部分对文章做出了总结。
2、模拟模型
本文模拟无限长的一维原子链在调节后的无缺陷原子表面的滑动过程来获得纳米摩擦学匹配度对摩擦的调节作用。滑移原子之间的相互作用采用Lennard-Jones势[11,12,13], 如下式:
其中, (i, j为标记不同的原子) ;0为T=0时的原子间的平衡距离;ε为势能阱深度。计算过程中截断到第三近邻原子, 晶格常数取α=0.99720[14]。
附加表面原子与衬底原子之间的相互作用以U (图1) 来表示, 由两个正弦势能复合而成如下式:
U (x) 是一个关于坐标的周期势函数, 周期为b, b=b1+b2;u1和u2分别为正弦周期势能的振幅的一半;b1和b2分别为两正弦周期势能的周期;定义势能的周期比参数, 即:b1= (1-UP) b, b2=UPb;势能振幅比为, 即:u1=u0, u2=UAu0 (其中u0为势能设置参数) 。
附加原子还受到一个复合高斯分布的随机力i, 随机力通过涨落-耗散机制与电子摩擦系数ηe相关[15], 且满足下式:
其中, 尖括号表示对时间求平均值;σ (t) 表示delta函数;k为Boltzmann常数;T为系统温度。根据涨落-耗散方程选取随机力i可以使系统温度保持在设定的温度。
在模拟过程中, 对所有的附加原子附加以恒定的拉力F, 经长时间驰豫系统进入稳定状态后求得质心速度, 由公式确定摩擦系数[14]。
附加的原子应满足如下的Lagevin方程[16,17,18]:
模拟采用对分子运动的数值积分得到不同状态下的摩擦过程。在模拟计算中, 为减少计算量, 提高效率, 对计算公式进行无量纲化处理。选用0, ε和m作为基本单位, 将其他单位进行换算, 例如时间单位变为。另外, 为了减少模拟尺寸的影响, 在模拟过程中采用了周期性边界条件, 所以在有限长度内附加原子数Nα和衬底原子数Nb之间地有如下关系:
其中, 比例常数C称为附加原子与衬底原子的匹配度。
为了得到一个稳定状态, 首先在附加拉力F=0的情况下, 时间步长△t=0.01t0, 对初始条件下的Lagevin方程进行模拟, 模拟过程进行30万步。模拟结果系统达到一个没有人为干扰的自然稳定的初始状态。然后对附加原子链中所有的原子附加拉力F进行模拟, 直到系统进入另一个稳定状态。在平衡状态下求解质心速度, 从而求出系统的摩擦系数。整个模拟过程设定在一定温度T下, 通过调节两势能周期比UP和振幅比UA, 附加拉力F等参数得到不同的模拟环境, 了解不同情况下的匹配度C对纳米摩擦的影响。
3、结果与讨论
首先, 给出在UA=UP=1极限情况下的摩擦模拟结果, 进行程序的正确性验证;然后计算不同势能周期比例和振幅比例随着匹配度变化纳米摩擦的变化。
3.1 一种极限情况UA=UP=1 (普通的正弦势能形式)
由前文可知, 当UA=UP=1复合准周期势能形式将变为一种简单的正弦形式, 也就是经典纳米摩擦分子动力学模型中所用的势能形式。
图2给出了附加拉力与质心速度关系模拟结果与文献[14]所得到的结果 (c, d为本文得到的模拟结果) 。通过对比发现本文的计算方法得到的摩擦系数会比文献[14]的要稍微大一些 (最大误差小于10%) , 但两种计算方法对摩擦系数的预测的趋势都是一致的。之所以出现误差的主要原因是在求解Lagevin方程时采用了不同的计算方法。但是因为本文的重点在于找出在复合势能形势下滑动摩擦所表现的规律, 所以两种计算结果的不同不会影响到摩擦规律的发现。
3.2 一般情况
为了得到纳米摩擦在一般情况下的规律性, 模拟过程中设置了相关可调参数, 如表1所示。
在KBT=0.1ε, 势能u0 (0.01ε～0.04ε) , UP (0.03～0.09) , UA (0.03～0.09) 情况下, 模拟C在 (0.03～0.09) 之间对摩擦系数的影响。图3给出了不同势能u0大小, 不同UP, UA情况下得到的摩擦系数的倒数随C的变化曲线。
由图可以看出, 最显著的特点是在不同摩擦条件下的摩擦系数倒数曲线随C的变化在C=0.5处出现了汇聚且出现最小值, 而在两侧则离散增加。当C=0.5时出现摩擦系数倒数最小值, 也就是出现了摩擦系数最大值, 即在完全匹配的情况下出现摩擦最大值。完全匹配处之所以出现摩擦最大值的原因在于此时所有的附加原子步伐统一, 能够同时受到衬底原子的拉力而加速移动, 又能够同时受到阻力减速, 相互之间无法将能量转化为势能保存, 而只能以热量的形式进行耗散。曲线的聚合同时说明对于纳米摩擦无论势能形式如何组合变化, 势能大小的如何变化, 在对摩擦的调制作用中界面晶格匹配度C对整个摩擦过程的调制作用最为明显, 在完全匹配的情况下摩擦最大。
由图可以发现另一个特点, 摩擦系数在C为0.5两侧随C的增大或者减少而迅速减小, 但在C为0.8时出现了转折点。此转折点在不同的条件下, 对摩擦的调制作用表现出不同程度的影响, 说明此转折点为复合势能情况下产生的影响, 而这种调节作用与势能周期比例UP和势能振幅比例UA之间存在非线性关系。转折点的出现原因为附加原子在不同的衬底复合势能阱中进行运动时受到衬底原子的拉扯作用产生了强烈的声子振动所带来的影响, 这种拉扯作用主要受到势能振幅比例与势能周期之间的比例关系的影响。当处于衬底势能最高点时, 附加原子将出现失稳并自动跳跃到下一势能最低点, 然后在此平衡位置剧烈震荡和激发声子, 从而使能量不可逆地以声子的形式耗散掉。当附加原子在势能阱较深, 而周期较小的势能阱中运动时受到的衬底原子拉力由最大值迅速变成了反向阻力最大值的过程中, 附加原子无法将拉力所带来的能量完全转化为原子之间的势能保存, 而只能以晶格振荡的方式转化为热的形式放出, 从而导致系统拉力作功成为摩擦产生热量过程, 进而实现能量以晶格振动的形式耗散, 摩擦系数增大。
模拟结果表明, 无论如何改变势能周期和势能大小比例组合势能, 均可以发现摩擦系数受到了匹配度的调制作用最大, 说明摩擦的主要原因来源于衬底与附加原子之间的晶格不匹配造成滑动过程中能量的耗散。
4、总结
本文在不同势能组合形式下, 不同晶格匹配度情况开展了纳米摩擦分子动力学模拟, 得到了一维纳米摩擦的一般性规律。模拟过程中主要设定了势能周期的分配, 势能大小的分配等几个重要参数, 通过固定不同参数模拟了匹配度对摩擦过程的调制作用。通过模拟发现在不同的摩擦环境中, 匹配度对纳米摩擦的影响出现了两个极值点, 一个为匹配度0.5, 此时纳米摩擦的主要调制作用由匹配度来决定, 其他因素对摩擦的调制作用不明显;另一个为匹配度0.8, 此时纳米摩擦的主要调节作用为势能振幅比例和势能周期比例调制, 这种调制作用实现对摩擦的控制作用, 为多原子的表面膜的摩擦系数控制提供了参考依据。通过控制多原子薄膜中不同原子种类和它们之间的晶格常数来控制薄膜的摩擦系数, 从而实现对摩擦磨损的合理控制。
摘要：本文在建立一维纳米摩擦模型的基础上, 对纳米摩擦现象进行了分子动力学模拟研究。研究结果表明, 表面接触界面之间的晶格匹配度对纳米量级的摩擦起到较大的调制作用, 而这种调制作用随着晶格匹配度的变化而产生调制作用的大小也不同。分析表明在完全匹配的情况下摩擦最大, 而其他因素例如界面接触势垒的大小、形式等与匹配度相比对摩擦的影响较小。
分子模拟 篇5
密度泛函与分子模拟计算介孔孔径分布比较
用巨正则系综 Monte Carlo模拟 (GCMC)方法和密度泛函理论( DFT)结合统计积分方程( SIE)计算了介孔材料的孔径分布.为比较这两种方法 ,以 77 K氮气在介孔活性碳微球中的吸附数据为依据 ,求出其孔径分布.在 GCMC模拟和 DFT计算中 ,流体分子模型化为单点的 Lerrnard-Jones球 ;流体分子与吸附剂材料之间的`作用采用平均场理论中的 10-4-3模型.在 DFT方法中 ,自由能采用 Tarazona 提出的加权近似密度泛函方法 (weighted density approximation,WDA)求解.结果表明 ,对于孔径大于 1.125 nm的介孔材料 ,GCMC和 DFT两种方法都可以用来研究介孔材料的孔径分布 ;对于小于 1.125 nm的介孔材料 ,不能用 DFT方法计算孔径分布( DFT方法本身的近似产生了误差) ,只能用分子模拟方法.
作 者：邵晓红 张现仁 汪文川 作者单位：北京化工大学化学工程学院,北京,100029刊 名：物理化学学报 ISTIC SCI PKU英文刊名：ACTA PHYSICO-CHIMICA SINICA年，卷(期)：19(6)分类号：O647.3 O641.2关键词：巨正则系综 Monte Carlo方法,密度泛函理论,孔径分布,吸附
畜禽保种群分子濒危标准模拟初探 篇6
世界自然保护联盟 (IUCN) 、联合国粮食及农业组织 (FAO) 、欧洲联盟 (EU) 、欧洲畜牧生产委员会 (EAAP) 等政府、非政府机构及各国学者在濒危评价标准方面做了大量工作[4]。然而从整体上看, 不同组织、不同国家的濒危评价标准存在差异, 标准不一, 具有很大的主观性, 这样就很难对物种的濒危状态做判定, 具体实施起来存在随意性缺陷。
遗传多样性分析是遗传资源濒危状态评价的一个新途径, 马月辉等[1]采用群体遗传学的原理, 建立了畜禽遗传资源受威胁的评价方法, 然而系谱记录原始数据的缺乏使得近交系数的使用受到限制, 另外近交系数只反映基因纯合位点, 不能有效反映群体基因频率的动态变化。
种群是否发生遗传变化很难在一两代或几代间表现出来, 利用计算机及试验动物做一些长期保种的模拟试验, 可以提高保种的预见性[5]。分子生物学的发展给畜禽保种工作提供了新方法, 研究拟根据分子标记特点, 结合分子遗传多样性评价指标, 运用群体遗传学知识, 利用计算机模拟中期、长期畜禽遗传规律变化, 研究不同濒危程度对畜禽遗传多样性的影响, 为建立畜禽濒危资源预警机制及合理规划保种方案提供依据。
1 材料与方法
1.1 模拟参数设置
模拟参数设置参照参考文献[1, 6], 以表1[1]中遗传资源受威胁程度为参考, 根据不同畜禽物种世代间隔, 分别计算不同畜禽保种群有效群体大小 (Ne) 。根据公式分别计算随机交配、随机留种及随机交配各家系等量留种模式下保种群体的公母数。式中:Nm、Nf分别表示保种群体的公母数。当随机留种为Nm∶Nf=1∶5时, 则Nm=0.3Ne, Nf=1.5Ne;各家系等量留种为Nm∶Nf=1∶5时, 则Nm=0.2Ne, Nf=Ne。
1.2 模拟试验的实现
模拟试验参照参考文献[8-9], 根据设置参数采用Fortran90自行编写的程序完成。根据前期试验结果选择期望杂合度 (He) 、观测等位基因数 (Ao) 为分子度量指标[6,9]。
2 结果与分析
按照不同保种群畜禽世代间隔计算其不同濒危等级下繁育公母数及Ne值, 统计各濒危等级条件下基础群及100年对应世代数时的He、Ao, 并计算这2个指标的下降幅度, 见表1。
由表1可知, 基于100年近交系数分子濒危标准, 不同畜禽保种群随机交配、随机留种模式及随机交配、各家系等量留种模式下, 对应严重威胁、威胁、最低威胁、潜在威胁及安全5个等级临界点。实际应用时, 可以利用微卫星标记检测的分子数据模拟100年, 分别计算基础群和100年对应畜禽世代时的He、Ao, 计算其下降幅度, 与表1中的标准相对照, 以评估畜禽保种群濒危状态。
3 讨论
等级的评价标准应兼顾科学性和实用性, 这给标准的制订带来了极大困难。科学性要求这类标准客观准确和精细, 尽可能使用定量而不是定性的依据;而实用性则强调标准的简单性和实用性。两方面的要求有时是相互矛盾的, 因此标准的制订常处于一种两难的境地[2]。
J.D.Ballou等[10]指出保种的两个主要目标:一是群体规模, 二是遗传多样性。鉴于当前难以找到评价畜禽保种群濒危程度理想的标准, 因此在实际应用中, 可以结合畜禽保种群群体规模及遗传多样性进行综合考虑。
%
注:3~4列, 每行中上行数据为 (R∶R) 随机交配、随机留种模式;下行数据为 (R∶F) 随机交配、各家系等量留种模式;a为He0n↓, b为Ao0n↓。
摘要：为了制订畜禽中期、长期保种方案, 试验基于F-100濒危标准等级初步建立了分子濒危标准, 研究不同濒危程度对畜禽遗传多样性规律变化的影响。结果表明:F-100濒危标准等级可行。
关键词：畜禽保种群,分子濒危标准,模拟研究
参考文献
[1]马月辉, 吴常信.畜禽遗传资源受威胁程度评价[J].家畜生态, 2001, 22 (2) :8-13.
[2]王智, 蒋明康, 秦卫华.中国遗传资源优先保护等级评价标准的构建[J].环境科学与技术, 2008, 31 (5) :49-53.
[3]GANDINI G C, OLLIVIER L, DANELL B, et al.Criteria to assess the degree of endangerment of livestock breeds in Europe[J].Livest Prod Sci, 2004, 91:173-182.
[4]IUCN.IUCN eed list of threatened animals[M].Gland:Switzerland, 2001.
[5]吴常信.畜禽遗传资源保存的理论与技术[J].家畜生态, 2001, 22 (1) :1.
[6]鲁云风, 王茜, 黄文波, 等.畜禽保种群体遗传多样性的模拟[J].中国农业科学, 2012, 45 (18) :3849-3858.
[7]吴常信.动物遗传学[M].北京:高等教育出版社, 2009.
[8]张勤.动物遗传育种中的计算方法[M].北京:科学出版社, 2007.
[9]LU Y F, LI H W, WU K L, et al.Dynamic change of genetic diversity in conserved populations with different initial genetic architectures[J].J Integrative Agriculture, 2013, 12 (7) :101-109.
牛津大学液晶分子取向模拟软件项目 篇7
液晶显示屏 (LCD) 在电视、手机和标牌等方面应用广泛。LCD装置性能的最优化在很大程度上取决于指向矢与光线的传播和偏振的互动方式。因此, 对指向矢配置进行精确和可靠的模拟对于LCD装置性能的最优化非常重要。
牛津大学的研究人员开发了液晶分子取向模拟软件。该软件使用关系程序在一维、二维或三维内来决定向列型液晶装置结构的指向矢配置。该软件可用来理解液晶装置的复杂结构, 研究3D的松弛过程, 研究指向矢浮动的动力学, 使指向矢分布具象化以及模拟指向矢结构的光学等。
钨合金力学性能表征分子动力学模拟 篇8
钨和钨合金具有高强度、高密度, 良好的抗腐蚀性、导电、导热性和热膨胀系数小等一系列优良物理、力学性能。在军事领域中, 如:卫星、飞行器、航空发动机、核潜艇、导弹和舰艇等装备的一些关键部件大量采用了钨及钨合金, 它们的机构设计和安全性能问题与所采用的材料物理化学和机械特性有着十分密切的关系。
国内外许多研究学者对钨及钨合金进行了大量的宏观实验研究, Zhou等[3]采用拉伸冲击实验装置对钨合金试件在不同应变率下的工况进行了实验研究, 且在对冲击拉伸后断裂试件扫描分析, 建立了含有应变率相关参数的应力-应变模型。刘海燕等[4]基于旋转盘式间接杆-杆型冲击拉伸装置, 对不同颗粒度细化钨合金材料在不同应变率工况进行了实验研究, 通过分析钨合金材料在弹性阶段存在非线性, 判断了在此阶段材料内部已经损伤、钨颗粒减小材料屈服极限提高以及细化钨合金材料具有应变率致脆现象。李娜等[5]通过实验研究了83W和89W材料各向异性在不同温度环境的动、静加载条件下, 应力-应变关系及失效问题。Pink等[6,7]研究了钨合金材料在不同温度下的拉伸性能, 并发现钨合金材料在一定温度下出现动态应变时效现象和负应变率敏感性问题。王世良等[8]研究了钨合金的合成以及相应的物理问题, 谷怀鹏等[9]研究了镍基单晶合金物理特性问题。王建伟等[10]采用分子动力学方法研究bcc-Fe不同空位浓度对中子辐照损伤的影响。以上这些专家与学者分别对钨合金材料的实验研究、相关合金物理特性和金属材料的分子动力学模拟等方面均取得了一定的研究成果, 并促进钨及其合金材料的进一步发展。
然而, 利用扫描电子显微镜对钨合金材料和钨合金结构破坏观察研究分析, 钨合金材料主要由钨颗粒构成, 而在不同加载工况下, 钨颗粒的破坏是造成钨合金整体失效的主要原因。钨颗粒的强度决定了钨合金的物理力学性能, 但是由于钨的高强度、高硬度, 在宏观试验中很难测得其实际的物理性能的常规参数。针对这一问题, 采用分子动力学方法研究钨颗粒的材料力学性能成为了一种有效的手段, 钨颗粒尺寸对宏观材料尺寸相对无限小, 而对微纳米尺寸又可认为是无限大, 分子动力学方法它可以根据原子间的相互作用计算整个原子系统的性能, 分析材料的纳观变形机理, 再现宏观力学行为, 同时又提供大量的微观信息。
分子动力学模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法, 是指对于原子核和电子组成的多体系统, 应用牛顿力学运动方程, 对原子核运动过程进行模拟, 记录运动过程中的轨迹, 最终利用统计力学方法计算一些相关物理量, 得出物质结构特性的一种理论计算方法。该模拟手段已成功应用于研究晶体熔化、晶体冲击、晶格畸变等方面[11], 并且随着计算方法的发展以及计算机性能的提高, 将被广泛应用于生物、化学、医药等领域, 有着广阔的应用前景。采用分子动力学方法, 对钨合金的拉伸塑性问题还未展开系统研究。
利用分子动力学模拟钨合金拉伸塑性力学行为, 对[100]晶向进行温度效应、尺度效应、应变率效应的模拟研究, 得到了钨合金微观结构变化图像和应力应变关系等相关塑性力学数据, 以及解决了在微纳米尺度材料实验中较难观测原子轨迹及提取相应数据测量的重要问题, 对于进一步了解微纳米尺度钨合金材料力学性能及其材料制备与应用提供了重要的依据。
1 钨合金的拉伸模型
1.1 嵌入原子势
嵌入式原子势 (EAM) 是Daw和Baskes[12]根据Hohen-kohn定理建立的势函数模型, 具有运算量小、运算效率高、适合金属间的作用等优点。
EAM原子势的总能量E可以表示为:
式中:φij表示距离为rij的原子i和原子j之间的对势能;Fi表示电子密度为ρi的原子i嵌入能;Mi (Pi) 为修正项;Pi为原子i处电子密度中高次项的贡献。
1.2 钨合金计算模型
图1为金属体心立方结构示意图, 晶格常数为a=b=c=0.3165nm[13]。图2是以体心立方结构为基础, 建立钨合金圆柱形超晶胞计算模型 (钨合金主体钨为体心立方结构) , 初始构型按理想晶格排列, x, y, z坐标轴分别对应[100], [010], [001]方向。圆柱半径为r=14a, 长为l=40a, 两端分别固定长度h=2a为夹具区域, 模型左端固定, 右端采用“velocity”加载方式进行拉伸, 共有49773个原子。为模拟有限长度的力学行为, x, y, z方向均取自由边界条件。
2 模拟方法
采用Ackland的EAM势函数[14]描述钨合金原子之间的相互作用, 它能准确的拟合出钨合金结合能、缺陷形成能、层错能等参数, 更能准确拟合钨合金的弹性性质。首先, 本模型采用Hoover-Nose控温系统, 在正则系宗 (Canonical Ensemble, NVT) , 即表示具有确定的粒子数 (N) 、体积 (V) 、温度 (T) 环境下, 对模型体系进行弛豫, 将模型控制在一定的温度使体系达到一个初步稳定的状态。再利用微正则系宗 (Microcanonical Ensemble, NVE) 。即表示具有确定的粒子数 (N) 、体积 (V) 、总能量 (E) 对体系进一步弛豫, 达到理想的加载环境。在此基础上, 采取与宏观实验相类似的加载方式, 对钨合金圆柱模型进行拉伸加载问题的分子动力学计算模拟。
3 钨合金体系模型弛豫
所谓模型弛豫, 就是体系达到一个稳定的状态, 以便获得准确的数据。图3为[100]晶向、截面半径r=14a、长度l=40a、温度300K下模型弛豫的动能与势能的变化情况, 由图3分析可知, 在Hoover-Nose控温系统下, 在体系进行弛豫初始时, 体系动能与势能上下波动幅度较大, 此时体系并有达到平衡状态。在经过30ps的弛豫以后, 体系动能在2000eV小幅波动变化, 而势能在-433000eV左右小幅度周期性波动, 体系基本上达到一个相对稳定的状态, 可以对其进行加载计算。
4 数值模拟与实验验证
4.1 钨合金塑性变形数值模拟研究
4.1.1 在拉伸过程中的体系结构变化
采用[100]晶向、截面半径r=14a、长度l=40a、初始温度300K。首先在NVT系宗下弛豫50ps, 再在NVE系宗下弛豫50ps。对模型右端加载区域施加0.25nm/ps沿x轴正向速度进行拉伸, 中间区域为自由运动, 每步输出一次热力学信息 (温度、压力、势能、总能量、x轴应力、x轴应变) 共运行100000步。图4为对钨合金试件拉伸加载作用下, 原子运动变化趋势 (基于OVITO软件[16]对计算数据进行图形绘制) 。
当0≤t<7.5ps时, 在加载速率初始阶段, 钨合金为完整晶体, 原子结构为规则的晶格排列, 拉伸变形为弹性变形。随着晶体拉伸持续进行, 原子之间产生了位错, 并伴随着滑移变形, 在模型表面出现轻微滑移台阶。
当7.5≤t<15ps时, 钨合金在拉力的作用下不断伸长, 发生了由弹性变形到塑性变形的转变。从图中可看出原子在 (110) 晶面、沿[111]晶向进行滑移, 随着拉伸的进行, 原子的运动越来越快, 开始出现部分金属键发生断裂, 当应变不断增加时, 位错不断增殖, 滑移持续发生, 中间段横截面发生改变并不断的减小, 并出现了颈缩现象。当拉伸进一步伸长时, 塑性变形主要集中在颈缩部位, 其他部分原子位置基本不变, 颈缩区域不断伸长, 导致滑移系的连续发生, 颈缩区域横截面不断减小, 颈缩现象更加明显, 大部分金属键已发生断裂。
当15≤t<28.5ps时, 随着拉伸的进一步进行, 颈缩区域原子明显减少, 原子相互作用力越来越小, 当金属键完全断裂, 晶体发生断裂, 表明纳米尺度下钨合金高应变率拉伸变形方式类似于宏观体心立方材料的拉伸变形。
4.1.2 钨合金高应变率拉伸过程变化曲线
计算原子i的应力分量采用如下公式:
式中:m表示原子质量;a, b表示x, y, z轴对称张量的六个分量;r表示两个原子之间的距离;F表示两个原子之间的作用力。第一项表示原子i的动能, 第二项是原子i周围Np近邻原子数经过n次循环的对能, 第三项是原子i周围Nb的键能, 第四项是原子i周围Na键角对应力的作用, 第五项是原子i周围Nd面角对应力的作用, 第六项是原子i周围Ni非健项对应力的作用, 最后一项是原子i内部约束力Nf对应力的作用。
图5为不同工况拉伸加载下铝合金应力-应变曲线。从图5 (a) 可以看出, 分子动力学模拟与宏观实验得到的应力-应变曲线相同的是初始为线性上升, 应力与应变成正比例关系, 材料为弹性变形阶段, 区别是分子动力学模拟得到的弹性极限到最高点后发生突降现象, 并产生塑性流动。钨合金模型随着加载速率的增加, 屈服强度逐渐增大。钨合金在0.010, 0.050, 0.100nm/ps的高应变率加载速率拉伸过程中发生了二次屈服, 而在0.001nm/ps时钨合金模型没有产生二次屈服。
在分子动力学模拟中, 不同拉伸速率下的钨合金模型采用了相同的弛豫时间, 因此, 在较低的拉伸速率下, 弛豫时间充分, 钨合金可以再次达到稳定的平衡状态, 而在较高的拉伸速率下, 弛豫时间不足以使钨合金体系达到稳定的平衡状态。因此, 在相对较高的拉伸速率下, 钨合金体系会产生二次屈服, 而在较低的拉伸速率下没有发生二次屈服。
而MD模拟弹性极限到达最高点后发生突降现象, 是由于在拉伸过程中原子之间持续产生位错, 进一步发生了滑移, 当达到钨合金屈服应力峰值时, 原子之间滑移不断加剧, 引起原子之间大量金属键发生断裂, 故导致“突降”问题。
从对图5 (b) 应力-应变曲线分析可知, 钨合金模型在温度300K、长度l=14a、[100]晶向以0.1nm/ps速率加载的条件下, 随着模型横截面半径的增大, 在弹性阶段钨合金的屈服强度也越大。这说明在高应变率拉伸过程中, 横截面的越大, 在单位时间内产生的位错密度越大, 位错之间的交互作用也就越强, 增加了位错运动的阻力, 导致了屈服强度的提高。
由图5 (c) 可以看到, 在温度300K、横截面半径r=14a、[100]晶向以0.1nm/ps加载的条件下, 随着钨合金模型长度的增加, 钨合金在弹性阶段屈服强度在增大。从图中可以看出, 不同长度钨合金模型的应力-应变曲线斜率相等, 说明弹性模量与模型长度尺寸没有关系。
从对图5 (d) 的应力-应变分析可知, 钨合金塑性变形对温度的变化具有很高的敏感性, 但应力-应变变化的趋势基本相同。首先, 是线性上升, 为钨合金的弹性变形阶段。其次, 随着应变的增加, 体系进入塑性变形阶段, 且体系应力呈现一定的波动性。而随着温度的上升, 其屈服强度逐步下降, 屈服后应力波动幅度逐渐降低, 与宏观拉伸变形相似。
这是因为, 随着温度升高, 金属原子振动更加剧烈, 从而金属键合力减小, 金属发生软化现象, 当温度升高到一定程度时, 金属会发生一定的黏性流动, 且温度越高, 则黏性流动的抗力越小, 其发生蠕变的范围也越大, 蠕变很大程度上导致金属屈服应力等非弹性变形对应的流动应力下降, 但提高了钨合金的塑性。
图5 (e) 表明, 不同晶向的高应变率拉伸加载钨合金材料具有不同的等效弹性刚度和屈服强度, [111]晶向钨合金模型屈服强度最大, [110]晶向次之, 而[100]晶向钨合金模型屈服强度最低。[100]晶向的钨合金呈现最大屈服应变, [111]与[110]晶向的屈服应变近似, 可近似为[100]晶向的三分之一。晶向对力学性能的影响这一特性对于纳米材料在不同领域的应用具有极其重要的研究价值。
4.2 钨合金拉伸实验验证
4.2.1 准静态拉伸验证
实验材料选取91钨合金细化颗粒材料为对象, 进行了不同工况的实验研究, wW∶mNi∶mFe=91.0∶6.3∶2.7。它是将10.0~15.0μm的细化预合金粉末经固相烧结得到的钨含量相同晶粒度不同的细观粒状复合材料。图6为将烧结后钨合金材料 (棒材) 经线切割加工后成扁平哑铃型试件。在实验前需要把试件实验段的圆弧部分沿纵向研磨光滑, 以减小试件拉伸段圆弧部分的应力集中现象。
图7为采用S-570扫描电子显微镜试验系统对91钨合金细化颗粒材料试件扫描后, 给出的钨合金微观组织像。由图7可知, 钨合金钨颗粒为半径不同, 但近似球形, 相对均匀的被黏结相包围。颗粒与颗粒之间的界面接触面较小, 黏结相和晶粒之间的界面在所有界面中占多数, 还可以看到有孔洞等微缺陷。
准静态单向拉伸实验在MTS810液压伺服材料实验机上进行, 采用位移控制速率方式加载, 在试件实验段安装引伸计测量应变, 通过计算机输出实验结果数据文件。表1给出了不同应变率下的钨合金材料准静态屈服强度。
图8为91在不同应变率下, 钨合金材料准静态拉伸实验的应力-应变曲线。由图8分析可知, 随着拉伸应变率的增加, 91钨合金材料的强度有一定的提高, 具有典型的应变率强化效应;在不同的拉伸应变率加载条件下, 钨合金材料应力-应变曲线主要呈现了弹性变形阶段, 且材料伸长率达到0.153%时, 试件发生断裂, 几乎没有塑性变形, 据此可以判断此颗粒度的细化钨合金材料为典型脆性材料。
对准静态91钨合金细化颗粒材料在准静态拉伸加载破坏的试件材料断口利用电镜扫描图像 (见图9) 。对图像观测分析, 可以看到大部分钨颗粒发生解理断裂, 颗粒的解理断面上是扇形花样, 解理台阶为扇形的肋。扇形花样的肋线的汇合点, 为解理裂纹的起源点。河流花样以扇形的方式向外扩展形成所谓的扇形花样。从钨合金材料断口处还可以观察到相邻颗粒的接触面在受外载时分开后留下的痕迹。说明钨合金材料中同样存在少量颗粒与颗粒之间的W-W界面。
从图9中钨合金材料断口的显微图像上较为光亮的部分是材料的黏结相, 黏结相断口呈现纤维状, 有韧窝, 是韧性断裂。钨合金材料在外受载时, 黏结相协调了变形并重新分配应力, 使强度较高的钨颗粒承担了主要载荷;继续增大外载到一定值, 钨颗粒之中产生微裂纹, 裂纹扩展, 钨颗粒被解理。钨颗粒的解理断裂是材料失效的主要原因。
4.2.2 动态拉伸实验
图10为旋转盘式间接杆-杆型冲击拉伸实验加载装置示意图。采用该装置对钨合金试件进行加卸载实验, 其中加载装置为直径1.4m, 转动线速率可以达到100m/s的飞轮。
91细化钨合金材料动态拉伸实验是在室温 (14℃) 条件下进行的, 分别在1×10-3, 200s-1和500s-1三种应变率工况下进行动态拉伸实验研究, 在每一工况下进行了四组实验。图11为在不同应变率条件下实验获取的钨合金材料动态拉伸的应力-应变关系曲线。
对图11分析可知, 不同应变率动态拉伸导致了钨合金材料应力-应变曲线有一定的差异, 反映了91钨合金颗粒材料具有应变率效应, 且随着加载速率的提高, 屈服应力相应增加, 但是由于材料自身的缺陷, 导致了测量的数值具有一定波动, 也降低了钨合金材料强度的下降。
为了更好地了解钨合金材料的力学性能及在高应变率下钨合金强度变化规律, 采用分子动力学方法对该问题进行了模拟研究 (通过研究分析模型形状对强度变化影响可以忽略) 。
表2为在计算模型截面半径r=14a、长度l=40a、温度300K, 加载方向采用[100]晶向下, 不同应变率加载相对应的钨合金屈服强度数据。
通过对表2的数据进行平均差值及数据拟合得到了钨合金拉伸的应变率及屈服强度的幂律拟合公式。如下:
图12为通过拟合的公式与钨合金的实验数据对比。对图12分析可知, 低应变率拉伸加载区域拟合的较好, 且相同的应变率条件下, MD模拟对应的屈服强度高于宏观实验的屈服强度, 分析认为, 钨合金是完美晶体没有任何缺陷, 而宏观实验材料内部很可能存在微小的空洞、微裂纹以及杂质等外在的因素从而降低了钨的屈服强度。MD模拟钨合金得到的屈服强度很高, 这一结论与Михаловский等[15]实验测得拉伸钨晶须的屈服强度24.75GPa近似, 这说明采取的钨合金计算模型、弛豫时间、加载方式是合理的, 拟合得到的幂律公式是准确的, 尤其对于在实验中很难对高应变率拉伸加载测试手段获得相应的塑性力学参数, 分子动力学方法显得具有更加重要实际意义。
5 结论
(1) 钨合金材料在相同晶向、截面尺寸、温度的条件下钨合金随着拉伸应变率增加, 屈服强度升高;且高应变率下的钨合金在屈服时对应的屈服强度更大, 并且容易发生二次屈服。
(2) 相同晶向、截面尺寸、应变率的条件下拉伸钨合金随着温度升高, 屈服强度不断降低;而相同晶向、温度、应变率的条件下拉伸钨合金截面尺寸增加, 屈服强度有所增加。
(3) 在相同温度、应变率、截面尺寸的条件下, 晶向效应的改变对钨合金材料屈服强度和伸长率具有较大的影响。
金属Be物态方程的分子动力学模拟 篇9
1理论方法
分子动力学模拟方法是在蒙特卡罗罗模拟方法之后而产生的一种确定性的模模拟方法。它用来研究解析方法不能解决的的复合体系的平衡性质和力学性质,在材料料科学和计算机科学交叉学科中占有很重重要的地位。在分子动力学模拟之前,首先先需要建立一个模拟体系,其次列出体系内内粒子的运动方程,最后通过求解体系内粒粒子的运动方程来研究体系与微观量相关的过程。
2势函数及理论模型
在计算机模拟中势函数的选择直接关系到模拟体系是否能够较为客观准确地描述原子间的相互作用,因此势函数的选择很重要。本文采用的是修正型镶嵌原子法(即MEAM势),此势函数考虑了角度对电子密度分布的影响因子,其对电荷密度呈球对称的假设做了改进,因此MEAM势比EAM更准确。
由于金属Be常温下在400GPa内是六角密排结构 (hcp相),所以本文模拟的是金属Be的hcp相的MEAM势函数的参数。本文的模拟体系是一个周期性超晶胞,其由2048个原子组成,用速度verlet算法来完成体系的动力学方程的数值积分,计算选用的时间步长是1fs。模拟选用等温等压系综(NPT系综),温度和压强利用Nose-Hoover方法维持。
3结果与分析
3.1 T=300K 时的物态方程
利用分子 动力学的 方法,模拟了T=300K时hcp相的金属Be的物态方程。由于理论上计算的金属Be在T=300K时,400GPa左右时才 会发生相变。因此,我们每50 GPa模拟一次,模拟的压强变化范围是0到400GPa,每次统计的结果都是体系经过了充分的弛豫达到了平衡状态后的数据,模拟后的体积和压强的关系,用BirchMurnaghan物态方程(BMEOS方程 )即
进行了拟合 , 其中B0, B0'分别为体弹模量及其一阶导数。图1是拟合的金属Behcp相的等温(T=300K)物态方程。从图中我们可以看出在150GPa内 , 文献 [1] 中的实验数据与我们计算的结果吻合的很好,其实验数据是用冲击波约化的方法得出的结果;但文献 [2] 中的实验数据与本文的结果偏差稍大,其实验数据是采用静力学压缩的方法得出来的。研究人员通过对静力学压缩的实验数据进行了一些分析,发现静力学压缩的物态方程比其他实验的物态方程明显偏软 , 这可能是因为静力学压缩试验中测量的压力不是流体静水压造成的。而使用红宝石压力标度可对静力学压缩实验结果进行修正,这样可以减少实验数据之间的误差。由分析可知静力学压缩的实验结果有较大的误差,我们的计算结果与冲击波方法的实验值符合的很好。
3.2 体弹模量及其分析
由模拟计算T=300K时的体积与压强的关系,利用BMEOS方程拟合出了常温时六角秘排结构的金属Be的体弹模量B0 =113 GPa及其一阶导数B0' =3.97。本文计算值与其他值的比较如表1所示。由表可以看出,我们采用的MEAM势计算的金属Be hcp相的体弹模量与利用静态方法研究的实验值以及利用冲击波方法研究的实验值的差别都很小,与利用密度泛函理论计算结果误差分别为0.8%,体弹模量的一阶导数也界于实验值和理论值之间,差别很小。由数据我们可以看出室温时我们计算的结果无论是与实验值还是理论值符合的都很好,这也证明了采用MEAM势进行的分子动力学模拟的准确性。同时这为以后采用分子动力学方法对金属Be的其他性质的模拟计算奠定了基础。
4 结论
分子模拟 篇10
关键词：纳米薄膜,界面滑移现象,分子动力学
0 前言
目前, 微机电系统 (MEMS) 已经在产业化的道路上飞速发展, 而纳机电系统 (NEMS) 的研究仍处于起步阶段[1], 其中涉及到的润滑薄膜是近来研究的重点之一。在微米尺度下, 气体薄膜的流动已不遵守宏观流动规律。相对于气体分子, 液体分子间的距离更小, 微米尺度下液体流动的研究结果大多与宏观流动规律相吻合[2]。但是在纳米尺度下, 液体流动体系的尺度已经接近于分子的平均自由程, 所以常规的连续介质假定不再适用;另外由于尺度的减小, 影响流动的各种作用力 (如惯性力/粘性力、表面力/体积力等) 的相对重要性发生了变化[3,4]。
目前热门的生物芯片实验室 (Lab-on-a-Chip) 技术能够对单个生物分子进行检测和分析, 而液体常常是化学、生物和生理过程的实验环境, 以及能量、信息输运的载体, 其中涉及到极其重要的纳米尺度流体力学问题[5,6]。如固液界面的常规处理都是将界面附近流体的滑移速度认为是零, 在宏观尺度下, 这样的假设是合适的, 但是在流场本身尺度为纳米级别时忽略纳米级别的速度滑移是不可取的。另外由于流场表面积与体积比的增大, 与之有关的传热、传质过程及表面效应的作用大大增强, 其中的表面张力作用会影响到固液界面处的分子密度分布, 表面的吸附作用会使通道内流体流动的阻力损失较常规尺寸增大几十倍[7]。这些原本在宏观流场中忽略的因素都会凸现出来, 使得纳米尺度的流动规律发生变化, 从而影响芯片纳米通道内质量、动量和能量的输运性质。
对滑移现象的研究常以Couette流动为模型 (见图1) [8]。这种流动广泛存在于转子与马达底盘之间、梳状驱动器静止臂与活动壁之间以及轴承之中[9], 因此在纳米尺度下, 探索Couette流动是否遵守宏观流动规律非常具有研究价值, 对MEMS和NEMS中涉及的润滑薄膜设计起着指导作用。
由于纳米尺度的特殊性, 在实验中研究液体Couette流动具有较高难度, 目前该方面实验研究发表的文献较少, 如Haifeng Li等[10]采用多层纳米粒子图像测速技术测量出了Couette流动中距壁面80nm内颗粒的速度。在纳米尺度下, 液体Couette流动的研究方法主要集中在分子动力学模拟上。
本文在对国内外该方向研究成果进行总结的基础上, 介绍了纳米尺度下液体Couette流动的分子动力学方法, 着重介绍了分子动力学方法中涉及到的算法、系统模型、边界条件、势能函数以及目前研究者在这些细节上的不同选择, 并对他们所得出的结果做了总结。最后在目前研究现状的基础上, 对分子动力学在该方向未来的发展做了展望。
1 分子动力学方法
分子动力学 (MD) 诞生于1957年, Alder等[11]用MD研究了硬球系统的相变理论, 开创了利用MD方法来研究物质行为的先河。从此MD登上了科学发展史的舞台, 随着计算机技术的发展, MD被广泛用来研究材料科学、生物科学、纳米科学等不同领域。
分子动力学方法是建立在力场和牛顿力学基础上的计算方法, 与蒙特卡罗计算方法 (Monte Carlo method) 相比, 分子动力学方法中的粒子运动有正确的物理依据。此方法的优点为精准性高, 可同时获得系统的动态与热力学统计资料, 并可广泛地应用于各种系统及各类特性的探讨[12]。
Navier曾预言的液体与表面间“很小的滑移”是否存在及其发生机制的研究近几年受到了国际上微纳尺度流体力学领域的广泛关注[13]。学者们针对该问题开展了一些实验, 认为液体在固体表面的速度滑移是存在的, 并且与表面的浸润性质有关[14,15,16]。但是实验条件的不同导致实验结果差别较大, 这对机制的探讨造成了很大的困难。使用分子动力学方法探索纳米尺度的流体力学问题, 也证实了界面存在的滑移现象[17,18,19,20,21,22,23,24,25,26,27,28,29,30], 并且更好地揭示了现象背后的微观机制和规律, 所以在探索微观世界时, 分子动力学方法更具有优势。
由于计算机计算能力的限制, 分子动力学方法的研究尺度一般仅限于纳米水平。模拟开始于系统的N个分子, 按给定的密度值确定分子间距离, 并根据速度的概率分布函数和已知温度确定分子的初始速度, 从初始时刻出发, 沿时间轴以一定的步长对分子运动方程进行数值积分, 直到系统达到稳定[31]。
典型的模拟过程包括3个主要步骤[8]: (1) 确定初始条件和几何形状; (2) 确定分子间作用势; (3) 积分牛顿运动方程。
第一步中, 在一个二维或者三维的规则格子中放置一系列分子, 根据一定温度下的麦克斯韦尔分布, 给每一个分子赋一随机速度。第二步需要确定一个势能函数以描述分子或原子间的相互作用。势能函数可能是模拟一般流体的模型, 也可能是通过与实验比较, 迭代得到的用来模拟特殊情况的模型[9]。第三步是对运动方程进行积分, 是整个计算的核心。在分子动力学计算中需要解牛顿运动方程以得到速度和位置, 而粒子运动方程是一组常微分方程组, 求解时要采用有限差分法。由于作用力计算所占的工作量非常大, 故算法的选择非常重要。目前分子动力学模拟的算法主要有以下几种: (1) 预测-校正法[17,18,19,20,21]; (2) 跳蛙算法[13,22,23,24,25,26,27]和Verlet[28,29,30]。预测-校正法建立在对泰勒展开式进行微分的基础上, 而非直接由牛顿运动方程而来, 所以需要对预测的物理量进行校正。Verlet算法对速度的计算式为:
undefined
其中包括1/ (2δt) 项, 因为实际计算中通常选取很小的δt值, 所以计算速度时容易产生误差。跳蛙算法就是为矫正此缺点而提出来的。相对于预测-校正算法和Verlet算法, 跳蛙算法使用简单、计算存储量小、精度较高, 因此被广泛采用。
2 模型建立
2.1 系统模型与边界条件
关于Couette流动系统的模型设置如图2所示, 其中深色的为壁面原子, 浅色的为流体原子。两壁面平行于XY平面, 分别以相同的速率沿X或者Y方向作反向运动, 通道内便形成Couette剪切流动。在边界条件的设置上, Z方向为壁面边界条件, X和Y方向为周期性边界条件, 这种边界条件下, 当任一粒子移出计算域, 则立即有一粒子由相对的方向移入计算域[32,33], 如图3所示。这样的限制条件使得系统中的粒子数维持恒定, 密度不变, 符合实际的要求[12]。
2.2 固体壁面
为了真实反映流体的流动规律, 需要合理地构造固体壁面。分子动力学模拟中关于壁面原子的设置, 有以下几种常用的方法。 (1) Rigid-atom壁面[18], 以晶体结构来设置原子位置。若模拟的系统已有实验的X射线衍射晶体结构, 则建议以此结构为计算的起点。若无实验的结构可参照, 则通常以均匀分布的面心立方 (Face-centred cubic) 结构为起点[12], 将其固定在晶格点不动[13,18,20,23,24,25,26,27,28,29,30]。 (2) D.Einstein固体壁面, 基于Einstein固体理论, 固体原子在其晶格点附近以相同的频率作简谐振动[34,35], 简谐振动由连接原子与晶格点之间的虚拟弹簧力进行模拟[17,19,21,22], 相比第二种设置方式, 第一种设置方式简单, 计算量小, 并且在界面现象的模拟上能够取得与第二种方法相同的结果。
2.3 势能函数
除了壁面的构造方法之外, 流体分子之间以及流体与壁面之间的势能函数也是影响分子动力学模拟结果的重要因素。最常用的势能函数为Lennard-Jones作用势:
undefined
式中:r为分子之间的距离, ε为能量参数, σ为距离参数, r-12项为排斥项, r-6项为吸引项。具体流体分子之间的LJ势能作用参数可以从一些文献中查得, 流体与壁面之间的LJ势能作用参数可依据Lorentz-Berthelot原则[36]得到。
对于单原子液体分子, 大多数学者都采用常用的LJ势能函数[13,19,21,23,24,25,27,29,30];有些在吸引项和排斥项前乘系数, 以得到不同的排斥或者吸引作用[17], 有些对LJ势能函数添加修正项, 修正后的势能函数被称为移位形式的LJ势能函数:
undefined
添加的修正项可以保证在截断半径之外势能作用为零[18,20]。
对于长链分子, 有些采用FENE (Finitely extensible nonlinear elastic) 势能, 其势能函数为:
VFENE=-κ/2log[1- (r/r0) 2]
如Marek Cieplak等[21,28]对不同长链分子中原子之间的作用采用常用的LJ势能函数, 而对于同一长链分子中原子之间的作用则采用了2种不同的势能函数, 分别是Confining势能:
undefined
和FENE势能:
undefined,
以此来比较不同势能函数对模拟结果的影响。A. Jabbarzadeh等[22]对不同分子内的原子以及同一分子内被3个以上原子隔开的原子之间的势能作用采用移位形式的LJ势能函数, 同时对长链分子内存在的化学键伸缩、键角弯曲、二面角扭曲, 分别采用以下势能函数:
undefined,
undefined,
undefined。
对于液体分子与壁面原子之间, 有些研究者同样也采用常用的LJ势能函数[17,19,27,28,29,30];有些研究者为了研究不同壁面浸润性情况下的流场性质, 则采用了对LJ势能函数中的吸引项或者排斥项乘上系数[13,21,23,24,25];有些也采用移位形式的LJ势能函数[18,20,22]。
在势能函数的选择上, 对于长链分子, 由于化学键的存在, 需要考虑化学键伸缩势能、键角弯曲势能等的作用, 计算量较大。对于单原子分子, 基本上都选择常用的LJ势能函数或者其修正后的形式, 采用该函数在对粒子运动方程求解时非常方便, 同时又不人为引入假设, 是目前文献中应用最广泛的势能函数。
目前的模拟结果表明:
(1) 在Couette流动中, 液体分子在通道内的密度呈层状的涨落分布, 如图4所示, 这种分布与壁面性质有关, 与壁面的剪切速率无关[13,18,21,22,23,24,25]。
(2) 通过改变液体与壁面的势能作用可以产生滑移、无滑移和负滑移现象[13,19,22,23,24,25], 如图5所示。
(3) Couette流动中存在临界速度剪切率, 在该剪切率以下, 固液界面的边界条件基本上与Navier模型一致, 滑移长度仅与壁面性质有关, 减小εwf或者增大壁面原子密度都可以使滑移长度变大;超过临界剪切率, 滑移长度会迅速增大[20,29,30], 如图6所示。
(4) 在不同浸润性作用下, 滑移长度会随温度的升高而出现略微增大和减小2种情况, 如图7所示[25]。
(5) 滑移率随着膜厚的减小而迅速增大, 如图8所示。这主要是因为膜厚减小到分子尺度时, 流体发生类固化相变, 在剪切作用下固态薄膜就会被撕裂从而导致滑移。因此这种滑移率也可以作为衡量薄膜固化程度的定量判据[37]。罗彬宾等[26]在研究膜厚对固液比的影响时也得出类似的结论。目前的研究结果也有存在争议的地方, 如曹炳阳等[13,23,24]认为, 随着固液界面接触角的增大, 滑移长度会一直增大;而Roman S. Voronov等[29]则认为随着接触角的增大, 滑移长度不会一直增大, 其过程中会出现一段减小的区域。
3 结语