当前位置：文档库 › 材料设计发展计算

材料设计发展计算

文章编号:100520299(2000)0420057206

材料设计的发展新趋势———材料设计计算方法

冯武锋,王春青,张　磊

(哈尔滨工业大学现代焊接生产技术国家重点实验室,黑龙江哈尔滨150001)

摘　要:计算机计算技术和量子理论等领域的发展,使材料信息数据系统、热力学方法、分子轨道量子理论

乃至神经网络方法、拓扑几何学都应用于材料设计中.采用计算方法进行材料设计具有经济、高效及预测功

能的特点,该方法将成为材料设计发展的重要新趋势.

关键词:材料设计;计算方法;分子轨道

中图分类号:TG409 文献标识码:A

C alculation during material design,a ne w tendency in material design

FEN G Wu2feng,WAN G Chun2qing,ZHAN G Lei

(National K ey Lab of Advanced Welding Production Technology,Harbin Institute of Technology,Harbin150001,China)

Abstract:Presents the application of database system of material information,thermodynamic method, quantum theory about molecular orbital,NN method and topology in material design with the development of computer technology and quantum theory as an important tendency in material design.

K ey w ords:material design;calculation method;molecular orbital

通过大量的实验,寻找规律性、确定最佳性能的材料组分是研究人员常用的材料设计方法.但这种方法极其费时、费力并带有很大的经济消耗.基于科学发展对材料设计的高效、经济并富有预测性的要求,通过计算方法进行材料设计成为必然的趋势.诚如1956年Slater所言[1]:“我不明白为什么材料学家们如此忙于通过实验确定多元合金的组分、合金的结构和相图.根据原子序数我们可以知道元素的原子结构,又有量子力学定律及电子计算机,可以很快计算出相关的结果…….”

但是,评价体系总能量,即使是在绝对零点也很难.一个孤立原子的能量大约在106eV.而每个原子的内聚能在110eV数量级.所以必须有一种方法能够精确到10-6以上[2](10-9范围内),并具有高速运算能力的计算机硬件基础.基于上述

收稿日期:1999-07-21

作者简介:冯武峰(1974-),男,博士生.原因,各种计算设计方法直到近几年才有较大的发展,本文将分类简述.

1　半经验模式材料设计计算方法

自1996年起,日本工程制造中心着手建立L PF(Linus Pauling File)数据库.该数据库覆盖了合金、金属间化合物、陶瓷、矿物等全部无机物的有关结构、衍射、组成及本征性质等信息.大约10万种相关出版物的实验确定或计算所得的数据、至少20万个有关结构、衍射、固有性质的数据项以及3.5万个自1900年记录的结构数据项[1]都将存储在该知识信息系统中.以此高精确度的强大材料数据库为基础,寻找规律.从原理上讲,最终只要输入原子序数以及考虑到可能的化合物组成,将能够实现对材料性能的预测.

L PF研究人员认为,基于L PF这样的强大数据库,可以建立一个知识-信息体系,通过计算有效地预测、开发新材料.

第8卷　第4期材　料　科　学　与　工　艺 Vol.8　No.4 2000年12月 MA TERIAL SCIENCE&TECHNOLO GY Dec.　2000

2　非经验材料设计计算方法

2.1　热力学辅助材料设计

北京科技大学以Redlich 2K ister 多项式及化

学计量复合模型确定了Ag -Sn 二元相模型,采用由Lukas 等人开发的软件进行了系统优化[3],计算所得的合金相图与热力学实验数据相符.对二元系统的理论分析和实验研究,也可以延拓到三元合金领域.芬兰T M K orhonen 根据热力学第二定律,通过最小吉布斯能计算进行了Sn -In -Ag 相平衡计算.由二元延拓到三元的三相组织热力学模型的Muggianu 表达式[4]为

G =

∑3

i =1

x i

+R T

∑3

i =1

x i

ln (x i

)

∑2

i =1∑3

j =i +1

x i

∑2

n =0

(x i -x j )n

+x 1x 2x 3L 123(1)

式中:G 0i 为纯i 元素吉布斯自由能;n 为相状态;

x i 、x j 为i 、j 的摩尔分数;L n

ij 为n 相中i 、j 相互作

