锆及锆合金是核能利用的重要材料,具有热中子吸收截面小、热膨胀系数小、热导率高以及与核燃料相容性好等特点,这使得锆合金拥有优异的核性能。所以核工业的发展与锆工业的发展息息相关,锆成为一种重要的战略材料,被誉为“原子时代的第一金属”[1]。基于以上特点,锆及其合金被广泛应用于核反应堆的燃料包壳中, 在反应堆运行期间处于高温高压高辐照的条件下[2],尤其是在长时间服役期间应力应变和辐照缺陷协同作用下,引起微观结构的变化继而导致蠕变、膨胀、硬化、脆化等结构失效,从而导致材料宏观性能的退化[3-5]。自间隙原子是辐照损伤效应之一,目前为止,对于α-Zr的自间隙原子的研究大多停留在无应变条件下。除此之外,对于密排六方(hexagonal closepacked,HCP)晶体结构形式的α-Zr,沿各个晶向原子排列存在差异,决定了α-Zr的自间隙原子迁移行为的各向异性。本文以α-Zr单晶为研究对象,在300~600 K模拟温度下,沿[1120]、[1100]和[0001]3个晶向分别施以1%~3%的应变,利用分子动力学方法系统讨论自间隙原子的迁移行为。
1 缺陷团簇扩散模拟方法本文利用经典的分子动力学软件LAMMPS[6]进行α-Zr的模型建立和在应变条件下缺陷团簇迁移的模拟计算,利用OVITO[7]软件观察缺陷团簇的迁移现象,并计算出缺陷团簇的扩散系数继而拟合得到迁移能。对于α-Zr体系,本文采用嵌入原子势(embedded atom method, EAM)[8]描述体系中各个原子之间的相互作用力。
模拟计算需要设定α-Zr材料体系的晶格参数、弛豫温度、拉伸方向、拉伸变形比例以及时间步长和步数等参数。首先利用LAMMPS软件构建α-Zr体系模型,接着在正则系综(NVT)下进行充分弛豫达到热力学平衡,弛豫温度分别为300、400、500和600 K。图 1为弛豫温度300 K下的α-Zr体系模型,其中,模型A中x轴方向对应[1120]晶向,尺寸为20.89 nm×13.39 nm×12.62 nm,包含147 456个Zr原子;模型B中y轴方向对应[1100]晶向,尺寸为12.93 nm×21.27 nm×12.67 nm,包含145 920个Zr原子;模型C中z轴方向对应[0001]晶向,尺寸为11.64 nm×13.43 nm×21.10 nm,包含153 600个Zr原子。无应变条件下缺陷团簇的扩散模拟本文使用图 1中模型C进行,应变条件下的模拟则是通过分别对这3个模型的x轴、y轴、z轴方向施加拉伸应变后进行,应变比例为1%、2%和3%,α-Zr泊松比为0.35[9]。
![]() |
Download:
|
图 1 热力学平衡状态下的α-Zr模型 Fig. 1 α-Zr model in thermodynamic equilibrium state |
本文利用N.de Diego等[10]的方法在体系中间位置添加4个Zr间隙原子,建立典型的BS型自间隙锆原子以形成缺陷团簇,如图 2所示,在微正则系综(NVE)条件下,模拟计算整个体系在4 ns时间内的微观结构变化情况,同时利用Voronoi cell[11]方法判断体系内的间隙原子并跟踪间隙缺陷质心的运动情况。
![]() |
Download:
|
图 2 α-Zr体系添加自间隙原子构建缺陷团簇 Fig. 2 Defect clusters constructed by adding self interstitial atoms in α-Zr system |
图 2(a)为平面示意图,图 2(b)是体系达到稳定后间隙原子的排布情况,其中较深的原子代表处于完美晶体位置上的锆原子,较浅的原子代表间隙原子。在图 2(a)中,较大的原子是间隙原子,其他较小的原子对应着图 2(b)中较深的原子。
2 缺陷团簇扩散分子动力学研究结果 2.1 缺陷团簇的扩散现象在300 K温度条件下,在无应变和3个晶向([1120]、[1100]和[0001]晶向)均施加3%应变的各情况下,缺陷均呈现出团簇结构。观察缺陷团簇的质心坐标轨迹,结果如图 3所示,其中图 3(a)~(e)是无应变,时间点分别为0、1、2、3、4 ns,模型竖直方向对应着[0001]晶向,发现缺陷团簇在系统中心呈小范围的三维扩散;对于拉伸应变方向为[1120]晶向,如图 3(f)~(j)所示,时间点相同,模型竖直方向对应着[1120]晶向,发现缺陷团簇的扩散径迹是沿[1120]晶向的一维扩散。
![]() |
Download:
|
图 3 300 K下无应变与沿[1120]晶向施加3%拉伸应变时缺陷团簇扩散运动的对比 Fig. 3 Comparison of diffusion motion of defect clusters under no strain and 3% tensile strain along [1120] direction at 300 K |
为了更直观地观察缺陷团簇质心的运动情况,分析300 K温度条件下无应变与沿[1120]、[1100]和[0001]晶向施加3%拉伸应变下缺陷团簇的质心随时间迁移的径迹,结果如图 4所示。
![]() |
Download:
|
图 4 300 K温度下在无应变与沿各方向施加3%的应变条件下缺陷团簇的质心迁移径迹 Fig. 4 Centroid migration path of defect clusters at 300 K without strain and with 3% strain applied in three directions |
图 4中x、y、z坐标方向分别代表[1120]、[1100]和[0001]晶向,点(0,0,0)为模拟系统的中心点,长度单位为Å。无应变情况下,如图 4(a),缺陷团簇呈位于系统中心小范围的三维无规则扩散;当3%的拉伸应变作用在[1120]晶向时,如图 4(b),缺陷团簇的扩散呈现出沿[1120]晶向的一维扩散,扩散径迹沿x方向的坐标变化范围为-40~77Å;当3%的拉伸应变作用在[1100]晶向时,如图 4(c),缺陷团簇的扩散呈现出(0110)晶面上的二维扩散;当3%的拉伸应变的作用在[0001]晶向时,如图 4(d),缺陷团簇的扩散呈现出沿[0001]晶向进行一维扩散,但是范围极小。
根据图 4可知,在α-Zr中发生相同拉伸应变的情况下,对比[1120]、[1100]和[0001]晶向,缺陷团簇沿[1120]晶向的扩散最明显,这与Whting等[12]使用分子动力学方法模拟α-Zr中自间隙原子迁移行为得到的结论相似。Samolyuk等[13]提出在α-Zr中间隙原子的一维扩散倾向沿HCP晶体基面的密排方向,即[1120]晶向。
2.2 缺陷团簇的扩散系数与迁移能利用扩散系数与温度的相关性能够计算缺陷团簇的迁移能。其中,扩散系数的计算公式为:
$ {D_N}\left( T \right) = \frac{{\overline {{R^2}} \left( T \right)}}{{2{n_d}t}} $ | (1) |
式中:R2 (T)是缺陷团簇质心坐标的均方位移;nd是迁移运动的维数;t是模拟时长。为了获得较为准确的扩散系数,又不至于使计算时间过长,本文将总模拟时长平均分成了160份,计算每份的均方位移并结合式(1)考虑平均值继而得到扩散系数:
$ {D_N}\left( T \right) = \frac{1}{{160}}\sum\limits_{i = 1}^{160} {\frac{{\overline {R_i^2} }}{{2nt/160}}} $ |
扩散系数与迁移能和温度之间的函数关系为:
$ {D_N}\left( T \right) = {D_{0, N}}{\rm{exp}}\left( { - \frac{{{E_{\rm{m}}}}}{{{k_{\rm{B}}}T}}} \right) $ | (2) |
式中:Em是迁移能;kB是玻尔兹曼常数;T是温度;D0, N是扩散率预系数。可以看出温度越高,原子动能越大,扩散系数随之呈指数增加。在高温下发生的与扩散有关的过程,温度是最重要的影响因素[14]。
为了研究缺陷团簇的扩散现象,在300~600 K温度区间与各晶向施加拉伸应变情况下,得到的扩散系数如图 5所示。在相同拉伸应变情况下,随着温度的上升,扩散系数呈现出上升趋势,表明原子的热运动能力随着温度的上升而增强,缺陷团簇的扩散能力增强;在相同温度情况下,拉伸应变作用在[1120]与[1100]晶向时,如图 5(a)和图 5(b),扩散系数随着应变量的增加而上升,这说明沿着[1120]与[1100]晶向的拉伸应变促进了缺陷团簇的扩散,对于拉伸应变作用在[0001]晶向的情况,如图 5(c),缺陷的扩散系数随拉伸应变量的增加无明显变化规律。
![]() |
Download:
|
图 5 300~600 K模拟温度下沿各方向施加0~3%应变条件下缺陷团簇的扩散系数 Fig. 5 Diffusion coefficient of defect clusters at 300~600 K simulation temperature with 0~3% strain applied in three directions |
迁移能即扩散激活能,与扩散系数息息相关,是扩散系数与温度的函数。式(2)指出,迁移能的数值即为扩散系数以e为底取得的对数与1/kBT构成的线性关系中斜率的绝对值。本文取300、400、500和600 K共4个温度点分别对在0应变和沿[1120]、[1100]和[0001]晶向施加1%~3%的单向拉伸应变模拟计算得到的缺陷团簇在α-Zr单晶模型中的扩散系数进行拟合,迁移能的拟合结果如表 1所示。
![]() |
表 1 迁移能与预扩散系数 Table 1 Migration energy and pre diffusion coefficient |
无应变条件下,α-Zr中缺陷团簇的迁移能为0.044 eV。当沿[1120]晶向对α-Zr系统模型施加应变时,1%~3%的应变均使缺陷团簇的迁移能降低,且不同应变下数值接近,在0.026~0.028 eV。拉伸应变方向为[1100]时,应变使得缺陷的迁移能降低,应变量为3%时的迁移能相比2%的情况小幅上升可能是由于拟合带来的误差导致。而应变方向为[0001]的情况与前2种不同,随着应变的增加,迁移能随之增加,当应变量达到3%时的迁移能相比无应变时增加了超过一倍,由无应变情况下的0.044 eV增加至0.110 eV。
2.3 缺陷团簇的形成与形成能按图 2所示的方法向模型中添加自间隙锆原子后,使系统在NVT系综下保持添加原子前的温度充分弛豫并使模型达到能量稳定的过程中,每个原子都因为热运动发生相对位置的变化,模拟中后添加的锆原子与模型中原有的原子相互挤压,系统内部逐渐形成较稳定的缺陷团簇,在温度为300 K时,观察到形成的缺陷团簇如图 6所示。
![]() |
Download:
|
图 6 建立4个BS型间隙原子后在300 K弛豫后系统达到热平衡后形成的缺陷团簇 Fig. 6 Defect clusters formed after the system reached thermal equilibrium after 300 K relaxation after 4 BS type interstitial atoms were established |
本文计算了缺陷团簇在无应变与沿[1120]、[1100]和[0001]这3个相互正交的晶向拉伸应变1%~3%的情况下的形成能。其中应变作用在[1120]与[1100]时,缺陷的形成能小于零,当应变量达到3%时,形成能甚至达到了-86 eV,这是不应当的,并会在后续的研究中加以修正。0~3%的拉伸应变作用在[0001]晶向时α-Zr中缺陷团簇的形成能如图 7所示。无应变时,缺陷的形成能为9.299 eV,这与N.de DIEGO等[10]计算得到的8.98 eV相差不多,随着应变量的上升,缺陷团簇的形成能逐渐下降并呈线性关系,拟合得到的形成能与沿[0001]晶向应变量的线性关系为:
$ {E_f} = - 0.673x + 9.345 $ | (3) |
![]() |
Download:
|
图 7 沿[0001]晶向0~3%拉伸应变情况下α-Zr中缺陷团簇的形成能 Fig. 7 Formation energy of defect clusters in α-Zr under 0~3% tensile strain along [0001] direction |
为了研究缺陷团簇在α-Zr中扩散方向的差异性,本文计算了在模拟温度为300 K,分别沿[1120]、[1100]和[0001]施加1%~3%的拉伸应变情况下,缺陷团簇沿着3个晶向的扩散系数,计算结果如图 8所示。
![]() |
Download:
|
图 8 300 K施加0~3%应变情况下缺陷团簇沿各晶向的扩散系数 Fig. 8 Diffusion coefficient of defect clusters along each crystal orientation under 0~3% strain at 300 K |
由图 8可见,在不施加应变的情况下,缺陷团簇沿[1120]和[0001]晶向的扩散系数接近,略高于沿[1100]晶向的扩散系数,扩散系数比值为0.357∶0.256∶0.387;当沿[1120]晶向发生拉伸应变时,随着拉伸应变量的增加,缺陷沿[1120]晶向的扩散系数逐渐增大,沿[1100]和[0001]晶向的扩散系数则逐渐减小,当拉伸应变量达到3%时,缺陷沿3个方向的扩散系数比值为0.965∶0.011∶0.024,沿[1120]晶向的扩散系数相比另外2个方向的扩散系数高了近2个数量级,说明在这种情况下缺陷的扩散主要沿[1120]晶向进行;当沿[1100]晶向发生拉伸应变时,随着拉伸应变量的增加,缺陷在[1120]和[1100]晶向的扩散系数有上升的趋势,[0001]晶向的扩散系数降低,当应变量达到3%时,沿3个方向的扩散系数比值为0.188∶0.797∶0.015,虽然沿[1100]晶向的扩散系数占全部扩散系数的比值很高,但仍不可忽略沿[1120]晶向的扩散;当模型沿[0001]晶向发生应变时,应变量小于2%时,沿3个方向的扩散系数接近且无明显变化规律,当应变量达到3%时,[1120]、[1100]和[0001]晶向的扩散系数均减小,且沿[0001]晶向的扩散系数明显高于另外2个方向,沿3个方向的扩散系数比值为0.243∶0.241∶0.516。
综上所述,无应变下,缺陷团簇在α-Zr中扩散为三维扩散;当发生拉伸应变时,[1120]晶向的拉伸应变对扩散的影响最大且为一维扩散,[1100]晶向的拉伸应变对扩散的影响次之,[0001]晶向的拉伸应变对扩散的影响最小。这些现象与α-Zr的晶格结构相关,由于α-Zr的HCP晶体结构,导致辐照缺陷在3个方向上扩散的差异。
3 结论1) 无应变下缺陷团簇的扩散现象是三维的,且扩散范围较小;[1120]晶向的拉伸应变对扩散的影响最明显且为一维扩散,[1100]和[0001]晶向的拉伸应变对扩散的影响不明显,原因是HCP晶体中[1120]晶向相比[1100]和[0001]晶向的原子排列更加密致。
2)[1120]和[1100]晶向的拉伸应变会减小缺陷团簇的迁移能,而[0001]晶向的拉伸应变会增加缺陷团簇的迁移能。
3) 缺陷团簇的形成能随着[0001]晶向的拉伸应变的增大而减小。当拉伸应变方向为另外2个晶向时本研究尚未得到较科学的缺陷形成能数值。
以上结论能够为Zr材料的辐照损伤实验研究及性能改性提供理论参考。此外关于α-Zr缺陷团簇的扩散机制研究,还需要在间隙原子的数量和应变方式等方面进行改进和补充。
[1] |
刘建章. 核结构材料[M]. 北京: 化学工业出版社, 2007.
( ![]() |
[2] |
KURTZ R J, HEINISCH H L. The effects of grain boundary structure on binding of He in Fe[J]. Journal of nuclear materials, 2004, 329-333: 1199-1203. DOI:10.1016/j.jnucmat.2004.04.262 ( ![]() |
[3] |
CHOI Y, INOUE H. Crystallographic texture change of pilgered zirconium alloy tubes with heat treatment[J]. Materials science forum, 2007, 558-559: 1379-1382. DOI:10.4028/0-87849-443-x.1379 ( ![]() |
[4] |
GAUMÉ M, BALDO P, MOMPIOU F, et al. In-situ observation of an irradiation creep deformation mechanism in zirconium alloys[J]. Scripta materialia, 2018, 154: 87-91. DOI:10.1016/j.scriptamat.2018.05.030 ( ![]() |
[5] |
CHOI S I, KIM J H. Radiation-induced dislocation and growth behavior of zirconium and zirconium alloys-a review[J]. Nuclear engineering and technology, 2013, 45(3): 385-392. DOI:10.5516/NET.07.2013.035 ( ![]() |
[6] |
PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of computational physics, 1995, 117(1): 1-19. DOI:10.1006/jcph.1995.1039 ( ![]() |
[7] |
STUKOWSKI A. Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool[J]. Modelling and simulation in materials science and engineering, 2010, 18(1): 015012. DOI:10.1088/0965-0393/18/1/015012 ( ![]() |
[8] |
DE DIEGO N, SERRA A, BACON D J, et al. On the structure and mobility of point defect clusters in alpha-zirconium: a comparison for two interatomic potential models[J]. Modelling and simulation in materials science and engineering, 2011, 19(3): 035003. DOI:10.1088/0965-0393/19/3/035003 ( ![]() |
[9] |
BRYUKHANOV A A, MANZHIKOV A V, TARASOV A F. Estimating Poisson's-ratio anisotropy for sheet zirconium[J]. Soviet atomic energy, 1988, 64(1): 80-81. DOI:10.1007/BF01124010 ( ![]() |
[10] |
DE DIEGO N, OSETSKY Y N, BACON D J. Mobility of interstitial clusters in alpha-zirconium[J]. Metallurgical and materials transactions A, 2002, 33(3): 783-789. DOI:10.1007/s11661-002-0145-y ( ![]() |
[11] |
DU Qiang, FABER V, GUNZBURGER M. Centroidal Voronoi tessellations: applications and algorithms[J]. SIAM review, 1999, 41(4): 637-676. DOI:10.1137/S0036144599352836 ( ![]() |
[12] |
WHITING B J, BACON D J. Computer simulation of the migration of self-interstitial atoms in alpha-zirconium[J]. MRS online proceedings library, 1996, 439: 389. DOI:10.1557/PROC-439-389 ( ![]() |
[13] |
SAMOLYUK G D, BARASHEV A V, GOLUBOV S I, et al. Analysis of the anisotropy of point defect diffusion in hcp Zr[J]. Acta materialia, 2014, 78: 173-180. DOI:10.1016/j.actamat.2014.06.024 ( ![]() |
[14] |
胡赓祥, 蔡珣, 戎咏华. 材料科学基础[M]. 3版. 上海: 上海交通大学出版社, 2010. HU Gengxiang, CAI Xun, RONG Yonghua. Fundamentals of materials science[M]. 3rd ed. Shanghai: Shanghai Jiao Tong University Press, 2010. ( ![]() |