网络首发时间: 2016-06-22 08:23

稀有金属 2016,40(10),1002-1007 DOI:10.13373/j.cnki.cjrm.xy15020803

高压高温条件下铍弹性性能数值模拟

李婷婷 何力军 李国龙 梁森 许德美 钟景明

宁夏大学宁夏光伏材料重点实验室

西北稀有金属材料研究院稀有金属特种材料国家重点实验室

摘 要：

基于密度泛函理论(DFT)的平面波赝势方法对密排六方结构(hcp)铍(α-Be)的弹性性能进行了第一性原理计算研究。利用Material Studio软件的CASTEP模块并采用广义梯度近似(GGA)完成计算。分别得到了铍单晶体的弹性常数,铍多晶体的体积模量、剪切模量、杨氏模量、泊松比等参数对外加静水压变化(0~100 GPa)和对环境温度(0~1300 K)变化的变化规律。还利用第一性原理数值模拟方法对体心立方结构(bcc)的铍(β-Be)的弹性性能进行了计算。计算结果与文献实验数据符合度高。结果显示:各弹性参数值随压力的增长而单调增加,升温过程的影响作用则相反;描述铍晶体变形耦合性的C_(12)和C_(13)对压力或温度变化比其他弹性参数更为敏感;加压过程使铍晶格结构参数c/a值增加并趋向理想密排值,而在温升过程中这个值则持续降低;加压过程是一个原子密堆度提高的过程,据此对铍高压hcp-fcc(面心立方结构)相变的可能性和实现形式进行了讨论。

关键词：

铍;密度泛函;弹性性能;高压;高温;

中图分类号： O614.21;O641.1

作者简介：李婷婷(1990-),女,宁夏银川人,硕士研究生,研究方向:稀有金属及材料力学性能;E-mail:lttkuaile-2009@163.com;;何力军,教授;电话:13895008252;E-mail:helijun2003@sohu.com;

收稿日期：2015-02-10

基金：国家自然科学基金项目(11062010,51474177);宁夏自然科学基金项目(NZ14038)资助;

Numerical Analysis on Elastic Properties of Beryllium under High Pressure or High Temperature

Li Tingting He Lijun Li Guolong Liang Sen Xu Demei Zhong Jingming

Ningxia Key Laboratory of Photovoltaic Materials,Ningxia University

State Key Laboratory of Special Rare Metal Materials,Northwest Rare Metal Materials Research Institute

Abstract：

Based on plane wave pseudo-potential in the framework of density functional theory( DFT),first principle calculations were performed to investigate elastic properties of beryllium for hexagonal close packed( hcp) structure( α-Be). The generalized gradient approximation( GGA) was adopted in concrete calculation in using CASTEP code of Material Studio program. The elastic constants of monocrystalline beryllium,the bulk modulus,shear modulus,Young's modulus and Poisson's ratio of polycrystalline beryllium were calculated under the condition of hydrostatic pressure in the range of 0 to 100 GPa. The same work mentioned above was also finished under the condition of environment temperature varying from 0 to 1300 K. First principle numerical simulations were also performed to investigate elastic properties of beryllium for body centered cubic( bcc) structure( β-Be). Calculation results were in good agreement with those experimental data obtained from references. Results showed that the elastic constants increased monotonically with pressure increasing,but decreased with temperature increasing monotonically; C_(12) and C_(13),describing deformation coupling along different directions,were more sensitive than the other elastic parameters to the changes of pressure or temperature; pressure increase made the value of beryllium lattice structure parameter c/a increase and approach to the ideal close packed value,while in temperature increasing process,this value deceased steadily. Pressure increase was a process of raising the density of atom packing,based on this understanding,the possibility and actualizing form of hcp-fcc( face centered cubic) phase transformation under high pressure were discussed.

Keyword：

beryllium; density functional theory; elastic properties; high pressure; high temperature;

Received： 2015-02-10

