分子动力学模拟

分子动力学模拟（精选8篇）

分子动力学模拟 第1篇

1 分子动力学模拟方法

分子动力学是一套分子模拟方法, 该方法主要是基于牛顿运动方程计算模拟系统中所有粒子的运动。模拟过程中, 每个粒子在全部其它粒子所提供的经验势场作用下按牛顿定律运动, 其受到的合外力满足:

其中mi和ri为第i个粒子的质量和位置, Fij为第j个粒子作用在第i个粒子上的分子作用力, 其表达式为:

本文选用常用的Lennard-Jones (LJ) 双体势能函数, 具体表达式为:

其中:ε, σ分别为能量和长度常量;rij代表粒子i和粒子j之间的分子间距。

对Ar-Ar和Cu-Cu间的LJ势能参数可在文献中查到, 其中氩和铜交互部分的势函数参数由Berthlot混合法则求出:

其中:下标1代表流体Ar, 下标s代表固体颗粒Cu。

导热系数由Green-Kubo公式求出, 表达式为:

热流密度Jq的表达式为:

过剩能量Ei的表达式为:

其中, V是体积, KB是Boitzmann常数, T是温度。

2 模拟结果及分析

2.1对纯液氩的模拟

首先对液氩的导热系数进行模拟计算, 以验证模拟可靠性。模拟温度设为86K, 系综选正则系综 (NVT) , 粒子初始时按照面心立方结构FCC (Face Centered Cubic) 布置, 时间步长为2fs。模拟总步数60万步, 其中前10万步为弛豫过程, 以使系统达到平衡。截断半径取2.5, 对纯液氩包含的总氩分子数取4000个, 氩的LJ参数和分别取1.654E-21J和0.3405nm, 铜的和分 别取65.625E-21J和0.23377nm, Ar-Cu间的和分别为10.4184E-21和0.28714nm。模拟计算后得到液氩的导热系数为0.1264W/ (m K) 。

2.2对Cu/Ar纳米流体模拟及分析

对Cu/Ar纳米流体, 选取Cu纳米颗粒体积浓度为1%, 对应纳米颗粒直径为1.8nm, 用48个Cu原子代替54个氩原子, 模拟总原子数为3994个。通过计算得到纳米流体的的导热系数为0.143W/ (m K) , 相对于纯氩流体来说, 导热系数提升为原来的1.13倍。

通过跟踪模拟过程中氩原子的运动轨迹发现, 模拟过程中, Cu纳米颗粒对周围的液体分子存在吸附作用, 被吸附的液体分子在随后的模拟过程中不再离开纳米颗粒, 而是随着纳米颗粒一起运动, 且被吸附的液体分子排布不再是杂乱无章的, 而是倾向于固体的有序排布。

液氩的径向分布函数具有“短程有序, 长程无序”的特点, 而晶体Cu则体现出“长程有序”。当在纯氩基液中加入Cu纳米颗粒后, 纳米流体的径向分布函数在中和了晶体Cu的“长程有序”特点后, 表现出了既有液体的“短程有序”, 又有类似于固体的“长程有序”特点。

图1为在400ps, 700ps和1000ps时刻纳米颗粒表面吸附层的液体分子数密度分布。其中纵坐标p表示数密度 (N/m3) , 横坐标d (nm) 表示距纳米颗粒表面的距离。由图可知几个时刻的数密度分布大致类似, 都是在约0.23纳米处出现第一个波峰, 且峰值大小差异不大, 说明紧靠纳米颗粒表面的液体分子在被Cu纳米颗粒吸附以后就不会再离开;随后在约0.45处出现第二个峰值, 且峰值的大小随着模拟时间的推移逐渐减小, 说明在模拟过程中, 由于Cu-Ar的分子作用力远远大于Ar-Ar的分子作用力, 靠近Cu纳米颗粒的Ar分子逐渐被吸附在Cu粒子周围, 使得离纳米颗粒稍远的范围内液体分子的数密度值降低;当Cu-Ar间达到一定的距离后, 纳米颗粒对较远处的液体分子几乎没有影响。

3 结论

通过对Cu/Ar纳米流体导热系数的分子动力学模拟, 发现在纯氩基液中添加浓度为1%纳米颗粒后, 即可将导热系数提升为原来的1.13倍。在纳米颗粒表面会形成液体吸附层, 由于该吸附层的排布类似于固体, 使得纳米流体中的这部分流体的性质呈现出类似于固体粒子的现象, 吸附层的导热系数远远大于纯基液的导热系数, 从而导致纳米流体的有效导热系数提高。

参考文献

分子动力学模拟 第2篇

超临界NaCl水溶液的分子动力学模拟?

采用分子动力学模拟的方法对超临界NaCl水溶液的.微观结构进行了研究.模拟发现在所研究超临界条件下,密度的变化比温度的变化对超临界NaCl水溶液的微观结构影响更大.温度及密度对Cl--H2O径向分布函数的影响比对Na+-H2O径向分布函数的影响要大.超临界条件下,各gNa+-Cl-在0.261nm处出现峰值,表明Na+、Cl-之间发生了离子的缔合.超临界条件下,随温度增加,缔合作用增强;随密度增加,缔合作用减弱.本文工作为建立可适用于超临界条件下的电解质热力学模型提供了依据.

作 者：周健 朱宇 汪文川 陆小华 王延儒 王钧  作者单位：周健,朱宇,陆小华,王延儒,王钧(南京化工大学化工学院,江苏省化学工程与技术重点实验室,南京210009)

汪文川(北京化工大学化学工程学院,北京100029)

刊 名：物理化学学报  ISTIC SCI PKU英文刊名：ACTA PHYSICO-CHIMICA SINICA 年，卷(期)：2002 18(3) 分类号：O645.164 关键词：氯化钠   微观结构   超临界流体   分子动力学   分子模拟   电解质溶液   离子缔合  

分子动力学模拟 第3篇

相变储能技术是利用相变储能物质相变时吸热和放热来储存和释放能量的技术。相变储能技术在太阳能、建筑节能和工业余废热等方面都有广泛应用[1]。石蜡是常用的有机相变储能物质,主要由直链烷烃混合而成,通式为CnH2n+2。在实际应用中,相变储能物质应具有与工况匹配的相变温度,研究表明有机相变物质的混合比例不同,相变温度也不同。因此,研究有机相变物质混合比例,有助于相变储能技术的进一步应用。目前确定有机混合物质相变温度的方法主要有差示扫描量热法 (DSC) 和动态热机械分析法(DMA)[2],但这些方法所用设备昂贵,不利于大范围推广。分子动力学模拟(MD)方法是一种广泛使用的计算机模拟方法。近年来,不少学者采用MD方法,利用双体分布函数(PCF)和自扩散系数研究了金属相变的微观结构和热物性[3,4,5]。也有文献报道了对有机小分子的相变性质进行模拟的成果[6,7,8]。本研究采用MD方法,根据体系自扩散系数和比体积随温度变化的关系来确定各体系的相变温度,并根据十七烷二面角的变化研究了在熔化过程中体系微观结构的变化。

