铝相变临界压力的分子动力学计算Word格式.docx
- 文档编号:21612842
- 上传时间:2023-01-31
- 格式:DOCX
- 页数:15
- 大小:264.03KB
铝相变临界压力的分子动力学计算Word格式.docx
《铝相变临界压力的分子动力学计算Word格式.docx》由会员分享,可在线阅读,更多相关《铝相变临界压力的分子动力学计算Word格式.docx(15页珍藏版)》请在冰豆网上搜索。
航空、建筑、汽车三大重要工业的发展,要求材料特性具有铝及其合金的独特性质,这就大大有利于这种新金属铝的生产和应用。
铝的应用极为广泛。
铝是一种非常年轻的金属,从1944年到1988年的44年间,世界铝年产量增加了19倍,铝发展如此之快,说明铝本身有很多优点,因此铝不断在各个领域中代钢代铜代木也代陶土,铝就发展的一个重要形式就是“取而代之”。
在铝短暂的发展过程中,需求起了巨大的推动作用,如二战时期由于制造飞机,推动铝急剧发展,1940年到1943年铝的产量由23万吨猛增到146万吨。
而1960~1970年间铝由446万吨增加至1020万吨,其原因是建筑业广泛以铝代钢引起的。
汽车工业在发达国家是重要的支柱工业,当前汽车工业激烈竞争中一个新的方向即汽车“轻量化”-“以铝代钢”是必然趋势。
随着铝在现代生产生活中的大量应用,对铝材料性能的研究也逐渐深入,铝在各种特殊环境下的性能特征成了研究的重点,尤其是铝在高压下的相变研究得到了广泛的关注。
相变是指当外界约束(温度或压力)作连续变化时,在特定条件(温度或压力达到某定值)下,物相突然发生改变。
在通常条件下物质有固、液、气三种状态,分别称为固相、液相和气相,在极高温高压下还存在等离子态。
从广义上讲,构成物质的原子(或分子)的聚合状态(相状态)发生变化的过程均称为相变,如从液相到固相的凝固过程、从液相到气相的汽化过程。
金属和陶瓷等固态材料在温度和压力改变时,其内部组织或结构会发生变化,即发生从一种相状态到另一种相状态的转变,这种转变称为固态相变[1]。
对于相变的研究,一般情况下可以从热力学、结构变化以及相变动力学几方面进行分类。
其中热力学方面主要从相变临界条件方面考虑;
结构变化方面主要研究相变的产物;
而相变动力学主要指的是相变的机制等。
本文研究的铝的相变属于固态相变的范围。
在许多情况下,材料的性能决定于材料中的某些特定结构和微观组织结构和分布。
这种组织结构包括晶体的形状和分布(晶粒和金相组织)、构成的晶体结构及其中的晶体缺陷、电子组态等,当然也包括它们当中的组织缺陷。
因此研究固态相变对控制金属、合金以及某些非金属材料性能有极为重要的理论和实践意义。
引起固态相变的原因非常多,其中很重要的一种就是高压引起的相变。
在理论研究方面,人们认为铝在高压下存在结构相变。
在20世纪80年代Mcmahan和Moriarity等[2]用GPT(TheGeneralizedPseudopotentialTheory)和LMTO(TheLinearMuffin-Tin-Orbitals)两种方法分别计算了铝的面心晶格结构(fcc)、六角密堆积结构(hcp)和体心晶格结构(bcc),得到fcc-hcp以及hcp-bcc高压结构相变时的临界压力值分别是360GPa和560GPa(GPT),120GPa和200GPa(LMTO);
Boettger和Trickey[3]用LCGTO(TheLinearCombinationofGaussian-TypeOrbitals)方法采用LDA近似后得到fcc-hcp以及hcp-bcc高压结构相变时的临界压力值分别是(205±
20)GPa和(565±
60)GPa,同时得到fcc-bcc高压结构相变时临界压力值是330GPa。
Lam和Cohen[4]利用从头计算赝势平面波(PW)的方法得到fcc-hcp及hcp-bcc高压结构相变时的临界压力值分别是200GPa和400GPa,同时得到fcc-bcc高压结构相变时的临界压力值为300GPa。
Mishrab和Chaturvedi[5]经过计算后认为fcc-bcc高压结构相变时的临界压力值为275Gpa。
可见,由不同的理论方法计算而来的相变点的差异还是比较大的。
实验上对于铝相变的研究也有很多。
1988年,Greene[6]等人用DAC(Diamond-Anvil—Cel1)方法研究了铝在高压下的相变,在当时的实验条件,压力加载最大只能达到220Gpa左右。
在这种实验条件下,Greene等人发现在小于219GPa(V/V0<
0.5)时fcc结构的铝并没有发生相变。
随着实验手段的不断进步,在2006年,Akahama[7]等人用DAC(Diamond-Anvil—Cel1)方法将压力值提高到300Gpa附近并考察了铝在这一加载过程的结构稳定性。
在这一工作中,Akahama等人发现铝在217±
10GPa附近发生了fcc-hcp的结构相变。
这是人们首次在实验中观察到了铝在高压下的结构相变现象,并且相变临界压力与Lam和Boettger的理论计算工作符合的非常好,这在一定程度上激发了人们对铝高压相变的研究热情。
事实上,要充分了解相变的微观机制和建立相应的理论,就必须从微观层面上研究结构相变,对结构相变做微观观测。
几十年来,虽然结构相变过程及其相变机制受到了广泛重视,但由于结构相变过程是一个突变过程,转变往往在ps量级完成,因而在实验上很难详细观察到相关的相变机制,尤其是微观(原子)尺度的变化机制。
相变机制往往是通过相变时一些宏观量的变化特征推测得到的。
这一情况随着近年来一些实验技术的飞速发展有所改观。
比如,利用超快激光技术来研究结构相变。
利用飞秒激光产生的纳秒以下脉宽的X光可以在材料发生结构相变时实时测量晶体结构的变化,这一技术在研究铁的冲击相变中得到了成功应用。
研究结构相变机制的另外一个非常重要的工具就是分子动力学计算机模拟。
分子动力学计算机模拟是研究分子、原子系统的有力工具。
它具有沟通宏观与微观的作用,其模拟结果对理论分析和实验都是有力的补充,因此它在物理学、化学、生物学和材料学等许多领域中得到广泛应用[8-11]。
当实验研究方法不能完全满足研究工作的需求时,用分子动力学计算机模拟可以提供实验上尚无法获得或很难获得的重要信息;
虽然分子动力学计算机模拟不能取代实验,但可以用来指导实验,并验证某些理论上的假设,以此促进理论和实验的发展。
特别是结构相变中有许多与原子有关的微观细节,在实验中不容易获得,而在计算机模拟中却可以方便地得到。
这种优点使分子动力学计算机模拟在研究结构相变时显得非常有吸引力。
随着计算机计算速度与计算容量的不断提高,现有的大规模并行分子动力学模拟可以处理多达百万的原子系统,在模拟结构相变的过程中,能够获取原子位置变化的详细信息,模拟结果可以和实验结果相校验,因此成为在微观尺度研究结构相变的有力工具。
本文即是用分子动力学计算机模拟方法研究了铝的相变临界压力等高压相变规律。
一.研究方法简介
1.分子动力学
(1)分子动力学简介
研究材料的性质需要更深入的认识材料的物理特性和物理规律,这就要求研究领域必须从传统的宏观层次逐步转入细观乃至微观层次,以便考察材料的微观作用机理。
受目前实验手段在时间和空间分辨率上的限制,现在还无法通过传统的实验方法直接获取原子尺度上的信息。
分子动力学较好地解决了这个问题。
分子动力学(MolecularDynamics,MD)是一门结合物理,数学和化学的综合技术。
分子动力学是一套分子模拟方法,该方法主要是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系综中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。
分子动力学(MolecularDynamics,MD)模拟是指对于原子核和电子所构成的多体系统,把其中每一个原子核视为在全部其它原子核和电子作用下运动,在此基础上求解运动方程。
通过分析各粒子的受力情况,用经典或量子的方法求解每个粒子在某一时刻的位置和速度,确定系统中各粒子的运动状态,然后使用一定的统计方法计算出系统的力学、热力学、动力学等性质。
分子动力学模拟具有沟通宏观特性与微观的作用,对于许多在理论分析和实验观察上都难以观察的现象做出一定的解释,因此它在物理学、化学、生物学和材料学等许多领域中得到广泛应用。
分子动力学是由Alder和Wainwrigh[12]在1957年首创的,他们首先在硬球模型下,采用经典分子动力学模拟方法研究了气体和液体的状态方程,取得了尝试性的成功。
但是由于计算机速度和内存的限制,早期模拟的空间尺度和时间尺度都受到了很大的限制。
计算机技术的飞速发展和多体势函数的提出与发展,为分子动力学模拟技术注入了新的活力。
1972年,Lees和Edwards[13]首先将分子动力学模拟应用到了非平衡的研究,进一步扩展了分子动力学模拟方法的应用范围。
它提供了一种研究高温高压下的液体和固体平衡特性的切实可行的方法,它不仅有助于理论研究,而且有助于解释实验中出现的现象,它是研究复杂的凝聚态系统的有力工具,是对物理系统的确定的微观描述,这个系统可以是少体系统,也可以是多体系统,并可通过哈密顿量、拉格朗日量或牛顿运动方程来描述。
由于分子动力学模拟是用运动方程来计算系统的性质,这一技术既能得到原子的运动轨迹,又能像做实验一样进行各种观察,其结果既有系统的静态特性,也有系统的动态特性。
对于平衡系统,可以在一个分子动力学观察时间内作时间平均来计算物理量的统计平均值;
对于非平衡系统,发生在一个分子动力学观察时间内(一般为几十个ps)的物理现象都可以用分子动力学进行仿真。
特别是许多与原子有关的微观细节,在实际实验中无法获得,而在分子动力学模拟中可以方便清晰地得到,这对理论和实验无疑是有力的补充。
(2)分子动力学原理
分子动力学(MolecularDynamics,MD)模拟是指对于原子核和电子所构成的多体系统,通过经典力学计算系统中各个原子的运动轨迹,然后使用统计方法计算出系统的力学、热力学、动力学性质的方法。
分子动力学方法中有两个基本假设,所有的分子动力学的模拟计算都是在这两个基本假设下进行的。
这两个基本假设是:
1、所有粒子的运动都遵循经典牛顿运动定律;
2、所有粒子间的相互作用都满足叠加原理。
这两个假设意味着分子动力学虽然是在原子层次研究问题,但它忽略了量子效应,是一种近似计算模型。
在分子动力学中,首先将由N个粒子构成的系统抽象为N个相互作用的质点,每个质点具有相应的坐标、质量、及成键方式,按体系温度根据Boltzmann分布随机给定各质点的初始速度,然后根据所选用的势函数来描述质点间的相互作用。
接着依据牛顿力学、根据所选用的势函数来计算各质点的加速度及速度,得到在系统演化时各时刻各质点的坐标和速度,设定时间间隔对这些计算结果进行保存。
最后可以对这些结果进行结构、能量、热力学、动力学等属性特征的分析,得到感兴趣的计算结果,从而得到系统在演化过程中的结构、能量、温度等属性的变化特征。
其优点在于系统中粒子的运动有正确的物理依据,准确性高,并且可同时获得系统的动态与热力学统计信息,可广泛地用于各种系统及各类特性的探讨。
缺点是模拟移动时间间隔不能过长,通常为1飞秒。
按时间步长为1飞秒计算,粒子移动103次模拟时间仅为1纳秒,尽管如此,对纳米尺度的系统进行模拟,分子动力学有着无可比拟的优势。
分子动力学模拟的实际步骤可以划分为四步:
1、设定模拟所采用的模型,给定初始条件;
使得体系粒子位置、速度要与模拟系统相容。
2、计算作用到每个原子上的力。
这一步是模拟的关键,在这其中要看所用到的势函数是否合理,能否精确地描述体系的性质。
3、趋于平衡的计算过程;
在这一步中我们要用一些系统控制方法使得体系符合模拟所要求的压力、温度等。
并且在程序的设计中算法要求精确、高效、简洁。
4、宏观物理量的得出及感兴趣量的提取;
分子动力学计算的结果是体系中粒子的位置、速度、加速度等量,我们需要从结果中分析全部或局部体系的状态,统计出相应的宏观量和体系特征。
(3)分子动力学模拟的基本方程
分子动力学假定原子的运动服从确定的描述,原子的运动和确定的轨迹联系在一起。
分子动力学方法中描述原子的运动方程主要有两种:
拉格朗日运动方程和哈密顿运动方程。
1.拉格朗日(Lagrangian)运动方程
对于由N个自由度的相互作用的质点构成的运动系统,拉格朗日方程为:
,i=1,2,…,n(1-1)
其中q为广义坐标,代表质点的空间位置,
为质点位置对时间的导数
由n个原子构成的模拟系统的拉格朗日方程在笛卡儿坐标系下可写为:
,i=1,2,…,n(1-2)
对于互不影响的粒子系统,可取拉格朗日量:
(1-3)
式中
为第i个粒子的质量,L为系统的动能,即系统中所有原子的动能之和。
如果考虑到原子间的相互作用,L可改记为:
(1-4)
式(2-5)中等式右端的两项分别代表系统的动能和势能。
代入拉格朗日方程中可得系统的牛顿运动方程:
,
(1-5)
即为原子i所受的内力,即由于系统中其它原子的作用而在原子i上体现出的合力。
在分子动力学模拟中,经典拉格朗日方程常被用来计算原子的运动,更适合于求解系统运动的过程。
在计算中能够施加外部载荷如外力、约束、边界条件等,在这些条件下通过求解原子的速度和位置,得到单个原子的运动轨迹,描述系统的运动过程,从而反映在系统整体特征上的变形、缺陷等。
2.哈密顿(Hamiltonian)运动方程
采用广义坐标和动量的形式来描述粒子的运动,则可以得到哈密顿形式的运动方程,求解多粒子系统的状态和演化过程。
含n个粒子的保守系统的系统哈密顿量在笛卡儿坐标系中可以写为:
(1-6)
式中粒子的动量
,哈密顿正则方程为:
(1-7)
(1-8)
给定系统的初始状态(初始位置和速度),则可以求解上述方程,得到系统中原子在任意时刻的位置和动量,得到运动轨迹,由统计平均得到系统的热力学表征。
保守系统的哈密顿正则方程适合于求解系统的热力学状态,如温度、热流动等。
另外,在实际模拟中,哈密顿运动方程可以在保守系统的基础上,通过改变系统状态变量,与外界进行能量交换,构成不同系综来进行模拟计算。
(4)分子动力学方法——程序流程
分子动力学模拟过程与真实情况下的实验过程非常相似。
在真实实验的进行过程中,首先要做的工作就是准备实验研究需要的材料,然后把实验仪器和测量仪器与被研究材料连接起来进行实验,在实验进行过程的某一时刻或时间段,测量我们感兴趣的物理量,进来得到被研究材料的某些特性。
与这一过程相似,在分子动力学模拟中也遵循同样的方法。
一个典型的分子动力学模拟的步骤如下:
1)读入指定运算条件的参数(如初始温度、粒子数、密度、时间等)。
2)体系初始化(即选定初始位置和速度)。
3)计算作用于所有粒子上的力。
4)解牛顿运动方程。
这一步和上一步构成了模拟的核心。
重复这两步直
至我们计算体系的演化到指定的时间长度。
5)中心循环完成之后,计算并打印测定量的平均值,模拟结束。
具体流程图如图1所示。
图1分子动力学模拟流程图
2.势函数
在分子动力学模拟计算中,原子间的相互作用由势函数来具体描述。
因此,势函数的选取是分子动力学模拟的核心问题,它的准确性对计算结果的精度影响极大,直接决定了模拟结果正确与否。
势函数一般表示为粒子坐标r的函数形式,它给出体系能量与粒子坐标的关系。
如果这个函数对体系的描述越准确,那么用它来计算所得到的体系的性质就与真实情况越接近。
(1)对势
在分子动力学模拟计算的早期,分子动力学模拟所采用的势函数一般是对势模型(Pairpotential)。
对势模型认为原子之间的相互作用是两两之间的作用,并认为系统能量为各粒子能量总和。
对势中一些参数主要通过宏观实验参数来确定,其中有平衡点阵常数、弹性常数以及空位形成能、内聚能和层错能等。
它们在特定的问题中均有各自的优越性,可以比较好地描述除金属和半导体以外的几乎所有无机化合物。
对势虽然形式简单,并且不能很好的描述金属及半导体的特性,但它的提出对势函数的发展及分子动力学的发展起到了不可替代的作用。
比较经典的对势模型有Lennard-Jones势、Buckingham势、Morse势及附加势等。
对势模型仅考虑原子对之间的作用力,无法考虑非金属材料化学键的方位变化,鉴于此,人们开始考虑粒子间的多体作用,构造出多体势结构,用它来精确描述更为复杂的系统。
(2)多体势
多体势于20世纪80年代初期开始出现,1984年,Daw和Baskes[14]首次对金属提出了唯像近似的嵌入原子(embedded-atommethod,EAM)。
该理论的基本思想是:
在一系列相互作用的原子中,一个原子对能量的贡献是该原子所在处电子密度的函数。
具体的操作方式是把晶体的总势能分为两部分:
一部分是位于晶格点阵上的原子核之间的相互作用对势,另一部分是原子核镶嵌在电子云背景中的嵌入能。
Baskes教授的领导下研究小组已经对元素周期表中常用的42种元素测出了EAM所必须的物理参数,该项研究为原子多体势理论的应用做出了极大的贡献。
在EAM模型中,系统总能量可表示为:
(1-9)
其中F是嵌入能,第二项是对势项。
为系统中所有其它原子的核外电子在i处产生的局域背景电子密度之和,表示为:
具体到用于模拟计算某确定晶格模型的时候,上式中的各参数必须根据晶格模型具体包含的原子及原子间的可能结合方式来确定,在确定这些函数的具体形式和相关参数时,通常不再过分关注他们所代表的实际物理意义。
实际应用中在这些函数形式中将包含一些可调参数,通过调整这些参数,使所得模拟计算结果和实验结果尽可能的一致。
在Daw和Baskes公布了原型EAM理论不久,1984年,Finnis和Sinclair[15]基于束缚能带理论发展了一种在数学上等同于EAM的势函数。
即将嵌入能函数设为平方根形式。
Acldand[16]等在此基础上通过拟合金属的点阵常数、弹性常数、压强体积关系、空位形成能、聚合能给出了适用于Cu,Al,Ni,Ag的多体势函数。
其中(4-7)式中的多体项及对势项分别为
(1-10)
在上式中,当
时,
;
当
。
为常数,且有
具体的取值因物质不同而不同。
基于EAM势的势函数还有很多。
由于这种势函数很好地描述了金属原子间的相互作用,因而成为描述金属体系最常用的一种势函数。
总体来说,分子动力学模拟中要求所采用的势函数能够尽可能详细地、真实地描述原子间的作用力,但是目前MD所采用的多体势是用验的或半经验的方法建立起来的,用其计算所得的理论数据与实验数据比较起来总会有一些差距,这在一定程度上影响和制约了分子动力学的发展。
因此,要使分子动力学技术在材料性质及规律的模拟方面有进一步的应用和发展,需要更新更优良的数值计算方法,同时,更有效的模拟理论也亟待出现。
本文所采用的势函数的具体形式及参数由Adams[17]等人给出,在Adams等人的工作中,对EAM势参数进行了拟合,拟合结果很好地满足了铝的晶格常数、结合能、体模量、弹性系数、声子谱、熔点、熔化潜热等,拟合结果表示在表1中,从表1中可以看出,该工作的拟合结果与已有的大量实验数据吻合的很好。
随后,在用分子动力学方法模拟铝的特性的工作中,该势函数被广泛地应用。
因此,在本文中采用这一势函数来描述铝在高压下的相变特征。
二.计算过程
1.搭建样品
分别利用分子动力学程序搭建样品:
(1)晶格常数为4.032Å
,每个坐标方向有20个原子的面心立方(fcc)结构的样品;
(2)晶格常数为3.147Å
,每个坐标方向有30个原子的体心立方(bcc)结构的样品;
(3)晶格常数a=2.854Å
,c=4.652Å
,每个坐标方向有20个原子的六角密排(hcp)结构的样品[18]。
表1Adams拟合的势函数计算的数据与已有实验数据的比较
2.用分子动力学程序计算冷能曲线
利用分子动力学程序,将三个样品沿三个方向等比例压缩,压缩至原来体积的0.42倍,即压缩至平衡晶格常数的0.75倍。
通过修改分子动力学程序,控制样品压缩过程,在压缩过程中记录样品的能量,压力等随体积的变化,绘制冷能曲线,冷能曲线即单原子能量和单原子体积的关系曲线。
图2fcc结构铝、bcc结构铝、hcp结构铝的冷能曲线
图3fcc结构铝、bcc结构铝、hcp结构铝的冷能曲线
图2、图3是本文计算的fcc结构铝、bcc结构铝、hcp结构铝在高压下的冷能曲线。
随着压强的增大,从图2所表示的冷能曲线可以看到,随压力的增大,铝在压缩过程中会发生有fcc-hcp的相变。
而从图3所表示的冷能曲线可以看到,随压力的进一步增大,在压缩过程中铝还会发生hcp-bcc的相变。
3.由冷能曲线判定铝的高压相变
从冷能曲线上可以判断铝在高压下有无相变发生及发生相变时相变点所对应的压力。
在相变点,以下两个条件同时满足:
1)两相的压力相等,
2)相变点的焓相等。
由条件1)得到
(2-1)
而通过条件2)得到
(2-2)
变
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 相变 临界压力 分子 动力学 计算