金属铍(Be)的原子序数为4,相对原子量为9,位于元素周期表的第二主族(ⅡA),在地壳中的含量极少,现实环境中以hcp(密排六方)结构,所谓α-Be存在。由于其优良的物理力学性能:低密度(1.85 g·cm^-3),高刚度(303 GPa),高导热系数(200 W·m^-1·K^-1),高比热(1825 J·kg^-1·K^-1),低热膨胀率(11.4×10^-6K^-1),高熔点(1558 K),高热中子散射截面(6.9 bam)等,被应用于包括航空航天、核反应堆、X射线、仪表材料、光学镜体材料等在内的诸多重要工程领域 ^[1] 。

在铍研究领域,随着基础理论的完善和计算方法的进步,数值模拟日渐成为一种对实验研究进行指导和补充的有效研究手段,如对金属铍在非极端条件下力学行为的有限元模拟 ^[2] 。另一方面,铍往往也在极端条件的环境中服役,因此其在高压或高温条件下的力学行为成为目前铍研究领域的热点之一。本文利用第一性原理方法,分别计算了hcp结构铍晶体(单晶、多晶)在高压、高温环境下的弹性参数,并预测了相变体bcc(体心立方)结构铍晶体(单晶、多晶)的弹性常数、弹性模量、泊松比等力学指标。

1 计算方法

基于Material Studio软件的CASTEP模块完成计算,具体计算方案为:对铍分别设置一系列不同晶格常数的原胞模型,并计算与之相对应体系的总体能量E和原胞体积V,利用Birch-Murnaghan状态方程(EOS)对E-V进行函数拟合,以能量最低状态为平衡状态,对应原胞体积和尺寸为平衡体积和优化晶格参数,并以此状态各参数作为后续弹性性能计算的基础,这是第一性原理计算最核心的部分。利用C_ij=(?²E/?ε_i?ε_j)/V得到弹性系数C_ij(式中ε_i为应变分量)。再基于弹性常数由Voigt-Reuss-Hill统计平均方法,获得多晶体的各弹性常数:体积模量B、剪切模量G、杨氏模量E和泊松比ν。

主要计算设置:采用密度泛函理论的赝势平面波方法及周期性边界条件,电子波函数通过平面基组展开,由超软赝势(USP)来描述离子实与价电子之间的相互作用势,采用BFGS算法对结构进行几何优化,广义梯度近似交换相关函数为GGA-PBE形式;对hcp铍体材料的截断能设为360 e V,倒空间中布里渊区积分采用Monkhorst-Pack方法,K点网格设置为17×17×10,能量收敛精度为1.0×10^-6,力精度为0.01 ^[3] 。

2 计算结果

2.1 hcp结构弹性特征

弹性常数C_ij代表36个分量,由于对称性,独立分量最多可以有21个。如前所述,通常状态下铍属于hcp结构(即α-Be),有5个独立的弹性常数即C₁₁,C₁₂,C₁₃,C₃₃,C₄₄,加上非独立的C₆₆(C₆₆=(C₁₁-C₁₂)/2),6个参数决定铍的全部弹性行为和特征。

2.2 高压条件下的弹性性能

按照计算方法,首先得到平衡状态下的晶格参数a,c,零压体积弹性模量B₀,以及零压体弹模量对压强的一阶偏导数B'₀,结果见表1。可以看出,计算结果合理。