1 判断依据及势能模型

自扩散系数D被广泛用于描述体系是否发生相变[9,10],其数值大小反映了体系内粒子扩散行为的强弱[11]。自扩散系数D随温度变化出现转折时,即可判断体系发生了相变,从而确定体系的相变温度。

根据爱因斯坦扩散定律,自扩散系数D可表述为均方位移MSD随时间的变化,其关系式如下:

undefined

undefined

式中:MSD是粒子位移平方的平均值;undefinedi(t)与undefinedi(0)分别表示粒子i在t时刻和零时刻的位置矢量。

比体积是说明物质在某一状态下分子疏密程度的物理量,也常用于判断体系是否发生相变[12,13,14]。当比体积随温度明显变化时,体系疏密程度剧烈变化,即体系发生了状态的变化。

本研究采用ORGANIC势函数,该势函数包含了分子间作用力L-J力场和分子内作用力OPLS力场。

体系势函数可以表示为:

Etot=Einter+Eintra (3)

分子间作用势:

undefined

分子内作用势:

Eintra=Eb+Eθ+Eφ (5)

式中:Eb表示键伸缩能,Eθ表示键角弯曲能, Eφ表示二面角扭转能。具体势函数及参数见表1。

2 模拟方法

采用分子动力学方法进行模拟研究。模拟压力为常压,采用周期性边界条件,截断半径为0.77nm。各模拟体系初始密度均为1g/cm3,各体系简称及组成见表2。

先将初始密度均为1g/cm3的体系置于NTP系综下弛豫10ps,温度为100K,然后将弛豫后的体系置于NTV系综下弛豫10ps,温度仍为100K。重复上述操作,弛豫时间总计为40ps,使体系充分弛豫达到平衡。图1为十七烷弛豫前后体系分子的分布结构示意图。

最后将各体系弛豫结构作为起始结构,在NTP系综下对各温度进行模拟,温度范围为310～270K,温度间隔为5K。模拟总时间为200ps,时间步长0.2fs。为了确保数据的可靠性,从100ps开始统计所需参数,共获得5001个数据点。

3 结果与讨论

3.1 十七烷、十五烷及十七烷/十五烷混合体系相变点的判断

模拟得到Hep体系和Pen体系在298K时的密度分别为0.7878g/cm3和0.7814g/cm3,与文献[15]的实验值吻合较好(如表3所示),证实了模拟方法的可靠性。

图2和图3分别为Hep和Pen体系自扩散系数D和比体积随温度变化的关系。从图2可知,Hep体系自扩散系数在温度小于294K时变化不大,在温度大于294K后出现了明显的增大,说明Hep体系中粒子的扩散性随温度的升高而加强,特别是在温度大于294K后,扩散能力急剧变化,即体系发生了相变;从图3可知,Hep体系的比体积在296K处出现了明显的转折,表明在该温度点体系疏密程度发生了变化,即出现由固态向液态的转变。在Hep体系中,比体积和自扩散系数随温度变化有几乎相同的规律,二者的转折点均在295K附近,这与十七烷实验相变温度295K[15]一致。同理可以确定Pen体系的相变温度为284～283K,这也与十五烷的实验相变温度283K[15]一致。

图4和图5分别为8Hep+2Pen、6Hep+4Pen、4Hep+6Pen、2Hep+8Pen体系自扩散系数及比体积随温度变化的曲线。根据图4和图5可以确定8Hep+2Pen、6Hep+4Pen、4Hep+6Pen 及2Hep+8Pen二元混合体系的相变温度分别为285～283K、281～280K、280K、282～281K。

结合上述结果绘制了十七烷/十五烷二元混合体系的相变温度随十五烷分子个数变化的关系曲线,如图6所示。

在十七烷中加入一定量的十五烷后,所形成的二元混合体系的相变温度远低于十七烷的相变温度,甚至比十五烷的相变温度还低,最低可至280K。随十五烷分子个数的增加,混合体系的相变温度逐渐接近十五烷的相变温度。这说明十五烷的加入对十七烷的相变温度有明显的影响。闫全英等[16]通过DSC方法研究了低温相变石蜡,发现当组成二元混合体系的纯物质的相变温度相差大于10K时,混合体系的相变温度将发生较大的改变,与本实验的研究结果类似。随着十五烷分子个数的不断增加,所形成的二元混合体系的相变温度并无太大变化,波动范围为3K。由图6可知,对于十五烷/十七烷二元混合体系,理论上可以获得280～296K温度区间的不同比例的混合物。

3.2 十七烷熔化过程微观结构的变化

以十七烷为例分析了纯物质在熔化过程中微观结构的变化。由热力学理论可知,物质的熔化过程是熵增大、有序度降低的过程。分子有序度与其构象有着对应关系。直链烷烃由于碳碳单键的自由转动可以形成许多不同的分子构象,表征这些分子构象最直观的参数就是由相邻4个碳原子所形成的二面角(Dihedral angel)。当直链烷烃形成180°的二面角时,整个分子表现为全反式构象,此时分子体系有序度最高。

根据MD模拟得到的分子内坐标,统计相邻4个碳原子形成的二面角,结果见图7。从图7可知,温度较低时,十七烷的二面角与全反式时的180°很接近,此时的结构可看作是十七烷熔化前的稳定结构,其有序度最高;随着温度的升高,二面角逐渐偏离180°,在310K时已经接近160°。因此,十七烷分子在升温过程中其微观结构不断发生变化,分子有序度不断降低。

4 结论

(1)采用分子动力学模拟方法和ORGANIC势函数研究烷烃大分子体系的相变行为是可行的。

(2)十五烷/十七烷二元混合物的相变温度为280～296K。可以通过调节组分比例,得到不同的相变温度,实现有机相变储能的要求。

(3)十七烷在熔化过程中其微观结构不断变化,有序度降低。

摘要：采用分子动力学模拟(MD)和NTP系综方法对十七烷、十五烷及十七烷/十五烷混合体系的相变行为进行了模拟研究。根据体系自扩散系数和比体积随温度的突变关系,获得的十七烷、十五烷相变温度值与文献中的实验结果吻合较好,并确定了不同混合体系的相变温度。二元混合体系的相变温度可根据组分比例进行调节,可满足一定温度下相变储能的应用要求。通过分析熔化过程中十七烷分子二面角的变化获得了十七烷分子微观结构变化的特征。在熔化过程中十七烷分子有序度降低。

分子动力学模拟 第4篇

1 模型和方法

1.1 模型描述