用能;L 123为三组分相互作用能.通过计算确定了可在接近Sn/Pb 钎料共晶温度范围进行再流钎焊的钎料合金,其成分w (B )/%为:3.5Ag 、20In 、76.5Sn ,实验证实计算结果正确.

瑞典皇家技术学院开发了Thermo 2Cacl 热化学软件,Seung W ook Y oon 等人将其用于四元合金的相平衡计算中,进行无Pb 钎料Sn -Bi -In -Zn 系统的热力学辅助合金设计及性能评价[5].用扫描差分测热(DSC )的结果与热力学计算的相转变及熔化温度值相比较,二者结论很一致.最终确定了Sn -Bi -5In -6Zn 合金作为Sn/Pb 替代钎料.2.2　分子轨道理论计算方法材料设计

分子轨道理论是量子力学的一种近似计算法,能弥补价键理论的不足[6],准确描述分子内部结构.该方法计算复杂.其中,从头计算法和DV -X α方法较多地应用于材料科学领域.2.2.1　从头计算法材料设计

从头计算方法是在分子轨道理论三个近似的基础上,仅利用普朗克常数、电子静止质量和电量三个基本物理常数而不借助任何经验参数,求解多电子体系的量子理论全电子计算方法.原则上,只要合适地选择基函数,自洽迭代的次数足够多,就一定能够得到接近自洽场极限的任意精确解[7].从头计算法大大优于半经验的计算方法,可以得到各类体系的电子运动状态及其有关的微观信息,能合理地解释或预测原子间的键合、分子的结构、化学反应的过程、物质的性质.但由于计算体系的复杂化,从头计算法的计算量与占用计

算机存储空间急剧增加,不断地对计算技术提出新课题.另一方面,从头计算方法的误差主要来自于非相对论近似和轨道近似.目前主要采用组态相互作用方法校正这些误差,但这种方法显著增加了对计算机内存和机时的需求,对较大分子难于实现.

2.2.1.1　材料电子结构及力学性能、热性能计算

美国G erbrand Ceder 教授应用从头计算法设计固态聚合物电解电池锂阴极材料[8],计算、评价其它过渡金属氧化物作为锂阴极的潜力.通过预测这些新化合物的结构和形成能,很快就能评价出其电池电压及夹层容量,不经实验验证就可选择合适材料.

日本Y oshitada Isono 建立材料原子簇模型,变化键长和键角,采用从头计算法评价立方氮化硼的力学性能、热性能,找出晶格常数与总能量的关系曲线(含金刚石),如图1所示.取能量最小值确定键长、键角,分别为0.156nm 和107.5

°.而实验测得的键长、键角分别为0.156nm 和109.47°.在此基础上计算了立方氮化硼的弹性模量和剪切模量及热传导、膨胀性能.

图1　总能量与晶格常数(a )的关系[9]

Fig.1　Relationship between general energy and lattice con 2

stant [9]

美国加利福尼亚大学进行了空位及其对材料性能影响的电子结构计算.以掺杂-空穴的间隙为函数,对Al 中诸多掺杂元素,应用从头计算自洽电子结构方法计算了总能量.将线性圆球分区近似、轨道方法、原子球型近似以及区域密度近似一起应用于一个含32个原子的超晶胞几何模型中[10].研究发现,在第二近邻晶格距离位置,掺杂与空穴的相互作用对电迁移和扩散性能有十分重要的影响.图2给出了Al 中不同掺杂元素掺杂-空穴作用能与其相距位置的关系.Al 基合金掺杂

?85?材　料　科　学　与　工　艺第8卷

与空穴相互作用的有关扩散趋势的计算结果与实际观察的结果基本一致

.图2　掺杂-空穴作用能E A 与位置的关系

[10]

Fig.2　Relationship between action energy and distance of

adulterant and cavity [10]

G erbrand Ceder 以从头计算法预测复杂氧化

物的结构与性能,在材料合成之前不需要输入实验数据就可知道结构特性并评价其性能,已成功地预测了简单材料的相稳定性[11],并对氧化物晶界结构和缺陷能进行预测[12].研究简单结构晶粒边界与基体缺陷能的差异、导致缺陷在边界上选择偏移的机理,并确定边界电子结构.2.2.1.2　表面、界面结构计算