对模型施加静水压,计算在此条件下的铍弹性性能,结果如图1。由图1可知,随着压力的增加,计算的所有弹性参数C_ij,M_i(即B,G,E),B/G,ν均单调增加。C₁₁,C₃₃的绝对值在C_ij中增长最多,但基于对各自在0 K时的弹性常数C_ij0的比较,由图1(b)可知,C₁₂,C₁₃的相对值增长最为显著。这说明,随着压力的增加,铍不同方向弹性行为的耦合度在增长,这一点可以通过泊松比ν的增长得到印证,见图1(f),众所周知,泊松比是描述不同方向弹性行为耦合度的物理量。在图1(c)中给出了B的计算值与实验数据 ^[7] 比对,无论数值还是趋势,符合度还是令人满意的。需要说明的是,在材料高压力学行为研究中,往往进行的是材料体积V变化对外加压力p响应的所谓p-V曲线测量实验,基于文献 [ 7] 给出的是p-V数据,利用B=-(dp/d V)/V计算得到各压力下的体积模量B。结合图1(a,b),还能看出C_ij的变化还呈现出分组的特征,即C₁₁与C₃₃,C₁₂与C₁₃,即C₄₄与C₆₆分别具有相似的变趋势,尤其是C₄₄与C₆₆,在图1(b)中图线几乎完全重合。

表1 α-Be晶体原胞晶格参数、体弹性模量及模量一阶导数Table 1 Calculated and experimental crystal parameters,bulk modulus and pressure derivatives of modu-lus ofα-Be  下载原图

表1 α-Be晶体原胞晶格参数、体弹性模量及模量一阶导数Table 1 Calculated and experimental crystal parameters,bulk modulus and pressure derivatives of modu-lus ofα-Be

图1 计算铍弹性性能随压力变化情况Fig.1 Calculated elastic parameters of beryllium versus pressure

(a)Elastic constants C_ij;(b)Relative values of elastic constant C_ij/C_ij0,C_ij0elastic constant at 0 GPa;(c)Elastic modulus M_i,compared with experimental data of Ref.[7];(d)Relative values of elastic modulus M_i/M_i0,M_i0elastic modulus at 0 GPa;(e)B/G;(f)Possion's ratioν

由图1(e)可知,压力外载可以提高铍的B/G。B/G是一个与材料塑形潜力密切相关的指标,其实质反映了材料在弹性状态下体积改变与形状改变的竞争关系,其值越大,则G相对于B越小,形状较之体积改变越容易发生,表现出的塑性越好。塑性性能优越的金属如铅、金、银等的B/G值在通常条件下(室温、标准大气压)分别达到5.19,6.18和3.46 ^[8] 。铍则仅在0.8左右,是常见金属中是最低的,与此对应的,铍具有强的脆性,室温延伸率一般不超过5%。

2.3 高温条件下的弹性性能

本节对各温度下铍弹性性能的计算,计算依据文献 [ 9] 中给出的在不同温度下测得的铍的晶格常数c和a。计算结果如图2。

与加压过程相反,总体而言,随着温度的升高,本文计算的所有弹性参数C_ij,M_i(即B,G,E),B/G,ν均单调降低。由图2(a)看,C_ij与实验值 ^[10] 符合程度较高,造成差别的主要原因在于计算的依据参数 ^[9] 与实验数据 ^[10] 对应的材料参数并不相同。还可看到C_ij分组变化特征还是很明显。与0 K时的弹性常数C_ij0比较的相对值变化仍以C₁₂,C₁₃的表现更为突出(图2(b))。这说明,铍不同方向弹性行为的耦合度也对温度变化产生敏感。图2(c)给出了弹性模量B,G,E随温度的变化及与实验数据的对比,计算数值和趋势与测量值符合度较高。图2(d)是在升温过程中各个弹性模量相对于0 K时弹性模量M_i0的变化,对比图1(d)可知,温度引起的模量值相对变化不仅明显小于压力作用效应,而且各模量的相对变化趋势具有很高的相似性。图2(c,d)还表明,铍的弹性模量在高温下衰减有限,显示出优越的高温力学性能。