无定形二氧化硅的初始构型包含1200个氧原子和600个Si原子,任意排列在一个边长为30.1的立方盒子里,其密度为2.2g/cm3,该密度与无定形Si O2实验测得的密度相当。原子间的相互作用采用Born-Mayer-Huggins[1]势能来描述。

其中,Aij和ρ分别为短程相互作用参数;Zi,Zi是指点位偏移电荷;βij代表长程相互作用参数。这些参数列于表1中。模拟采用周期性边界条件(PBC),5.5的球形截断半径。

1.2 模拟方法

采用正则系综(NVE)和微正则系综(NVT)的分子动力学模拟方法,模拟熔融-淬火[2]方法降温。模拟体系首先在8000K下,进行(20ps NVT+8ps NVE)分子动力学模拟,使固体完全熔化;然后逐渐降温经过一系列退火过程(6000K-4000K-2000K-1000K-300K)得到室温下(300K)低能量的无定形二氧化硅骨架结构。在300K下再进行20ps的NVT,并统计得到无定形二氧化硅的径向分布函数和角度分布。

在模拟池z-方向的中间切割无定形二氧化硅材料得到破裂的无定形二氧化硅表面。在切割的二氧化硅中插入一个50的真空层,以消除上下表面的相互作用,然后进行升温柔性化操作:首先采用MD模拟技术(20ps NVT+8ps NVE)在2000 K进行柔性化模拟,然后逐步(1000K-800K-500K-300K)冷却至室温300K得到柔性化的无定形二氧化硅表面。

2 结果和讨论

2.1 无定形二氧化硅骨架结构图

采用分子动力学模拟,得到一个边长为30.1的无定形二氧化硅立方形模拟晶胞。经过熔融-淬火方法得到的无定形二氧化硅骨架结构如图1所示。红色代表氧原子,黄色代表四面体。

2.2 径向分布函数

对得到的无定形二氧化硅在室温300K下再进行20ps的NVT分子动力学模拟,并统计得到它的径向分布函数,如图2所示。我们的统计结果与实验[3]测得的结果一致,说明我们的模拟结果是可靠的。图中g Si O(r)在1.6处出现一个突出的峰,这个距离对应于Si-O键的键长。

2.3 角度分布

图3给出了O-Si-O的键角分布情况。在图3中宏观的无定形二氧化硅的O-Si-O键角分布在109°附近出现最高峰,并且在峰的两侧呈现对称分布。

2.4 无定形二氧化硅表面

直接破裂得到的无定形二氧化硅表面如图4(a)和4(b)。其表面非桥氧(NBO)浓度为3.31个/nm2。对破裂表面进行升温柔性化得到表面如图4(c)和4(d),其表面NBO浓度为0.99个/nm2。另外表面还存在一些不饱和硅、二元环、三元环等缺陷,其表面真实的羟基化Si-OH)浓度比NBO浓度要大一些。

无定形二氧化硅表面的性质与热处理的方法和过程有关,实验研究给出了在不同条件下得到的不同性质的表面。对于水热合成方法得到的无定形二氧化硅表面的羟基(Si-OH)浓度一般在4-7个/nm2之间。这与直接破裂得到的亲水性表面可以比较。实验干氧氧化方法得到的无定形二氧化硅表面羟基化浓度(2.6个/nm2)。在本文的模拟中因为表面柔性化温度高于硅玻璃的转换温度得到的非桥氧浓度与实验中干氧氧化得到憎水的表面相当。

3 结论

采用微正则系综(NVE)和正则系综(NVT)的分子动力学模拟,熔融-淬火的方法得到了无定形二氧化硅结构。计算得到的不同原子的径向分布函数和取向分布和文献一致,说明得到的无定形二氧化硅结构合理。对该结构直接切割得到无定形二氧化硅表面,接近实验中水热合成法得到的亲水的无定形二氧化硅表面;采用分子动力学模拟对表面进行升温柔性化得到的无定形二氧化硅表面,接近于真实实验干氧氧化得到憎水的表面无定形二氧化硅表面。

参考文献

[1]Rosenthal,A.B.;Garofalini,S.H.J.Am.Cream.Soc.1987,70,821.

[2]Webb,E.B.;Garofalini,S.H.J.Non-Cryst.Solids.1998,226,47.

分子动力学模拟 第5篇

1 材料与方法

1.1 AgSnO2触点样品的反应合成与加工制备

以Ag粉、AgSn合金粉、氧化剂为原料, 在冷等静压机中压成素坯。素坯在真空反应合成烧结炉中, 烧结成SnO2含量为10% (质量分数, 下同) 的AgSnO2复合材料锭坯 , 其反应合成工艺为:室温450℃ (保温lh) 850℃ (保温3h) 随炉降至室温;真空度<210-2Pa。将锭坯在温度为650～850℃、压力为1000～1200MPa的条件下挤压成线坯, 经拉拔成ϕ=1.35mm, 真应变ε=11.7的丝材。用YFC-16冷墩复合触点机设备制备得到AgSnO2触点样品。

1.2 电接触实验

电接触实验条件为:直流、阻性负载, 电压18V, 对应电流为20A, 闭合力0.8N, 频率为1.2Hz, 触点间距1.30mm, 实验次数:10000次。比较阳极、阴极转移后的触点表面形貌及成分变化, 分析触点内部发生的微观变化机理。

1.3 分析测试

采用SEM, TEM和EDAX对经过电接触实验的AgSnO2触点样品进行显微组织和元素的观测。

1.4 分子动力学模拟建模

模型建立:本工作采用改进的Lennard-Jones (L-J) 双体势能模型[12,13,14], 系统模型为长方形控制盒, 如图1, 2所示, 其中x=50nm, y=200nm, z=50nm。建立Ag, SnO2的晶胞结构, 将Ag原子、SnO2分子置于计算对象中, 详细参数设置如表1所示。对电接触过程中熔池的各组分形态进行分析, 选定对Ag (液) +SnO2 (固) 模型和Ag (气) +SnO2 (发生分解) 模型进行计算;考查Ag (液) +SnO2 (固) +Ag (气) +SnO2 (发生分解) 等各物相对熄弧的过程过程的影响, 同时, 在每次计算之前均给予充分的弛豫步数20000步, 使模型状态更接近实际, 该过程中对有关物理量进行实时统计和记录, 当系统内产生明显的气液界面, 气固界面时, 即可判断灭弧过程正在发生。

以两相共存AgSnO2电接触复合材料反应合成烧结态显微组织为基础, 氧化锡均匀弥散分布在银基体中;反应合成方法制备的材料界面清洁, 同时由于电接触过程中, 熔池区域小, 因此电弧作用的界面区域内可忽略由于浓度梯度引起的界面的运动, 设置为Ag, SnO2两相的直接接触界面;考虑熔池的大小, 对含4000～6000原子的模型进行计算。

2 结果与讨论

2.1 电弧侵蚀表面形貌特征