K im F Ferris 进行了TiO 2的(110)、

(100)面上小分子吸附体的电子结构计算.保持电中性选择分子簇模型,该模型中包含了直接原子和第一近邻原子的相互作用.计算结果表明,在TiO 2的(110)、(100)面上水比甲醇易被吸附[13],其吸附

难易取决于表面上质子输运反应能量阈值的大

小,阈值愈小,愈容易吸附.如(110)面上对水吸附的能量阈值为6.28kJ /mol ,对甲醇的吸附阈值为11.72kJ /mol ,水比甲醇易吸附.但该模型中忽略了更长范围的相互库仑作用,且有关能量阈值的解释只是一种定性认识.

法国A Est ève 应用从头计算法计算总能量,研究了热氧化第一阶段在Si (110)表面氧的吸附.利用最小总能量确定了不同氧、硅原子比时Si (110)表面上氧的最优、次优吸附位置[14],如图

3所示.并利用Monte Carlo 方法模拟了表面的形成及氧化膜生长过程.其中,总能量计算有Super 2cell 、C W 及CPA 三种常用方法[15].

1998年美国真空科学技术杂志报导了有关Si (111)面上单层Sb 不同结构的计算.应用从头计算法确定了Si (111)表面单分子层Sb 的表面原

子结构.由计算得出具有相同最小能量的四种结构,进而推断出在Si (111)表面有不同结构的Sb 共存[16],其中有一种新结构为实验所没有发现的.

图3　不同O ,Si 原子比的(110)Si 面截面图[14]

Fig.3　Si (110)section of various proportion of O and Si

atoms [14]

2.2.1.3　合金相图预测

根据振动自由能经从头计算可以精确预测合金相图.理论及实验都表明[17],合金系统模型中引入振动自由能对精确定量预测合金相图很重要.该方法既是给定当量合金种类结构的弹性模量量子力学计算、计算振动频率的晶格动力学技术,也是评价或预测合金相图的热力学工具,并可综合从头计算法和统计热力学、考虑局部化学环境的作用精确计算合金相图[18].2.2.2　DV -X α分子轨道法材料设计

X α计算方法是Slater 在1951年首先提出

的[19],将电子交换作用能用其统计平均来近似,代入Hartree -Fock 方程,得到X α方程.DV -X α

方法基于Hartree -Fock -Slater 近似,即使是对较大尺寸的原子簇系统也能提供相当精确的电子结构.此方法中,Slater 局域电子交换作用势能包含了电子间的相互交换作用.V XC 由下式[20]给出:

V XC =-3α[3ρ(r )/8

π]1/3

(2)

式中:ρ(r )为局域电子密度;参数α为定值0.7.计算应用了自洽电子近似,采用离散采样方法对矩阵元素积分.哈密尔敦函数矩阵元素及重叠积分都由随机采样计算得出.分子轨道由原子轨道的线性组合而成.对合金元素而言,采用的原子轨道要么是1s -n p (Li 、Be ,n =2;Na ,n =3;K 、Ca 、Ti 、V 、Cr 、Mn 、Fe 、Co 、Ni 、Cu 、Sn 、G a 、G e ,n =4;Y 、Zr 、Nb 、Mo 、Ag 、Cd 、In 、Sn 、Sb ,n

5;La 、

Ce 、Pr 、Nd 、Sm 、Gd 、Tb 、Dy ,n =6),要么是1s -n d (Al 、Si ,n =3)[21]

.与从头计算法相比较,DV -X α是计算方法的改进.其精确度一般略逊于从头计算方法,但因其采用了交换势的定域密度泛函近似,故不必计算大量的多中心积分,计算量仅为从头算法的1/100,是计算含重原子分子和固体原子簇的常用

95?　第4期冯武锋,等:材料设计的发展新趋势———材料设计计算方法

方法之一.

2.2.2.1　合金机械性能预测及合金设计

1994年R Ninomiya在基于电子结构的分子轨道计算基础上提出了预测镁合金机械性能的量化方法.在这个预测中引入一个新参数M k(合金元素的S轨道能级)描述合金的性能.图4为原子簇MMg18模型.在六方系Mg的晶格中,每一个Mg原子周围有12个最近邻Mg原子和6个次近邻Mg原子.原子簇中心的Mg原子被不同的合金元素M所替代.晶格参数为a=0.3203