图2(f)给出的是关于泊松比ν对温度的变化情况。本文计算表明ν是关于温度的减函数。比较相关的实验数据,文献 [ 11] 给出的在20~300 K间ν也是随温度升高而下降,但其分布数值区间远小于本文计算值;文献 [ 12] 给出在300~900 K间ν的分布数值区间与本文计算值相符,但其表现为温度增函数的特征。由式(10)知,ν的数值可由B,G的测量值确定。由于B与G数值接近且均随温度升高而变化平缓,因此,二者测量值的些许浮动,都会引起引起ν发生显著的改变,其数值范围仍处于室温下金属铍泊松比的常见范围0.02~0.08内。

图2 计算铍弹性性能随温度变化情况Fig.2 Calculated elastic parameters of beryllium versus temperature

(a)Elastic constants C_ij,compared with experimental data of Ref.[10];(b)Relative values of elastic constants C_ij/C_ij0,C_ij0the elastic constants at 0 K;(c)Elastic modulus M_i,compared with experimental data of Ref.[12];(d)Relative values of elastic modulus M_i/M_i0,M_i0the elastic modulus at 0 K;(e)B/G;(f)Possion's ratioν,compared with experimental data

升温引起的B/G变化按本文的计算是下降的,但变化幅度很小,见图2(e)。另一方面,温度升高确实会引起铍金属由脆性向塑性的转变(升温至500 K以上,延伸率可提高一个数量级),显然,这种塑性性能改善与晶格结构的改变(即B/G的变化)关系薄弱,起关键性作用的因素主要在温度改变对铍内滑移系的影响。事实上,在室温条件下,基面滑移即(0001) 滑移占绝对优势地位,当温度升至500 K以上时,启动柱面滑移即(0001) 滑移的临界剪切应力(CRSS)已由室温时是基面CRSS的4~5倍降至与其相当 ^[13] 。于是,大量柱面滑移参与变形,铍的塑性性能显著提高。虽然此时其形变可以超过30%,但这种塑性并不是真正的或者说三维的塑性过程,基面和柱面滑移均属于面内滑移即a滑移,离面滑移即c滑移始终未能启动。这一缺陷也原则性地限制了铍塑性的改善。导致c滑移不能启动的原因在于铍特殊的晶格结构和价电子状态特征,这些因素发生效应的一个综合体现就是极低的B/G值。分析表明,提高环境温度,并不能改善B/G值,进而从根本上改善铍的塑性性能。

对温度因素的考虑体现在将该温度下的晶格常数作为计算依据,但通过CASTEP模块计算出的各弹性参数实质为对应晶格常数下温度在0 K的弹性参数(这是该软件的固有设定),至于温度对金属晶格振动、自由价电子运动本身的影响,并未涉及。由固体物理相关理论可知 ^[14] ,温度对自由价电子动能的影响甚小,金属自由价电子动能主要源于费米能,对铍而言,固态时其因温度产生电子动能变化较费米能低4个数量级。与实测结果符合程度较高的计算结果表明:对金属弹性行为的影响力可能主要来自自由价电子;在具体计算中将温度变化的影响仅以晶格尺寸变化作为体现是可以接受的。

2.4 温度效应、压力效应对弹性性能的影响对比

在以上分别分析的基础上,本节通过建立B/G与原子体积的关系以及泊松比ν与c/a的关系,对温度和压力对铍弹性性能的影响进行对比。总体而言,压力效应与温度效应是相反的。加压和升温均会改变晶格尺寸,进而原子体积发生改变,图3(a)综合二者影响给出了统一的关于B/G对原子体积变化的依赖关系,可知B/G是原子体积的减函数,图线简单,基本呈直线型。原子体积的减小能有效提高B/G值,这一过程与加压过程是相一致的,但升温发生的是反向效应。原子体积的减小,实质是原子密堆度的增加。越来越紧密的原子结构会使得金属继续被压缩的难度持续增长,因而体积模量会持续、显著地增加。相比而言,紧密的原子排列一方面也会导致晶体发生形状改变得困难,另一方面由此引起的持续减小的原子间距又有利于晶体形状改变的发生。两种因素竞争的结果,使得G值的增长远不及B值剧烈。因而B/G值在此原子体积减小过程中会持续增长。