反应合成AgSnO2材料阳极触点电弧侵蚀后的触点表面可观察到烧蚀痕迹和熔融金属流动现象, 如图3所示, 触点的表面形成了熔池, 熔珠, 气孔和坑洞, 可见阳极触点表面电弧侵蚀严重, 发生过剧烈气化、沸腾。这主要是由于在触点反复闭合、分断的过程中电弧在高温作用下, 触点表面形成的熔池温度分布极不均匀, 在熔池中心发生金属的熔化、气体的吸收;而在熔池的边缘却发生金属的凝固、气体的逸出, 形成了凝固后留下的气孔和孔洞。材料在发生阴极转移的情况下, 阳极表面的孔洞更多, 还可以观察到阳极表面有明显裂纹, 在AgSnO2内部组织中, Ag与SnO2的相界面两物相的膨胀系数相差较大和气孔是产生表面裂纹的根源。图4为电弧侵蚀前后形貌图, 发生侵蚀前氧化锡颗粒是弥散分布的, 说明模型1建立的基础是正确的, 发生电弧侵蚀后, 弥散分布的氧化锡在在晶界处聚集并明显长大。

对触电表面各点 (图3中对应各点) 进行了能谱分析, 如表2所示。从元素成分分析得知, A点为银流动冲刷的SnO2颗粒。B点为阴极触点发生损耗后在触点表面残留的小坑。C点为发生阴极转移后, 阳极触点表面气孔处的熔体, 成分为纯银, D点为银流动产生的富集区。从图3和表2得知, 反应合成AgSnO2材料在电接触实验条件下, 发生阴极转移后AgSnO2材料侵蚀程度大于阳极转移。

2.2 Ag (液) +SnO2 (固) 模型的熔池行为分析

图5所示为Ag (液) +SnO2 (固) 模型的直观3D原子运动示意图, 由图可见, 随着计算步数的增加, 银液对基体中团聚的氧化锡冲刷程度也加剧, 体系在运行500000步时, 银液在界面处把团聚的氧化锡颗粒分成了两部分, 这表明在银熔池中, 界面处银的扩散对于SnO2颗粒的冲刷作用明显, 可以将大团的团聚颗粒解离成较小的团聚颗粒。

团聚颗粒的流动、细化可以更好地增大熔池的黏性, 降低了材料的喷溅损耗, 最终提高反应合成AgSnO2电接触材料的抗电弧侵蚀能力。如果氧化物与基体的浸润程度不高, 这种冲刷作用带来的负面影响容易使氧化物在晶界和表面聚集 (见图4) , 从而增大接触电阻, 使与其相连的富Ag区熔焊趋势增大, 产生熔池、孔洞等表面特征。

在电弧作用下, 当AgSnO2触点发生电接触的表面累积的电弧能量达到氧化锡的分解所需的能量后, 此时溶池中银组分已经汽化。

2.3 Ag (气) +SnO2 (发生分解) 模型与实验分析

Ag (气) +SnO2 (发生分解) 模型的计算结果 (见图6) 直观示意出了反应合成AgSnO2触点材料中SnO2组分发生分解过程的3D原子运动过程示意图。

氧化锡分解成锡和氧之后, 氧会散失到空气中, 部分锡也将发生挥发。如下所示:

undefined

升华和蒸发都会导致气孔和孔洞产生 (图3) , 表2中D处与A, B两处相比, Sn, O量明显减少, 表明在电弧侵蚀作用下SnO2发生了分解, O随熔池的冷却而逸出, 而Sn则在电弧作用下挥发。模拟与实验结果一致。

通过对两种模型计算, 得到各原子的均方位移曲线, 表征粒子在各个时间步中的位移平方。如图7和图8所示, 该曲线的斜率表征了各元素的扩散系数。Ag (液) +SnO2 (固) 模型中 (见图7) , 体系中银的运动最剧烈, 氧和锡的运动非常缓慢, 并且几乎没有相对运动, 这表明第二相氧化锡未发生分解 ;Ag (气) +SnO2 (发生分解) 模型中 (见图8) , 氧化锡中锡和氧的均方位移很大, 且曲线基本呈直线上升的趋势, 这表明氧化锡已经发生了分解, 且形态已经到达气态。这些都会消耗电弧能量, 是发生灭弧的主要原因。

3 结论

(1) 反应合成的AgSnO2电接触材料在电弧侵蚀后, 触点表面呈现气孔、熔珠、熔池和裂纹形貌特征。在18V, 20A的工作条件下, 发生阴极转移后AgSnO2材料侵蚀程度大于阳极转移。

(2) 出现的电弧侵蚀形貌形成机理为:当起弧温度界于银基体的熔点和氧化锡的熔点之间时, 团聚的氧化锡颗粒一方面在银熔池中被不断冲刷、分散成较小的团聚颗粒, 增大了熔体的粘度, 另一方面使得颗粒在晶界或表面聚集;随着温度升高, 当熔池温度高于银的升华温度后, 氧化锡将发生分解, 分解的氧散失在环境中, 分解的部分锡也将发生挥发。这些物质状态的变化产生了气孔、熔池、裂纹等电弧侵蚀形貌特征。

(3) 通过分子动力学计算结果可知:可增多反应合成AgSnO2电接触材料的氧化物组分, 通过氧化物发生物态变化、化学分解过程来达到熄弧的目的, 改善氧化物与基体的浸润程度, 避免电弧侵蚀过程中Ag冲刷带来氧化物的大量聚集, 提高反应合成AgSnO2电接触材料的电接触性能。

摘要：研究电接触材料电弧侵蚀作用过程中灭弧过程的微观机理, 采用分子动力学和实验的方法对反应合成的AgSnO2电接触材料的熔池行为进行模拟。结果表明:根据起弧时的物相所构建的模型, 能够准确地模拟熔池特征, 熔池内部物质运动和存在状态是决定耐电弧侵蚀的关键因素。可通过增加反应合成AgSnO2电接触材料与基体浸润的氧化物组分, 以及氧化物发生物态变化、化学分解过程来达到熄弧的目的。

分子动力学模拟 第6篇

超强酸作为替代浓硫酸的催化剂的研究已引起人们的关注。超强酸按负载物和载体可分多种, 其活性和选择性同制备的原料 (前驱体) 、制备方法和助剂的添加有较大的关系[5,6,7]。研究较多的超强酸是SO4-2/ MxOy型金属氧化物超强酸[8]。

SO4-2/ MxOy体系的酸性与表面S-M配合物种类有关, 它是以金属氧化物为底物, SO4-2为助剂, 生成硫助金属氧化物, 而不是金属硫酸盐[9]。