nm,c =0.5199nm.图5为ΔM k平均值(ΔM k)与合金硬度的关系.图6为ΔM k与合金抗拉强度的关系.此方法在多组分Mg合金机械性能预测及高性能Mg合金设计中极其方便、有效.其中对合金ΔM

的定义[21]为:(x i为i元素摩尔分数)

ΔM

k =∑

x i|M k

-M k

|.(3)

图4　MMg18原子簇模型[21]

Fig.4　MMg18atom cluster model[21]

图5　ΔM k与合金硬度的关系[21]

Fig.5　Relationship betweenΔM k and hardness of alloy[21]

图6　ΔM k与合金抗拉强度的关系[21]

Fig.6　Relationship betweenΔM k and tensile strength of al2 loy[21]

M Morinaga等人用DV-Xα方法研究了Fe、Ni、Ti合金中元素的电子结构,得到20种合金元素的费米能级参数[22]:(1)S轨道能级;(2)键级;(3)离子度;(4)电子状态密度.以同样的方法预测了Al合金的机械性能[23],并用DV-Xα分子轨道理论建立合金设计方法[24],用M k的平均值反映合金机械性能.用此参数预测合金的机械性能与用ΔM k并没有本质的区别.

2.2.2.2　合金化对材料结构性能的影响计算

日本京都大学用DV-Xα分子轨道模型法计算氧化物表面电子结构,研究了ZnO(1010)的电子结构,研究发现该表面有特征价态,从静电势能形成角度讨论了表面电子状态[25].它的状态密度与实验所得到的U PS谱吻合得很好.

DV-Xα法计算还可用于分析合金化对Mg2Ni金属间氢化物电子结构的影响.簇模型为Mg2Ni1-x M x[26],应用重叠高斯函数计算了电子状态密度,计算得到Mg2Ni氢化物的能级结构、空间电子密度分布,由此确定了簇中的原子间键级及每一原子的离子度.研究表明,在纯Mg2Ni 氢化物中Ni-Mg原子间相互作用是最主要的,而且,H原子与Ni原子的键接要强于Mg原子的结合.因此,Ni-H及Ni-Mg原子间相互作用都直接影响了氢化物的相稳定性.当有合金元素M 替代Ni原子位置时,合金化使相互作用有很大变化.即使增加了M-H原子间相互作用,M-Mg 原子间的相互作用降低仍会引起氢化物的不稳定,导致实验中所发现的高生成焓.

日本神户制钢所材料研究实验室进行了二元铜合金电阻增值及固溶强化预测的初步研究,应用DV-Xα法进行了二元铜合金电阻增值和固溶

?材　料　科　学　与　工　艺第8卷

强化内聚能的理论计算.计算所得电阻增值与实验值一致.研究发现,纯铜簇与含有一个固溶原子的簇之间内聚能的差值与错合参数相关[27],而错合参数又与固溶强化有关.这一结论表明,一种不经实验数据就可以预测电阻增值和固溶强化的可行性方法.

DV-Xα法可以与从头算法结合应用于材料性能预测中.神谷列夫等人[28]应用DV-Xα法基函数进行Ab-initio Band计算,得到晶体性能.利用最小二乘法及径向分布函数计算得到了紧缩高斯基函数.基于不同的晶格参数及原子坐标进行总能量计算,确定了弹性常数和介电常数.

2.2.2.3　Ba TiO3与电极金属的DV-Xα法润湿性分析

DV-Xα簇方法还可以研究Ag或Al

Ba TiO3(100)表面吸附的润湿机理及化学键状态.Ba TiO3表面计算采用了Madelung势能函数,引入[TiO5]6簇,一个点电荷置于Ba原子位置.从计算所得能级结构中知道,在Ag或Al与Ba2 TiO3界面HOMO能级升高[29],这一升高和Ba2 TiO3与电极金属的润湿性能有关.

2.2.2.4　变成分比合金的电子结构拟合计算