加压和升温过程也会改变晶格形状,即改变c/a值,见图3(b)。总的看,泊松比ν与c/a成增函数关系,而且图线简单,基本呈直线型。升温导致c/a值降低,已经实测验证 ^[9] ;加压可引起c/a值的升高,这不单符合的本文计算结果,也有实验数据支持 ^[7] 。因此,升温过程降低了铍的泊松比,加压的效果正好相反。

在通常条件下,温度升至1534 K左右时,铍会发生相变,由hcp结构(α-Be)变成bcc结构(β-Be)。近年来关于铍的高压相变也是研究热点 ^[15,16] ,这些基于第一性原理计算的推测是高压导致的铍相变仍然是hcp到bcc结构的转变,但尚未得到任何实验研究的支持或否定。不过从加压与升温效果的区别看,hcp→bcc的高压相变推测并不成立,因为bcc的c/a值是1.414,且是一种较hcp结构更为松散的密堆结构。参考图3,升温过程可以降低c/a值,增大原子体积,降低原子堆积致密度,因此高温引起hcp→bcc相变。但加压效应恰与升温相反,此过程会引起升高c/a值,并使得原子堆积致密化,与bcc结构特征不符。hcp→fcc相变倒是与加压效应相协调的一种可能:fcc比铍通常的hcp结构更致密,fcc的c/a值是理想值即1.633。因此,由bcc向fcc相变,完全符合加压推进的结构转变方向。当然,这种可能性也如前述对hcp→bcc相变的推测,尚无实验结果予以支持或否定。

由于hcp结构与fcc结构在基面原子排布的相似性,因此hcp→fcc转变可能是逐步完成的,随着压力的增加,fcc结构层错变得容易发生,于是fcc通过层错逐步增多的方式逐渐实现hcp→fcc的转变,而非通常的一次性完成。

3 发生hcp-bcc相变后的弹性性能

bcc结构铍,即β-Be,自然条件下只能在距离铍熔点1560 K非常近的一个狭窄温度范围内存在,但因其可能存在的优于α-Be的塑性潜力而被关注 ^[17] 。bcc结构的独立弹性常数只有3个:C₁₁,C₁₂,C₄₄。计算中,除铍体材料的截断能为600 e V,K点网格设置为30×30×30外,其余程序及参数选择均同于前述对hcp结构的计算设置。在计算中温度分别设定为0和1534 K,前者利用软件优化功能获得匹配晶格参数为2.491×10^-10m,后者为β-Be存在条件下的实测晶格参数2.549×10^-10m ^[17] 。结果见表2。

图3 温度、压力效应对铍弹性性能影响的对比Fig.3 Comparison of temperature effect and pressure effect on beryllium elastics properties

(a)B/G versus volume of atom;(b)νversus c/a

与前节计算结果相比,bcc相铍的B/G值、及泊松比ν值明显大于hcp相的值。bcc相的G值仅为hcp相一半,这使其位错启动阻力显著减小。可以推测,bcc相铍的塑性性能确实可能大幅优于hcp相铍。如能有有效方法将该结构相固定保持到室温条件,将为铍材料的应用带来革命性推动。

表2 β-Be的弹性常数及泊松比Table 2 Elastic constants ofβ-Be and Possion's ratio  下载原图

表2 β-Be的弹性常数及泊松比Table 2 Elastic constants ofβ-Be and Possion's ratio

4 结论

总体而论,加压与温升过程对铍弹性性能的影响作用是相反的,加压将导致铍的所有弹性参数值上升及几何参数降低,增温则正好相反。对比计算结果与文献实验数据可知,将温度变化以晶格参数变化作为计算依据,对弹性性能的计算效果还是可接受的。与升温相比,加压效应对铍的塑性影响更为基础。对加压也能导致hcp→bcc相变的假说提出质疑,认为加压引起hcp→fcc相变更具合理性,并推测了一种实现的可能方式。各弹性参数的计算结果与文献实测值有较好符合,给出了其在温度或压力影响下的变化规律和趋势,具有理论和工程上参考价值。