超强酸是指酸强度比100%的硫酸还要强的酸, 即哈梅特酸度函数HO< -11.94的酸。由于其能在较温和的条件下活化酸催化反应, 还具有选择性高, 副反应少, 易和反应物分离, 对设备不腐蚀, 可重复使用, 废催化剂处理简单, 对环境污染少等优点, 而迅速成为催化剂领域中的研究热点, 并代替传统的酸催化剂广泛应用于有机反应中。超强酸作为一种新型催化剂材料, 是一种对环境友好的绿色化学催化剂[10]。

选用超强酸类催化剂有:SO4-2/ Fe3O4, SO4-2/ TiO2, SO4-2/ ZrO2, 利用hyperchem 7.5软件, 在opls力场下, 选择分子动力学 (MD) 方法将以上超强酸类催化剂用于模拟合成巯基乙酸甲酯, 通过对势能数据的分析, 得出巯基乙酸, 甲醇与超强酸类催化剂之间相互作用的适宜温度以及最佳催化剂。

1 实验材料和方法

利用hyperchem 7.5软件构建出SO4-2/ Fe3O4, SO4-2/ TiO2, SO4-2/ ZrO2, 甲醇和巯基乙酸结构, 然后分别保存为单体结构数据文件。

利用hyperchem 7.5软件, 在OPLS力场对上述一种单体结构、合并的两种单体结构, 以及合并的3种单体结构等3种单体各种可能的情况进行构象优化, 收敛标准为0.1kcal/mol (0.418 kJ/mo1) 。然后以优化后的结构为基础, 分别对各种可能的情况在OPLS力场进行分子动力学模拟, 模拟条件为:温度从270K到370K, 每隔10K为一个模拟温度;时间步长为0.5 fs (1 fs=10-15s) 。模拟时间共20.2 ps (1ps=10-12s) , 用0.1ps从0 K升到模拟温度, 在模拟温度进行20 ps模拟, 再用0.1 ps从模拟温度降到0K, 记录模拟过程势能数据[11]。

巯基乙酸和甲醇分子在发生反应或其他过程中, 其状态应具有所有可能的构型。因此为了得到能够确实反映模拟温度时的能量, 舍弃模拟温度进行前一半 (即0.1～10.1ps, 2万种状态) 的势能数据, 对模拟温度进行后一半 (即10.1～20.1ps, 2万种状态) 的势能数据进行平均, 作为该温度下的势能。这样既保持反应在模拟温度下有足够多的状态数, 同时又有效避免刚开始到达模拟温度时而反应可能还未达到平衡状态的影响。其相互作用能可由下式得到:

undefined

其中:E为体系的总能量;E1, E2为单体能量。相互作用能ΔE为负值时, 表示两者相互吸引作用大, 体系能量降低多。负值的绝对值越大, 说明相互作用越强。

2 结果与讨论

根据上述方法, 计算出不同温度下的巯基乙酸, 甲醇与超强酸类催化剂的势能以及它们相互作用能ΔE, 作出结合能温度关系图1-4并在下面分别进行讨论。

2.1 无催化剂时巯基乙酸与甲醇的相互作用分析

在不考虑催化剂作用的情况下巯基乙酸与甲醇相互作用时温度对相互作用能的曲线如图1所示。

从巯基乙酸与甲醇作用能图 (图1) 中, 可以看见270, 290, 300, 310, 360, 370K巯基乙酸与甲醇作用能数值均为负值, 而且在330K以后随着温度的升高, 巯基乙酸与甲醇作用能数值逐渐减小, 在370K这个温度点达到最低点ΔE=-34.91253086, 说明巯基乙酸与甲醇反应随着温度升高相互作用就越强。

2.2 考虑催化剂SO4-2/ Fe3O4, 巯基乙酸与甲醇的相互作用分析

SO4-2/ Fe3O4, 巯基乙酸与甲醇相互作用时温度对相互作用能的曲线如图2所示。

从SO4-2/ Fe3O4, 巯基乙酸与甲醇作用能图 (图3) 中可以看到:270, 290, 300, 340和370K的作用能数值为负值, 说明在这些温度下SO4-2/ Fe3O4, 巯基乙酸与甲醇的作用比其他温度下作用强烈, 且340K时SO4-2/ Fe3O4, 巯基乙酸与甲醇作用能数值最小ΔE=-16.00555627, 说明SO4-2/ Fe3O4, 巯基乙酸与甲醇作用最强烈, 即340K为SO4-2/ Fe3O4, 巯基乙酸与甲醇反应的最适温度。

2.3 考虑催化剂SO4-2/ TiO2, 巯基乙酸与甲醇的相互作用分析

SO4-2/ TiO2, 巯基乙酸与甲醇相互作用时温度对相互作用能的曲线如图3所示。

从SO4-2/ TiO2, 巯基乙酸与甲醇作用能图 (图4) 中可以看出:所有温度的作用能数值均为负值, 说明SO4-2/ TiO2, 巯基乙酸与甲醇作用强, 而且在310K时出现SO4-2/ TiO2, 巯基乙酸与甲醇作用能数值最小值ΔE=-35.59376049, 说明在此温度下SO4-2/ TiO2, 巯基乙酸与甲醇反应作用最强烈, 故而310K为SO4-2/ TiO2, 巯基乙酸与甲醇反应的最适温度。

2.4 考虑催化剂SO4-2/ ZrO2, 巯基乙酸与甲醇的相互作用分析

SO4-2/ ZrO2, 巯基乙酸与甲醇相互作用时温度对相互作用能的曲线如图4所示。

从SO4-2/ ZrO2, 巯基乙酸与甲醇作用能图 (图4) 中可以看到:除310, 330K和350K作用能数值为正值外, 其余温度作用能数值均为负值, 说明在这三点温度下SO4-2/ ZrO2, 巯基乙酸与甲醇作用不如其他温度下强烈, 同时在320K时SO4-2/ ZrO2, 巯基乙酸与甲醇作用能数值最小ΔE=-23.14282134, 在此温度下SO4-2/ ZrO2, 巯基乙酸与甲醇反应是最强烈的, 即320K是SO4-2/ ZrO2, 巯基乙酸与甲醇反应的最适温度。超强酸类催化剂特征性能见表1。

3 结论

(1) 在没有催化剂的情况下, 巯基乙酸与甲醇的反应作用缓慢, 加入了超强酸类催化剂后反应强烈。

(2) 通过巯基乙酸、甲醇与超强酸类催化剂作用能数据显示:超强酸类催化剂加入到巯基乙酸与甲醇的酯化反应中, 其反应的适宜温度范围为310～340K。

(3) SO4-2/ TiO、巯基乙酸与甲醇作用能数值ΔE的负值绝对值最大, 故而理论上相对于其他的超强酸类催化剂而言, 选用SO4-2/ TiO2较理想。

参考文献

[1]MATSUBARA Y, NAKAMURA T, YOSHIHARA M, et al.Anovel method for the preparation of thioformates and thiols byusing thioformimidates[J].Chemical&Pharmaceutical Bulletin.1982, 30 (9) , 3389-91.

