中国有色金属学报 2004,(02),255-261 DOI:10.19476/j.ysxb.1004.0609.2004.02.019
δ′相在外场下早期沉淀机制的计算机模拟
王永欣 陈铮 刘兵 马良 唐丽英 赵宇宏
西北工业大学材料科学与工程学院,西北工业大学材料科学与工程学院,西北工业大学材料科学与工程学院,西北工业大学材料科学与工程学院,西北工业大学材料科学与工程学院,西北工业大学材料科学与工程学院 西安710072 ,西安710072 ,西安710072 ,西安710072 ,西安710072 ,西安710072
摘 要：
基于离散格点形式的微扩散方程(Langevin方程)和非平衡自由能函数,编制了引入原子间相互作用能变化的Al3Li(δ′)相沉淀原子层面计算机模拟程序。该程序包容亚稳区到失稳区的全部温度、成分范围和孕育期至粗化的全过程,可以处理与时间相关的过程问题。开展了不同原子间相互作用势下原子图像、序参数的计算机模拟,进而探讨了最近邻原子相互作用能(W1)对有序相沉淀的影响机制。发现过渡区合金形核前出现短暂的等成分有序化阶段。探明随W1的增大,有序相沉淀的孕育期缩短,形核率增加,合金有序化速度和原子簇聚速度加快,在所研究的时间范围内达到的长程序参数和成分偏离序参数最大值增大。随W1的增大,有序相的析出呈现出失稳分解的特征。
关键词：
原子间相互作用势 ;沉淀 ;原子图像 ;序参数 ;计算机模拟 ;
中图分类号： TP391.9
作者简介： 王永欣(1971),男,讲师,硕士.;
收稿日期： 2003-04-30
基金： 国家自然科学基金资助项目(50071046);
Computer simulation on precipitation mechanism of δ′ phase in early stage in external energy field
Abstract：
An atomic-scale computer simulation program of δ′(Al3 Li) based on the microscopic diffusion equation and non-equilibrium free energy was firstly worked out, which can be applied to process the problem relating with time. The precipitation mechanism of δ′ was investigated by simulating the atomic pictures and calculating the order parameters, etc, and the variation in precipitation mechanism with the nearest interchange interaction energy(W1 ), particularly the precipitation behavior of 12%Li in Al-Li alloys was interpreted. It is found that the precipitation incubation period of δ′ is shortened, the nucleation rate is increased, and the rate of ordering process and atom clustering are increased while the nearest interchange interaction energy increases. Within the studied time period, the maximum values of long range ordering parameter and composition deviation parameter increase with the nearest interchange interaction energy increasing. The precipitation of δ′ occurs either by a congruent ordering process followed by, or by a non-classical nucleation mechanism. With increasing W1 , the precipitation of δ′ leans to decomposition.
Keyword：
interchange interaction energy; precipitation; computed microstructure; order parameters; computer simulation;
Received： 2003-04-30
计算机模拟作为科学研究
[1 ]
的重要手段已经应用于很多方面, 分子动力学方法和蒙特卡洛方法是其中应用最为广泛的两种方法。 分子动力学方法(MD方法)又可以分为2类, 其一为经典分子动力学方法, 其二是第一性原理分子动力学方法(CP方法)。 蒙特卡洛方法采用针对某种事件出现的概率建立一个随机模型, 然后制造一系列随机数以模拟这个过程, 最后进行统计性处理
[2 ]
。 在MD方法的处理过程中, 方程组的建立是通过对物理体系的微观数学描述给出的。 该方法可以处理与时间有关的过程问题, 因而可以处理非平衡态问题。
但是该方法的计算程序复杂, 计算量很大, 内存占用大。 蒙特卡洛方法是实现Gibbs统计力学的途径, 在处理基态、 平衡态或近平衡态、 均匀体问题时取得了很大成功。 随机模拟的方法的计算程序简单, 计算量小, 内存占用小, 但是该方法难以处理非平衡态问题。 相析出的早期阶段, 体系可能出现很大的浓度梯度, 可能从母相中产生随机分布的多个新相胚芽, 体系处于高度非线性、 高度非平衡态。 并且新相胚芽的形成过程恰恰是被关注的焦点问题。 此时MD方法(包括CP方法)以及蒙特卡洛方法的应用受到局限。
采用基于离散格点形式的微扩散方程(Langevin方程)和非平衡自由能函数的计算可以比较好地解决上述问题。 该计算模型在计算时进行了傅里叶变换和投影, 使得原本直接在实空间求解一个关于三维、 非线性问题的偏微分方程问题, 转变为在倒易空间求解一个关于二维、 线性问题的常微分方程问题, 使模拟计算量大为减小, 时空间尺度分别达到10-6 ～10-3 s, 10-7 ～10-6 m。 本文以二元铝锂合金为研究对象, 编制了引入原子间相互作用能变化的Al3 Li(δ ′)相沉淀原子层面计算机模拟程序, 开展了原子间相互作用势对原子图像、 序参数影响的计算机模拟。 在合金处理过程中, 施加外势场属于先进的处理手段, 从根本上讲, 外势场影响的是原子间相互作用势。 因此本文进而探讨了最近邻原子相互作用能(W 1 )对Al3 Li相沉淀的影响机制。
1理论模型
经典动力学理论将形核、 长大和粗化过程作为独立的阶段处理, 并在每个阶段中作了某些假设。 如假设沉淀相的原子结构已知, 新相与母相为均匀连续介质, 且被明锐相界面分开。 而这些假设是值得商榷的, 形核、 长大和粗化过程并不能截然分开, 应当作为互相重叠的伴随过程来处理
[3 ]
。 这使得经典形核理论的适用范围受到限制。
真实相界面在结构及成分上是具有一定厚度的过渡层。 Cahn和Hilliard称之为界面的弥散性质, 提出了非经典形核理论
[4 ]
, 并发展了匀相分解的连续模型, 即Cahn-Hilliard方程
[5 ]
。 但是该模型不能描述原子簇聚与有序化共存的过程, 如Ni-Al和Al-Li 两相合金中L12结构沉淀相从无序基体中析出的过程。 因此, 欲描述原子簇聚与有序化共存的过程, 必须将方程置于离散格点形式的基础之上。
离散格点形式的微观随机扩散方程是Cahn-Hilliard方程的微观形式, 又称微观Langevin方程, 它是在离散格点形式的微观扩散方程基础上, 添加一随机力项ξ (r t )而得到的。 它用原子占据晶格位置的几率描述原子组态和相形貌:
d Ρ ( r , t ) d t = c 0 ( 1 - c 0 ) Κ B Τ ? ∑ r ′ L ( r - r ′ ) ? F ? Ρ ( r ′ , t ) + ξ ( r ′ , t ) ? ? ? ( 1 )
上式右侧对晶体中所有晶格位置求和, L (r r r r T 为温度, K B 为玻尔兹曼常数, c 0 为基体平均浓度, F 为系统的总自由能, P (r t )表示原子在t 时刻、 占据晶格位置r
离散格点形式的微观Langevin方程解决了描述原子簇聚与有序化共存过程的问题
[6 ,7 ]
, 而且无须预先设定沉淀相的结构, 唯一需要输入的变量为原子间相互作用势, 计算机自动模拟可获得包括瞬时相在内的有序结构
[8 ]
。 同时该模型考虑了沉淀相的粗化过程
[9 ,10 ]
, 克服了传统粗化模型将第二相假设为由纯溶质原子构成, 且两相界面明锐, 以及第二相体积分数无限小等缺点, 充分体现了形核、 长大和粗化过程在时间上相互重叠、 伴随的特点。 对于Al-Li二元合金, P (r t )表示Li原子在t 时刻、 占据晶格位置r
对于式(1)的求解, 一般应经傅里叶变换后在倒空间中进行:
d ? Ρ ( k , t ) d t = c 0 ( 1 - c 0 ) Κ B Τ ? L ( k ) ? { ? F ? Ρ ( r , t ) } k + ξ ( k , t ) ? ? ? ( 2 )
式中 k
? Ρ ( k , t ) ? ? L ( k ) ? { ? F ? Ρ ( r , t ) } k ? ξ ( k , t ) 为与晶格位置坐标k F 可由平均场理论得到, 其中包括原子间有效作用能。
对铝锂合金而言, 其基体为面心立方结构, 亚稳析出相δ ′为L12 有序结构。 为减小计算量, 将三维空间在二维平面上投影。 经二维投影变换后, 式(2)变为
d ? Ρ ( k ′ , t ) d t = c 0 ( 1 - c 0 ) Κ B Τ ? L ( k ′ ) ? { ? V ( k ′ ) ? Ρ ( k ′ , t ) + Κ B Τ [ ln ( Ρ ( r , t ) 1 - Ρ ( r , t ) ) ] k ′ } + ξ ( k ′ , t ) ? ? ? ( 3 )
式中
? V ( k ) = 4 W 1 ( cos π h ? cos π k + cos π h ? cos π l + cos π k ? cos π l ) + 2 W 2 ( cos 2 π h + cos 2 π k + cos 2 π l ) + ? ?
式中 W 1 和W 2 分别为最近邻、 次近邻原子间有效交互作用能。
2模拟基本假设及方法
本文采用模型的基本假设如下:
1) 假设原子扩散是通过直接换位机制, 而非空位机制, 且不考虑位错等晶体缺陷的影响;
2) 假设溶质原子在格点的占位几率随时间的变化率与热力学驱动力成线性关系;
3) 假设动力学系数与原子占位几率无关, 即假设扩散系数与浓度无关;
4) 采用平均场理论计算自由能函数;
5) 由于Al-Li合金δ ′相与基体错配度极小, 因此忽略弹性畸变能对沉淀的影响。
上述假设在定量描述动力学特征时会引起一定的误差, 但不会影响沉淀过程发生的序列, 以及沉淀相的具体结构。 这也是该模型的一个重要的优点。 并且该模型无需预先设定新相结构与沉淀类型, 限定和分割沉淀阶段, 可同时描述原子簇聚和有序化过程, 获得清晰的原子图像演化过程。 为了研究最近邻原子间的相互作用能对δ ′相析出的影响, 本文再假设在外加势场的作用下, 仅最近邻原子间相互自由能发生变化, 而次近邻原子间相互自由能不发生变化。
根据文献
[
11 ]
的报道, 在192 ℃下, 二元Al-Li合金最近邻原子间相互作用能(W 1 )为64.70×10-18 J。 本文选取63.41×10-18 、 63.73×10-18 、 64.06×10-18 、 64.38×10-18 、 64.70×10-18 、 65.03×10-18 、 65.35×10-18 、 65.67×10-18 和66.00×10-18 J, 共计算了9个不同W 1 值的数据点。 模拟以Al-12%Li(摩尔分数)合金为对象, 在128×128格点上, 迭代计算10 000步, 随机噪声的施加强度和时间完全相同。
3模拟结果与分析
Al-12%Li合金在不同W 1 作用下的原子图像演化分别示于图1～图3。 溶质Li原子在各格点处出现的几率由灰度表示, 灰度越大表明Li原子出现的几率越大。 图像中黑白相间的形貌即为有序相, 灰度反差越大说明有序程度越大。 图1为W 1 取64.70×10-18 J时原子图像的演化过程。 在最初时刻, 合金处于完全无序状态, Li原子在各点的出现几率相同, 因此图像中各个格点的灰度相同。 随着迭代步数的增加, 从第900步至1500步, 观察到尺寸不等、 随机分布、 相互独立的有序相核心, 如图1(a)所示。 此后先后形成的有序相晶核相继长大。 图像通过灰度变化显示, 尤其是在有序相附近区域, 无序基体中的Li原子浓度显著下降, 如图1(c)所示。
图2为W 1 取63.41×10-18 J时原子图像的演化过程。 与W 1 取64.70×10-18 J时相比, 有序相形核过程类似, 但是迭代步数达到1 250步时才形核, 而且至10 000步时, 仅形成了一个稳定存在的晶核。
图3为W 1 取66.00×10-18 J时原子图像的演化过程。 可以发现此时有序相的形成过程与前述两种情况有较明显的区别。 首先, 仅迭代750步有序相便已经形成。 其次, 在首先形成的有序相逐渐长大的同时发现有新的有序相形成, 如图3(b)中A 处。 最后, 在个别部位有序相的形成是不完全独立的, 类似于小范围的单相有序畴, 相互毗邻的单相有序畴之间由反相畴界隔开, 如图3(b)中B 、 C 处, 随迭代步数的增加, 反相畴界发生无序化反应, 最终形成无序相界面。
比较图1、 图2和图3可以发现, 在有序相的形成和长大过程中, 最初形成的有序相内部格点的灰度较小, 灰度反差亦较小, 说明此时相内Li原子浓度不高, 有序度也不高, 随时间延长, 格点灰度及格点间灰度反差均逐渐加大。 另外, 有序相形成早期, 从无序的基体到有序相之间存在一个明显的灰度低的过渡带。 从上述δ ′相形成过程原子图像演化特征的共同点来看, 可以判断Al-12%Li合金在不同W 1 值下δ ′相的形成均属于非经典形核长大机制
[11 ]
。 不同的是, 随W 1 的增大, δ ′相形核数量明显增多, 核心由球形变成椭圆形。 在W 1 取66.00×10-18 J的情况下, 小范围单相有序畴的出现, 有序畴被反相畴界分割及反相畴界的无序化反应等, 表现出一定的失稳分解机制的特征, 说明随W 1 的增大, δ ′相的形成机制有向失稳分解过渡的趋势。
图1 Al-12%Li合金在W1取64.70×10-18 J时δ′相沉淀的原子图像
Fig.1 Computed microstructures for Al-12%Li alloy when W 1 =64.70×10-18 J (a)—900 steps; (b)—1 500 steps; (c)—5 000 steps; (d)—10 000 steps
图2 Al-12%Li合金在W1取63.41×10-18 J时δ′相沉淀的原子图像
Fig.2 Computed microstructures for Al-12%Li alloy when W 1 =63.41×10-18 J (a)—1 250 steps; (b)—2 500 steps; (c)—5 000 steps; (d)—10 000 steps
图3 Al-12%Li合金在W1取66.00×10-18 J时δ′相沉淀的原子图像
Fig.3 Computed microstructures for Al-12%Li alloy when W 1 =66.00×10-18 J (a)—750 steps; (b)—1 500 steps; (c)—5 000 steps; (d)—10 000 steps
图4为W 1 取64.70×10-18 J时, 有序相析出过程中晶核内部局域成分序参数和局域序参数随时间和位置的变化规律。 有序相的局域序参数和局域成分序参数是由溶质原子的占位几率换算得到的。 前者反映了合金晶胞尺度范围的有序化程度, 其取值从0到1, 取0表示完全无序状态, 取1表示完全有序。 后者反映了晶胞尺度内溶质原子簇聚的程度, 其取值从0到0.224, 取0表示不含溶质原子, 取0.224, 即δ ′相在192 ℃时的平衡浓度
[12 ]
。 由图1的分析可知, 迭代计算900步时, 在原子图像中可以观察到明显的有序相, 此时相内已经发生了明显的有序化反应, 见图4(a), 但是此时合金的成分却没有发生明显的起伏, 见图4(b)。 这说明在有序相形成的初期, 原子的簇聚滞后于有序化反应。 本文所计算的其它W 1 取值情况下均有类似的现象。 此时的有序相化学成分与平衡状态下的δ ′相的化学成分(22.4%, 摩尔分数)有很大的差别, 因此这种有序相称之为“非化学计量比有序相”。 这种非化学计量比有序相在整个δ ′相的析出过程中存在的时间极为短暂, 经换算大约为0.9 s。 若要在实际实验中观察到该非化学计量比有序相是非常困难的, 这也正体现了计算机实验的优势。 随迭代步数的增加, 局域序参数不断增大, 并最终达到1, 即达到完全有序化。 局域成分序参数亦增大, 但在本文计算的时间范围内, 始终没有达到δ ′相的平衡成分, 即没有得到“化学计量比有序相”。
图5所示为不同W 1 值下, 长程序参数和平均成分偏离序参数随时间的变化规律。 同局域序参数和局域成分序参数类似, 有序相的长程序参数和成分偏离序参数也是由溶质原子的占位几率换算得到的。 不同的是长程序参数反映了合金整体有序化的程度, 平均成分偏离序参数反映了合金整体溶质原子簇聚的程度。 在不同W 1 值下, 长程序参数和成分偏离序参数均由于模拟初期施加的随机力项ξ (r t )出现一个小的峰值, 在随机力项ξ (r t )撤除后, 两序参数下降。 该峰值的出现是计算中人为造成的。 本文中对于不同W 1 值的情况, 随机力项ξ (r t )的施加强度和施加时间完全相同。 在经过一定的驰豫阶段后, 体系开始析出有序相, 表现为长程序参数和成分偏离序参数随时间的延长逐渐升高, 此时体系开始自发分解。 但是, 在不同W 1 值下, 体系所需的弛豫时间是不同的。 可以看出随W 1 值的增大, 弛豫时间缩短。
图4W1取64.70×10-18 J时晶核内局域序参数和局域成分序参数的变化
Fig.4 Local order parameter and local composition parameter profiles across a nucleus while W 1 =64.70×10-18 J
(a)—Local order parameter; (b)—Local composition parameter; 1—900 steps; 2—1 500 steps; 3—5 000 steps; 4—10 000 steps
图5 不同W1值下长程序参数和成分偏离序参数随时间的变化
Fig.5 Variation in order parameter and composition deviation parameter with timesteps at different W 1
(a)—Order parameter; (b)—Composition deviation parameter
图5(a)所示为不同W 1 值下长程序参数随时间的变化规律。 可以看出经过弛豫阶段后, 长程序参数上升, 且在不同阶段表现出明显不同的速度, 说明合金有序化的速度不同。 当W 1 取64.70×10-18 J和66.00×10-18 J时, 长程序参数先快速升高, 而后减缓, 但取66.00×10-18 J时长程序参数增大速度明显较快。 当W 1 值取63.41×10-18 J时, 长程序参数则始终保持比较平缓的增大趋势。 图5(b)所示为不同W 1 值下成分偏离序参数随时间的变化规律。 可以看出经过弛豫阶段后, 成分偏离序参数平稳上升, 说明溶质原子簇聚速度始终相近。 总的来讲, 随W 1 值的增大, 长程序参数和成分偏离序参数增大速度加快, 说明合金有序化反应速度和溶质原子簇聚速度快, 并且在本文所研究的时间范围内, 长程序参数和成分偏离序参数所能达到的最大值亦增大。
图6所示为时间步长为10 000步下长程序参数和成分偏离序参数随W 1 变化曲线。 可以看出,长程序参数和平均成分偏离序参数均随W 1 增大而增大, 且W 1 值在63.41×10-18 ～64.06×10-18 J和65.35×10-18 ～66.00×10-18 J范围内, 两序参数增大幅度较大, 说明在该范围内最近邻原子间相互自由能对合金有序化反应速度和溶质原子簇聚速度影响较大。
图6 10 000步时长程序参数和平均 成分偏离序参数随W1变化的规律
Fig.6 Variation in order parameter and composition deviation parameter with W 1 at 10 000 timesteps
4结论
1) 基于离散格点形式的微扩散方程(Langevin方程)和非平衡自由能函数的计算能比较好地解决求解由高度非平衡、 高度非线性体系中产生随机分布、 多晶胚随时间演化的问题。
2) Al-12%Li合金在不同W 1 值下, 有序相的形成均属于非经典形核长大机制。
3) 随W 1 的增大, δ ′相形核数量明显增多, 孕育期明显缩短。 若W 1 较大, 会出现小面积的单相有序畴、 有序畴被反相畴界分割及反相畴界的无序化反应。 这说明δ ′相的形成机制有向失稳分解过渡的趋势。
4) 在不同W 1 值下, 在有序相形成初期有序化反应先于原子簇聚, 存在非化学计量比有序相, 并逐渐向化学计量比有序相过渡。
5) 在不同W 1 值下, 合金有序化反应速度和溶质原子簇聚速度均增大, 且随迭代步数的增加表现出不同的增大速率。 当W 1 值较小时, 有序化反应速度和溶质原子簇聚速度始终较低; 当W 1 值较大时, 有序化反应速度先快后慢, 而原子簇聚速度变化平缓。
6) 在不同W 1 值下, 合金有序化程度和溶质原子簇聚程度均随时间的延长逐渐增大, 在研究时间范围内, 合金有序化程度和溶质原子簇聚程度所达到的最大值随W 1 值的增大而增大。
参考文献
[1] HafnerJ.Atomic scalecomputationalmaterialsscience[J].ActaMater,2000,48:7192.
[2] 吴兴惠,项金钟.现代材料计算与设计教程[M].北京:电子工业出版社,2002.7382.WUXing hui,XIANGJin zhong.ModernMaterialsComputationandDesign[M].Beijing:ElectronicsIndustryPress,2002.7382.
[3] 哈森P.材料的相变[M].北京:科学出版社,1998.252253.HaasenP.PhaseTransformationsinMaterials[M].Beijing:SciencePress,1998.252253.
[4] HassenP.PhaseTransformationsinMaterials[M].Weinheim:VCH,1991,281.254257.
[5] ChenLQ,KhachaturyanAG.Computersimulationofstructuraltransformationsduringprecipitationofanorderedintermetallicphase[J].ActaMetaMater,1991,39:25332551.
[6] ChenLQ,KhachaturyanAG.Formationofvirtualorderedstatesalongaphase decompositionpath[J].PhysRevB,1991,44:46814684.
[7] ChenLQ,KhachaturyanAG.Kineticsofvirtualphaseformationduringprecipitationoforderedintermetallics[J].PhysRevB,1992,46:58995905.
[8] KupperT,MasbaunN.SimulationofparticlegrowthandOstwaldripeningviatheCahn Hilliardequation[J].ActaMetal,1993,42:18471858.
[9] OkudaH,OsamuraK.ComputersimulationofthekineticsofphasedecompositionwiththeL12typeorderinginanisinglatticesystematlowtemperature[J].ActaMetal,1993,42:13371343.
[10] PoduriR,ChenLQ.Non classicalnucleationtheoryoforderedintermetallicprecipitatesapplicationtotheAl Lialloy[J].ActaMetal,1995,44:42534259.
[11] SchmitzG,HaasenP.DecompositionofanAl Lialloy—theearlystagesobservedbyHREM[J].ActaMetal,1992,40:22092217.
[12] PoduriR,ChenLQ.Computersimulationofthekineticsoforder disorderedandphaseseparationduringprecipitationofδ′(Al3Li)inAl Lialloys[J].ActaMater,1997,45(1):245255.