哈尔滨工业大学对钎料合金采用DV-Xα法分子轨道计算,探讨了簇模型的选择原则[30],并以不同成分比的Au/Cu合金为研究对象,采用DV-Xα量子学计算方法分析了由合金化导致的晶格常数变化对电子状态密度的影响,发现晶格常数变化前后的态密度十分相近.在此基础上,提出变成分比合金计算的有序拟合、无序拟合及宽域有序拟合、狭域有序拟合方法.研究结果表明,对任意成分比例的Au/Cu合金,选取相应的权重因子,采用有序拟合、尤其是采用合适的狭域拟合,可以得到精确的合金电子结构.

2.2.3　神经网络与统计学方法联合材料设计

神经网络与统计学方法结合可以进行材料结构-性能描述.基于不同的模型方法得到神经网络与统计学联合的非线性连续回归方法(NLCR)进行材料的结构-性能预测.该方法包括[31]:一般最小二乘回归(OL S)、偏微分二乘回归(PL S)、组分回归(PCR)、脊回归(RR),以及一些非线性方法,如:单隐含层BP网、投影跟踪回归(PPR)及非线性偏微分最小二乘回归(NL PL S).NLCR 方法与OFBNN方法[32]相似,不同之处在于前者基于线性映射、后者基于非线性局域映射.

2.2.4　拓扑结构学材料设计

Linn W Hobbs教授将拓扑数学方法与计算算法应用于网状玻璃结构的可行性的研究中.目前正在开发完全或部分相连的四面体网状物[33],如SiO2、SiC、Si3N4、P2O5.

3　结束语

基于理论计算的材料设计计算方法具有省时、高效、经济的特点,尤其是对材料结构、性能的定量预测功能使其必将成为材料设计方法的主要发展趋势.而计算方法的难度在于物理模型的准确建立、量子化学方程求解的复杂性和精确度及符合实际情况的较大簇模型的计算量.基于合理的近似用高速计算机解薛丁谔方程已经成为可能.单电子势能极大地简化了多体相互作用.随着计算科学的发展、交叉学科的渗透,尽管有些问题十分复杂,最终都有可能逐步得到解决.

参考文献:

[1]V ILLARS P,ONODERA N,IWA TA S.The linus

pauling file(L PF)and its application to material de2

sign[J].Alloys and Compounds,1998,279:127. [2]CHEL IKOWSKY J R.Encydopedia of applied physics

[M].[s.l.]:[s.n.],1991.

[3]XIE Y,Q IAO Z.Thermodynamic reoptimization of

the Ag-Sn system[J].Phase Equilibria,1996,17

(3):2082250.

[4]KORHONON T M,KIV ILAHTI J K.Thermody2

namics of the Sn-In-Ag solder system[J].Elec2

tronic Materials,1998,27(3):1492158.

[5]Y OON S W,SON J R,L EE H M,et al.Thermo2

dynamics2aided alloy design and evaluation of Pb2free

solder,Sn-Bi-In-Zn system[J].Acta Mater,

1997,45(3):9512960.

[6]路琼华,朱裕贞,苏小云.工科无机化学[M].上

海:华东化工学院出版社,1998.

[7]廖沐真,吴国昆,刘洪霖.量子化学从头计算方法

[M].北京:清华大学出版社,1984.

[8]Ceder G.The design of a Li-cathode materials for

solid polymer electrolyte batteries[EB/OL].http∥

https://www.wendangku.net/doc/d718504927.html,/papers/data.html,1997.

[9]ISONO Y,KISHIMO TO H,TANA K A T.Evalua2

tion of mechanical and thermal properties of cubic

boron nitride by ab2initio calculation[J].Materials

Science Research International,1998,4(1):39244.

[10]STEME P A,EK J V,HOWELL R H.Electronic

structure calculations of vacancies and their influence

on materials properties[J].Computational Materials

Science,1998,10:3062313.

[11]CEDER G.Ab2initio prediction of structure and prop2

　第4期冯武锋,等:材料设计的发展新趋势———材料设计计算方法

erties of complex oxides[EB/OL].http://burgaz.

https://www.wendangku.net/doc/d718504927.html,/papers/data.html,1997.

[12]CEDER G.Prediction of structure and defect energies

at grain boundaries in oxides[EB/OL].http//

https://www.wendangku.net/doc/d718504927.html,/papers/data.html,1997.