[2]JIANG Hua, DIAO, Kai-sheng, PAN Ping-lai, et al.Study onthe stability of rhodium complex with uncoordinating donor andthe property of catalysis[OL].Chemical Journal on Internet[on-line computer file].2000.

[3]俞善信, 文瑞明, 龙立平.固体酸催化合成乳酸正丁酯[J].应用化工, 2000, 29 (4) :2-5.

[4]林昌志.巯基乙酸甲酯的合成研究[J].安徽化工, 2004, (3) :26-27.

[5]HINO M, ARATA K.Solid superacid form sulfuric acid and Zir-miun hydroxide for acetic with ales[J].J C S Chem Comm, 1980, (2) :851-853.

[6]ARATA K.Preparation of superacids by metal oxides for reac-tions of butanes and pentanes[J].Applied Catalysis A:General, 1996, 146 (1) :3-6.

[7]MORTERRA C, CERRATO G, SIGNORETTO M.On the roleof the calcination step in the preparation of active (superacid) sulfated zirconia catalysts[J].Catal Lett, 1996, 41 (1-2) :101-104.

[8]姚胜.SO4-2/MxOy型固体超强酸催化剂[J].化学通报, 1990, (2) :23-28.

[9]曾健青, 罗庆云, 王琴, 等.SO4-2/MxOy型固体超强酸催化剂的研究进展[J].石油化工, 1994, 23 (3) :191-197.

[10]汪显阳.固体超强酸SO4-2/TiO2-ZrO2催化合成对羟基苯甲酸正丁酯[J].精细石油化工进展, 2004, 5 (5) :42-45.

分子动力学模拟 第7篇

关键词：中子辐照,分子动力学,空位浓度,位移级联,间隙原子团

反应堆压力容器(RPV)是在高温、高压等强烈的辐照等恶劣条件下运行的,A508-Ⅲ钢是国产反应堆压力容器用钢,具有较高的抗中子辐照脆化性,设计寿命预计在40年以上[1]。辐照脆化是制约RPV寿命的重要的性能劣化机理。中子辐照使RPV钢的力学性能劣化,其劣化的程度由多种因素决定,如中子剂量、温度、中子的通量和钢材料中有害杂质(Cu,P,S等)的含量。为了确保核反应堆的安全运行,需要监测反应堆压力容器(RPV)材料中子辐照引起的脆化程度,并对整个寿期内RPV的辐照稳定性进行评估[2,3,4]。 辐照损伤会在反应堆压力容器钢中产生大量空位、间隙原子、Frank缺陷对、层错四面体等缺陷,影响压力容器钢的使用寿命。近年来,随着国家核电事业的发展,核电关键材料的损伤机理也得到了广泛的研究。如三维原子探针分析结果表明,富铜原子团簇的数密度为1023m-3数量级,富铜原子团簇的直径在1～3nm[5],初步建立了适用于国产RPV辐照脆化预测模型[6];Calder等用分子动力学模拟了bcc-Fe辐照级联的过程[7];Marian等用分子动力学模拟了 Cu在α-Fe的扩散和析出过程[8]。

本工作利用分子动力学(MD)模拟了0.0%～0.9%空位浓度下中子辐照级联过程,深入讨论分析了不同空位浓度对级联过程的影响。

1 模拟方法

采用改进的F-S势(Finis-Sinclair Functions)[7],对bcc-Fe晶体分别在X(100),Y(010), Z(001)方向上进行堆垛。图1为初始模型堆垛图,模型体积为(50a050a050a0)nm3,约250000个原子。晶格常数a0=0.287nm。PKA原子入射角度为[135]方向,体系初始温度为300K,时间步长选取1fs(110-16s),PKA能量为5keV。对模型分别随机去掉0%,0.1%,0.2%,0.9%的原子形成空位。形成空位浓度分别为0%,0.1%,0.2%,0.9%的初始模型, 分别对四个模型进行模拟。

图2为识别空位和间隙原子的示意图。在辐照级联的模拟过程中,产生大量的空位和间隙原子等点缺陷。为了观察方便,需要将晶格上的原子略去,只保留间隙原子和空位。实现方法是,以晶格位置为圆心,以0.3 a0为半径,做一个球体,如果在这个球内没有原子,这个晶格点即为一个空位(见图2(b));同样,以每个晶格点为圆心做一个球,如果某个原子位置均不在所有的球内,则这个原子是一个间隙原子(见图2(c))。将计算得到的数据进行以上处理,即可得到任意时刻的空位和间隙原子的图像。

2 模拟结果与讨论

当空位浓度为0%时,即纯铁的模拟中,辐照级联过程首先由初始碰撞原子(Primary Knocked-on Atoms,PKA)碰撞基体内Fe原子,

被碰撞的原子获得能量后继续碰撞下一个原子,从而产生位移级联过程。图3为纯Fe的辐照级联过程,PKA能量为5keV,入射方向为[135], 温度为300K。在位移级联过程中随着模拟时间的增加,点缺陷数量NF不断增加,当经过约0.5ps时,达到最大值NF=1632,然后部分点缺陷被空位湮灭,数量随之减少。经过3～4ps后逐渐稳定下来,10ps以后基本不再有太大变化,此时NF=60。从图3 (e),(f)观察到,在级联过程中产生的离位峰效应[9]使间隙原子向区域的外部扩展,空位向级联区域的中心不断聚集。公式(1)是诺盖特罗宾逊和特伦茨(Norgett Robinson and Torrents,NRT)用来估计辐照金属中辐照剂量(Displacements Per Atom, DPA)的标准。

undefined

式中:NNRT表示NF的值,k是近似等于0.8的常数。Ed是所有可能的晶向上的平均离位阈能(Ed值为35eV,是在所有可能的晶向上的平均值[7,9,10,11]),E为弹性碰撞产生位移级联的PKA能量,取E=5keV,由公式(1)算出NNRT的值为57。NF为最终由辐照产生的Frank缺陷数量即NF最终稳定下来的值,模拟中10ps后的NF值为60,可见与理论估计值相当接近。

在辐照级联的MD模拟中发现,中子辐照主要产生大量的空位和间隙原子等点缺陷,点缺陷的构成主要结构有哑铃型结构,哑铃型结构能量最低,因此在辐照过程中产生大量哑铃型结构,并且在模拟结束后,点缺陷基本以哑铃型结构存在。图4为间隙原子团和空位团的形成状态,某些空位会形成空位团,间隙原子会形成间隙原子团,如图4(a),(b)所示。在模拟过程中偶尔会产生间隙原子形成的八面体间隙原子团结构,见图4(c)。对于bcc晶体,四面体间隙大于八面体间隙,因而,自间隙原子只可能存在于四面体间隙位置[12],八面体间隙原子团由四个正六边形面和四个正三角形面组成,每个六边形面和三角形面的方向均为〈111〉方向,八面体间隙原子团共由12个原子堆垛产生。在对纯铁及空位浓度为0.1%的Fe的模拟中均出现了间隙八面体原子团结构。而在对空位浓度为0.2%,0.9%的Fe的模拟中没有发现间隙八面体原子团结构。可见,在低浓度空位下容易产生八面体间隙原子团,高空位浓度下,间隙原子团较难形成。这是因为空位浓度的增加使点缺陷中间隙原子与空位复合的几率增加,使间隙原子更容易因复合而湮灭,很难形成团簇结构。当八面体间隙原子团产生后一直到模拟过程结束,其结构均没有太大变化,可见八面体间隙原子团是以一个较稳定的结构存在。