参考文献

[1] Zhong J M.Research of Micro-Yield Behavior and Mechanism of Beryllium Metal[D].Changsha:Zhongnan University,2001.1.(钟景明.金属铍的微屈服行为及机理研究[D].长沙:中南大学,2001.1.)

[2] Wan S,He L J,Zhang W,Xu D M,Zhong J M,Wang Z H.Numerical analysis of micro elastic mismatch of polycrystalline beryllium[J].Chinese Journal of Rare Metals,2015,39(8):764.(万朔,何力军,张伟,许德美,钟景明,王战宏.多晶铍微观弹性失配行为的数值分析[J].稀有金属,2015,39(8):764.)

[3] Sholl D S,Steckel J A.Density Functional Theory[M].Hoboken:John Wiley&Sons,Inc,2009.49.

[4] Mackay K J H,Hill N A.Lattice parameter and hardness measurements on high purity beryllium[J].Journal of Nuclear Materials,1963,8(2):263.

[5] Timmerhaus K D,Barber M S.High Pressure Science and Technology(Volume 1)[M].New York:Plenum Press,1979.80.

[6] Benedict L X,Ogitsu T,Trave A,Christine J,Wu Philip A,Sterne,Eric Schwegler.Calculations of highpressure properties of beryllium:construction of a multiphase equation of state[J].Physical Review B,2009,79(6):064106.

[7] Nakano K,Akahama Y,Kawamura H.X-ray diffraction study of Be to megabar pressure[J].Journal of Physics:Condensed Matter,2002,14(44):10569.

[8] Guinan M W,Steinberg D J.Pressure and temperature derivatives of the isotropic polycrystal line shear modulus for 65 elements[J].Journal of Physics and Chemistry of Solids,1974,35(11):1501.

[9] Stonehouse A J,Beaver W W.Physical properties[A].Hauser H H.Beryllium:Its Metallurgy and Properties[C].Berkeley(CA,USA):University of California Press,1965.191.

[10] Goldberg A.Atom,crystal,elastic,thermal,nuclear,and other properties of beryllium[R].UCRL-TR-224850,Livermore(CA,USA):Lawrence Livermore National Laboratory,2006.

[11] Smith M F,Mc Donald J M,Watson R D,Whitley J B.Thermo mechanical testing of beryllium for limiters in ISX-B and JET[J].Fusion Technology,1985,8:1.

[12] Scaffidi-Argentina F,Longhurst G R,Shestakov V,Kawamura H.Beryllium R&D for fusion applications[J].Fusion Engineering and Design,2000,51-52:23.

[13] Martin A J,Ellis G C.The ductility problem in beryllium[A].The Metallurgy of Beryllium[C].London:Chapman&Hall.Ltd.,1963.3.

[14] Cao Q X,Lei T M,Huang Y X,Li Y F.Foundation of Solid State Physics[M].Xi'an:Xidian University Press,2008.160.(曹全喜,雷天明,黄云霞,李云芳.固体物理基础[M].西安:西安电子科技大学出版社,2008.160.)

[15] Guo Z C,Luo F,Cheng Y.Phase transition and thermodynamic properties of beryllium from first-principles calculations[J].Computational Materials Science,2014,84(1):139.

[16] Sin'ko G V,Smirnov N A.Relative stability and elastic properties of hcp,bcc,and fcc beryllium under pressure[J].Physical Review B,2005,71(21):214108.

[17] Aldinger F,Petzow G.Constitution of beryllium and its alloys[A].Donald W,Gilbert J L.Beryllium Science and Technology(Volume 1)[C].New York:Plenum Press,1979.235.