[13]FERRIS K F,WAN G L Q.Electronic structure cal2

culations of small molecule adsorbates on(110)and

(100)TiO2[J].J Vac Sci Technol,1998,A16(3):

179021793.

[14]ESTΙV E A,ROU HAN I M D,ESTΙV E D.First

principle calculation of oxygen adsorption on a(100)

silicon surface2first stages of the thermal oxidation

[J].Computational Materials Science,1998,10:942

98.

[15]TA KEUCHI N.First2principle calculations on Si

(111)[J].J Vac Sci Technol,1998,A16(3):

179021793.

[16]ABRIKOSOVL A,RUBAN A V,JOHANSON B,et

al.Total energy calculations of random alloys:super2 cell,connolly2williams and CPA methods[J].Compu2 tational Materials Science,1998,10:3022305. [17]Ceder G.Accurate predictation of alloy phase dia2

grams:the contribution of the vibrational free energy

[EB/OL].http//https://www.wendangku.net/doc/d718504927.html,/papers/data.

html,1997.

[18]COLON ET C,PASTUREL A.Ab2initio calculation

of thermodynamic data an phase diagram of binary

transition metal based alloys[J].Phase Equilibria,

1994,15(3):3302337.

[19]唐敖庆.量子化学[M].北京:科学出版社,1982.

[20]李俊清,河天敬,王　俭,等.物质结构导论[M].

合肥:中国科学技术大学出版社,1990.

[21]NOMIY A R,YU K AWA H,MORINA G A M,et al.

An electronic approach to the prediction of the me2 chanical properties of magnesium alloys[J].Alloys

and Compounds,1994,215:3152323.

[22]MORINA G A M,NASU S,J SAITO H A,et al.Al2

loying effect on the electronic structure of aluminum

[J].J Phys Condes Matter,1991,3:681726828. [23]MORINA G A M,K AMADO S.An electronic ap2

proach to the prediction of the mechanical properties

of aluminum alloys[J].Modeling Simul Mater Sci

Eng,1993,1:1512164.

[24]MORINA G A M,N INOMIG A R.An approach to al2

loy design with the aid of a molecular orbital theory

[J].Trans Mat S oc Jpn,1994,A16:1952199. [25]Ohba F,Matsunaga K,Tanaka I,et al.Electronic

structure calculation of ZnO surface[A].Papers of DV

-XαSeminar in Japan[C].[s.l.]:[s.n],1996.832

87.

[26]TA K AHASHI Y,YU K AWA H,MORINA G A M.

Alloying effects on the electronic structure of Mg2Ni

intermetallic hydride[J].Alloys and Compounds,

1996,242:982107.

[27]ITSUMI Y.A preliminary study for prediction of

electrical resistivities and solid solution hardening of

binary Cu alloys using DV-XαMethod[J].Papers

of DV-XαSeminar in Japan,1996.68271.

[28]ITSUMI Y,SA K AI E,DAIMON M.Calculation of

crystal properties based on Ab2initio based band calcu2 lation using DV-Xαbases[A].Papers of DV-XαSeminar in Japan[C].[s.l.]:[s.n.],1996.2902 293.

[29]SEKIN E R.Analysis for wetting properties of Ba TiO3

with electrode metals by DV-Xα[A].Papers of DV

-XαSeminar in Japan[C].[s.l.]:[s.n.],1996.

1072110.

[30]WAN G C Q,FEN G W F.Selecting cluster model in

alloy designing calculation with discrete variational Xαmethod[J].1999Acta Metallurgica Sinica(English

Letters),2000,13(1):84288.

[31]BA KSHI B R,CHA TTERJ EE R.Unification of neu2

ral and statistical methods as applied to materials struc2

ture2property mapping[J].Alloys and Compounds,

1998,279:39246.

[32]CHEN C L P,CAO Y,L ECLAIR S R.Materials

structure2property prediction using a self2architecting

neural network[J].Alloys and Compounds,1998,

279:30238.

[33]HOBBS L V.Topology of glass structure[EB/OL].

http:∥https://www.wendangku.net/doc/d718504927.html,/papers/data.html,1997.

(责任编辑　吕雪梅)

?材　料　科　学　与　工　艺第8卷

材料设计发展 计算

材料设计发展计算