图5为在不同空位浓度下辐照级联过程中NF随时间变化的关系。随着模拟时间的增加,NF逐渐升高,到0.5ps左右时达到峰值NF(max),随后逐渐减少,经过约4ps后逐渐稳定下来,10ps以后基本不再有太大变化。可见在相同能量的PKA级联下,空位浓度对级联过程表现出积极的影响。

图6为不同空位浓度下经过10ps后的状态。在相同能量的PKA级联下,空位浓度对级联过程表现出积极的影响。当空位浓度为0.1%时,经过10ps,Frank缺陷对为NF=38;当空位浓度为0.2%时,经过10ps,Frank缺陷对为NF=13;当空位浓度为0.9%时,经过10ps,Frank缺陷对为NF=6。由于增

加点缺陷的密度可加速金属原子和置换原子的扩散[13],所以空位浓度增加后,使最终剩余因辐照产生的点缺陷数量NF明显降低,从模拟过程中我们也观察到预置空位加快了级联过程中点缺陷的湮灭。随着空位浓度的增加,级联过程中产生的点缺陷湮灭的速度也加快并导致最终产生点缺陷数量NF值降低。

3 结论

(1)纯Fe在中子辐照级联过程中产生大量空位、间隙原子及空位团、间隙原子团。四面体间隙原子团是一种比较稳定的结构状态。辐照级联模拟过程中产生的离位峰作用使间隙原子向外部扩展而使空位向中间聚集。辐照级联中产生的点缺陷浓度随PKA原子能量升高而升高。

分子动力学模拟 第8篇

Pb-Cu偏晶合金具有优异的性能,在结构材料、电触头和超导材料方面具有广阔的应用前景[1]。关于Pb-Cu合金的制备、微观组织结构及性能已进行了大量的研究工作[2,3,4],研究表明,金属材料的宏观性能主要是由其微观结构决定的,而微观结构与合金凝固前的熔体结构及熔体组元的扩散性质有很大关系,因此研究Pb-Cu合金的液态结构和合金的扩散性质可以使科研工作者更好地理解以及充分利用Pb-Cu合金。

与固态相比,金属熔体中原子没有恒定的格点位置,其结构存在不稳定性和不确定性,难以用一个很好的图景来描述。近几十年来,经过众多学者的努力,已积累了金属熔体及其过冷态的大量实验数据,从而为在原子和电子层次上探讨液态金属和合金的微观信息及其动力学性质奠定了基础。目前的研究工作表明,分子动力学方法模拟在揭示金属熔体的结构演变、非晶倾向以及动力学性质计算等方面具有巨大的发展和广阔的应用前景[5,6,7]。

液态金属扩散过程的研究在材料科学、流体物理学和电化学等众多领域中起着极其重要的作用[8,9]。几十年来,液态扩散系数的测量方法一直是人们的研究热点,但是由于合金系众多和实验的复杂与困难,只测定了少数合金系的扩散系数[10,11]。自20世纪80年代以来出现了采用分子动力学方法进行金属单质自扩散系数的研究,并取得一定的进展。Belashchenko[12]模拟了1423～7400K时Cu单质的自扩散系数,1423K时模拟得到的自扩散系数为0.358Å2/ps,实验值为0.471Å2/ps[13]。Dalgic[14] 模拟了1273K时Ag单质的自扩散系数,得到的模拟值为0.3628Å2/ps,实验值为0.281Å2/ps[15]。Akhter[16]的研究表明,1500K以上Ag自扩散系数的模拟结果较好,但模拟Ag的熔点温度为1415K,与实验值1234K相差较大,在熔点处的模拟结果不符合实际情况。目前使用分子动力学方法仅限于对金属单质的扩散性质进行模拟,尚无二元合金扩散性质的研究报道。

由Darken方程可知,合金的互扩散系数取决于组元的自扩散系数和热力学因子,合金的热力学因子与合金组元的活度系数紧密相联,由此本研究提出了一种计算二元合金互扩散系数的新方法。首先采用分子动力学方法计算二元合金中各组元的自扩散系数,通过NRTL方程得到二元合金的热力学因子,进而结合Darken方程计算合金的互扩散系数。

本研究通过上述方法计算了Pb-Cu合金的互扩散系数,并模拟液态Pb-Cu合金的对相关函数、配位数及相关半径,用于分析合金的熔体结构及熔体结构对液态合金扩散性质的影响。

1 分子动力学模拟

分子动力学模拟的关键在于选择准确的势能模型来描述分子间相互作用。由X.W.Zhou提出的普适嵌入原子模型(GEAM势)[17,18]可以用于若干金属和合金的分子动力学模拟,是近年来常用于描述过渡族金属原子间相互作用的势能模型,其基本方程为:

undefined

式中:Φij(rij)为相距为rij的2个原子i、j之间的对势作用;Fi(ρi)为原子i的嵌入原子能;ρi可以表示为undefined为与原子i相距rij的原子j在i原子所在位置产生的电子密度;嵌入原子能Fi(ρi)可以表示为:

undefined

undefined

undefined

式中:Fni、Fi、Fo、η 为模型参数,可由结合能和体模量计算得到;ρe 为平衡时的电子密度。两体作用势为:

undefined

式中:re 为最近邻原子间的平衡距离; A、B、α、β为4个模型参数;κ、λ、m、n是有关截断半径的附加参数,一般情况下m=n=20,λ=2κ,α/β=1.875。

本研究的分子动力学模拟在容纳2048个原子的立方晶胞中进行,使用周期性边界条件。由于涉及的Pb、Cu两种元素均为面心立方结构,因此模拟初期阶段立方晶胞中Pb、Cu原子以面心立方结构随机排列。在目标温度(1323K)下Pb-Cu合金体系首先采用NPT系综运行3104个时间步,以确保合金体系冲破笼子束缚效应,充分驰豫达到平衡,而后转入NVT系综运行3104个时间步,用于收集所需数据。时间步长选择2fs。在NVT系综下进行数据收集的原因是:当采用NPT系综进行模拟时,由于原胞体系不恒定,体积变化较大,因此对金属原子的扩散系数产生影响,表现为均方位移(MSD)曲线不稳定,会出现在一定的时间段内MSD曲线突然减小的情况;当采取NVT系综进行模拟时,由于立方盒子体积恒定,MSD曲线线性增大。故采用NVT系综来模拟液态金属的MSD曲线。

2 液态合金结构的模拟

液态合金的微观结构对合金的性质有着重要的影响,本研究将使用对相关函数、合金配位数和相关半径来分析液态合金的熔体结构。

图1为模拟得到1323K时Pb20Cu80、Pb35Cu65、Pb65Cu35、Pb90Cu10 4种合金的对相关函数曲线。根据合金的对相关函数曲线计算得到合金的配位数列于表1。

表1中Ntotal为合金的平均配位数,Ni-j为原子i周围存在的j原子的偏配位数。从表1中可以看出,合金中同种原子间的偏配位数NPb-Pb+NCu-Cu远大于异种原子间的偏配位数NPb-Cu+NCu-Pb,表明熔体内部异种原子的成键数目远小于同种原子的成键数目,在熔体内部的局部范围内分别形成Pb团簇和Cu团簇,与文献[19]中Pb-Cu合金形成液相互不溶区的结论一致。同时从图1中可以看出,随着Pb含量的增加,合金gCu-Cu的第一峰高度显著增加,表明熔体中Cu团簇的有序度不断增大。为了进一步考察合金熔体内部的团簇现象,本研究计算了合金的相关半径及团簇的原子个数。相关半径rc在液态结构中代表原子团簇的尺寸,是原子团簇尺寸的下限。理论上,在r较大的地方双体分布函数g(r)=1,原子全密度函数ρ(r)=ρ0,液态中的原子呈统计分布,原子之间的相关性消失。相关半径rc的定义为:

g(rc)=1±0.02 (4)

表2列出了计算得到的Pb-Cu合金的相关半径和团簇内的原子个数。其中rc,alloy为合金的相关半径,Nc,alloy为合金团簇内部的原子个数。rc,Pb和rc,Cu分别为Pb组分和Cu组分的偏相关半径,Nc,Pb和Nc,Cu为Pb和Cu团簇内部的原子个数。

从表2中可以看出,在Pb35Cu65、Pb65Cu35合金中合金团簇内部原子个数远小于富Pb端Pb90Cu10合金和富Cu端Pb20Cu80合金,表明在中间浓度范围内Pb、Cu原子间易形成较大的无序区。当Pb含量为90%时,Cu团簇的相关半径和团簇内部的原子个数均急剧降低,表明在此成分下合金熔体中Cu团簇突然变小,个数增多,Cu团簇在熔体内部呈微观均匀分散分布。

3 Pb-Cu合金扩散性质的模拟

3.1 合金自扩散系数的计算方法

合金中各组分的自扩散系数可通过Einstein方程计算合金各组分的均方位移再对均方位移求斜率得到。

均方位移(MSD)为:

undefined

自扩散系数为:

undefined

式中:ri(t)是原子i在t时刻的位置,N是原子总数,D为自扩散系数,〈〉代表所有原子在整个时间段的系综平均。

3.2 合金互扩散系数的计算方法

根据Darken公式,二元合金体系的互扩散系数可表示为:

D=[Dundefined(1-xi)+Dijjxi]ft (7)

undefined

式中:fi为热力学因子,ai和γi分别为组元i的活度和活度系数,xi为组元i的摩尔分数, Dundefined、Dundefined分别为组元i、j在熔体中的自扩散系数。 可以看出只要得到合金的热力学因子fi,就可以计算二元合金的互扩散系数。本研究采用NRTL方程计算Pb-Cu合金的活度系数,进而计算合金的热力学因子fi。

NRTL方程是在Wilson的局部组成概念和Scott的双流体模型基础上导出的,该模型对于模拟二元、三元甚至多元合金的活度系数具有较高的精确度, 满足合金热力学性质计算的要求。对于i-j二元合金体系,其组元i活度系数表达式为:

undefined

式中:undefined为原子的相互作用能;αij=αji,使用Pb-Cu合金的无限稀活度系数通过牛顿迭代法拟合得到τij、τji、αij三参数, 采用式(8)、式(10)计算合金的热力学因子fi。

undefined

undefined

表3列出了计算得到的NRTL方程的三参数和热力学因子fi,用于计算互扩散系数。图2和表4列出了计算得到的Pb-Cu合金的自扩散系数和互扩散系数。

文献[20]中通过对Pb-Cu合金偏结构因子的计算,得到了合金中两组分的自扩散系数和互扩散系数(见表4)。通过对本研究计算值与文献值的比较可以看出,对Pb组分自扩散系数的计算两者较近似,但对Cu组分自扩散系数的计算相差较大。同时,文献[20]中没有考虑热力学因子对合金互扩散系数的影响。

从图2和表4中可以看出,Pb含量对Pb-Cu合金的互扩散系数有较大的影响,Pb90Cu10合金中Cu的自扩散系数和互扩散系数远大于其它成分的相应值,而Pb65Cu35合金中Cu的自扩散系数小于其它成分的相应值。

本研究从合金熔体结构的角度分析产生上述现象的原因。一般而言,扩散系数与原子间的异种键成键数有关。这是由于异种键易引起原子周围附近应力场发生畸变,易导致形成空位,并降低扩散激活能。同时熔体内部局部范围内组分的浓度差越大,合金的热力学因子越大,从而导致合金的互扩散系数变大。

从表1中还可以看出,Pb90Cu10合金异种原子间的偏配位数(NPb-Cu+NCu-Pb)在1323K时为7.6,远大于其它成分下异种原子间的配位数,说明在Pb90Cu10合金中Cu与Pb原子的成键数目较多,Cu被Pb原子包围,从而导致Cu的自扩散激活能降低,使Pb90Cu10合金中Cu的自扩散系数大于其它合金。此外,在1323K时Pb90Cu10合金中Cu的偏相关半径突然降低至8.05,团簇的原子数骤减到13,说明此时合金熔体中Cu原子形成较小的团簇,团簇的个数较多,各个Cu团簇较分散地分布在合金熔体中,易在熔体内部多个局部范围内形成较大的浓度差,导致热力学因子增大,从而使合金的互扩散系数增大。以上因素均导致Pb90Cu10合金的扩散系数增大。

从表4中还可以看出,在Pb65Cu35合金中Cu自扩散系数突然变得极小。由表2可知1323K时Pb65Cu35合金中异种原子间配位数(NPb-Cu+NCu-Pb)为4.72,小于其它成分的NPb-Cu+NCu-Pb值,表明此成分下Cu-Pb异种键成键数最小,导致合金Cu自扩散激活能增大,合金的Cu自扩散系数降低。通过以上分析可以看出,本研究对液态Pb-Cu合金互扩散系数的计算具有一定的合理性。

4 结论