基于修正球形双晶模型的金属Al晶界能分子动力学计算

张明亮，杨 亮，魏承炀，李赛毅, ，张新明,

(1. 中南大学 材料科学与工程学院，长沙 410083；2. 中南大学 有色金属材料科学与工程教育部重点实验室，长沙 410083；3. 中南大学 有色金属先进结构材料与制备协同创新中心，长沙 410083)

作为金属多晶材料内相邻晶粒间的交界区，晶界(GB)是材料显微组织结构的重要组成部分[1]。完整描述晶界的空间取向需要5个参数(或自由度)[2]，其中3个用来确定两相邻晶粒间的取向差，另两个表示晶界面取向。晶界因存在与其空间取向相关的原子错排而具有比晶粒更高的能量，单位面积晶界所具有的附加自由能即为晶界能[3]。晶界能作为晶界的重要特性之一，影响着材料再结晶和晶粒长大等过程中的组织演变，进而影响材料的力学及物理性能[4-10]。因此，有关晶界能的定量表征及其与晶界参数相关性的研究一直是材料科学领域的重要课题。通过实验定量测定晶界能时，所得结果受到材料纯度和晶界参数精度等的严重影响，目前仅获得了极少数的晶界能结果[3,11-14]。READ等[15]基于晶界的位错模型，推导出晶界能与取向差角的关系式，在理论上描述了晶界能受该晶界参数的影响，能够较好地应用于小角度晶界的能量计算。

近年来，随着计算机模拟技术的发展，基于分子动力学的晶界能计算方法逐渐成熟，为从原子尺度探究晶界能及其各向异性提供了有效手段。该类方法的基本思想是：首先根据晶界参数构造具有一定边界条件的双晶体系，然后选取合适的原子势描述原子间的相互作用，并在绝对零度下采用经典牛顿方程对体系进行弛豫。当能量降低至平衡时，晶界所具有的附加能量即为晶界能。用于构造双晶体系的模型主要有考虑周期性边界的块状(Block)模型和考虑非周期性边界的球形(Sphere)模型。块状模型需保证晶界面内两个方向的模型尺寸与晶体原子排布的周期长度相适应，而对于任意取向参数的晶界这点很难或无法满足。所以，该模型主要被应用于一些特殊取向的晶界，如重位因子较低的倾斜晶界和扭转晶界[16-22]。为了克服块状模型的局限性，LEE等[23]提出了球形模型。该模型采用非周期性边界(即自由表面)，可实现任意晶界的晶界能计算。但是，球形模型仅通过不同取向晶粒的刚性对接来确定晶界结构[23]，这种对接可能使晶界处的原子相距太近甚至重叠，从而使所构造的晶界处于非稳态或亚稳态，对应的晶界能将高于理论值。对于给定参数的晶界，其结构和能量都无法预知。因此，需要尝试构造大量初始晶界结构并比较它们的能量，其中能量最低的结构可以当成该晶界参数所描述的晶界结构，基于该结构所计算的附加能量才是晶界能。块状模型[21-22]通过引入原子删除操作来删除晶界处不合理原子，能够成功获取初始晶界结构。已有球形模型则尚未考虑原子删除操作，无法排除不合理初始晶界结构对计算结果的影响。

本文作者在基本球形双晶建模方法的基础上，引入已在块状模型中获得成功应用的原子删除操作来修正原始球形模型，并以金属 Al 〈011〉对称倾斜晶界(STGB)为例，对比分析原始球形模型、修正球形模型和块状模型所算的晶界能，考察原子删除操作对晶界能的影响。同时，将修正模型应用于Al〈001〉扭转晶界(TWGB)和〈011〉非对称倾斜晶界(ATGB)的晶界能计算，以探讨取向差角和晶界面倾斜角对晶界能的影响。最后，基于上述3类晶界的计算结果，考察晶界能与重位因子的相关性。

图1 任意晶界的球形双晶体系构造过程示意图Fig. 1 Schematic diagram showing construction of spherical bicrystal with an arbitrary GB

1 模型与方法

本研究借助 LAMMPS软件构造晶界能的计算模型[24]。图1所示为球形双晶体系构造过程的示意图。首先，在参考坐标系下构建一个直径为D的球形初始单晶(样品Ⅰ)，将初始单晶分别按照取向gA、gB进行旋转分别得到单晶A和B，使它们满足给定的取向差关系 g = gA；然后，沿给定的晶界面方位分别将A和B剖成两个完全相同的半球，保留不同侧的半球形晶粒 A′和 B′；最后，将 A′和 B′沿剖面刚性对接，组成一个球形双晶(样品Ⅱ)。在0 K下对单晶和双晶进行能量最小化弛豫，弛豫后体系的总能量分别为E1和E2。由于单晶与双晶表面原子的排布总体相同，二者所具有的表面能可认为相等[23]。所以，相对于单晶而言，双晶所具有的额外能量为晶界的能量，即晶界能 γ0为

式中：S是晶界面积( S = π D2/4)。

由图1可知，通过半球形晶粒A′和B′沿剖面刚性对接而构造的双晶体系，可能在晶界处存在相距过近或重叠的原子。从物理角度看，这些原子不能存在于实际材料的晶界结构中。块状模型的应用研究[21-22]表明，此类不合理原子的存在将导致晶界结构处于非稳态或亚稳态，从而使晶界能量异常升高，即所得晶界能将高于其理论值。因此，需合理删除这些原子，以获取稳态的晶界结构。

本研究参照块状双晶模型中的原子删除操作[21]来删除不合理的原子，并构造不同的初始晶界结构以获取给定参数晶界的结构和晶界能。该操作具体过程：先设定一临界原子间距(dc)识别晶界处不合理的原子对(即间距小于 dc的相邻两原子)；然后删除原子对中的任意一个来获取晶界结构，再通过改变dc以获取不同的初始晶界结构，最后选定最稳态的晶界结构(对应最低晶界能量的结构)为给定参数晶界的结构。对于具有面心立方结构的金属，晶体理想点阵中原子间的最近邻距离为02 /2a (a0为晶格常数)。为了避免删除体系中间距合理的原子，同时又尽可能多地获取不同的初始晶界结构，本研究在 0.005a0到 0.7a0之间以0.005a0为增量依次选取dc的取值。其中，对于晶界处间距小于dc同时又分别属于晶粒A′和B′的两个原子，选取不同原子进行删除将得到两种不同的晶界结构。因此，修正模型在原理上可构造晶界参数相同的初始晶界结构总数为2×(0.7-0.005)/0.005=278。显然，进行了原子删除操作的双晶和初始单晶的原子总数不再相等，与双晶具有同等原子总数单晶的能量也不再为E1。此时，对应单晶的能量为双晶原子总数N和平均原子势能e1的乘积。因此，利用修正模型所构造晶界的晶界能(dγ)计算公式为

式中：2E′为引入原子删除操作后的双晶的能量。

在计算某给定参数晶界的能量时要对上述278种初始晶界结构所对应的双晶体系进行能量最小化弛豫，并按照式(2)计算所对应的能量。其中，能量最低值被认为是该晶界参数下的晶界能，相应结构为全局最稳态结构，即给定晶界参数晶界的结构。

本研究利用修正模型考察了金属 Al中的〈011〉STGBs、〈011〉ATGBs和〈001〉TWGBs 3 类晶界，这些晶界对应的取向差角(θ)、晶界面倾斜角(Ф)和重位因子(Σ)见表1。所构造双晶模型的直径为D=10 nm。选定二次近邻修正嵌入式原子势[25-26]来描述原子的能量和相互作用力，并利用共轭梯度算法进行体系能量最小化弛豫。

表1 本研究中晶界所对应的取向差角、晶界面倾斜角及重位因子Table 1 List of misorientation angle, inclination angle and coincidence factor of GBs considered in this work

2 结果与讨论

2.1 原子删除对晶界能的影响

图2所示为分别采用原始球形模型、修正球形模型和块状模型计算得到的〈011〉对称倾斜晶界的晶界能(分别用 γ0、γd和 γb表示)随取向差角的变化。由图2可以看出，3种模型所得晶界能随取向差角的变化趋势整体上相同，都在70.53°和129.52°处形成两个显著的极小值，该趋势与实验所测Al的晶界能结果一致[3]。70.53°和 129.52°两晶界刚好对应(111)和(113)孪晶界，晶界处原子错排程度较低，所对应的附加能量也随之较低，即晶界能较低。图2的结果还表明，在16.10°、38.94°和109.47°附近，原始模型和修正模型所得晶界能在绝对值大小和随取向差角的变化趋势上存在明显差异。在这几个晶界处修正模型所得晶界能呈现出明显的极小值，而原始模型的结果却呈现出明显的极大值，且比前者的分别高出40%以上。修正模型所得结果与块状模型结果一致。

图2 采用原始球形模型、修正球形模型和块状模型计算所得〈011〉对称倾斜晶界的晶界能(分别用 γ0、γd和 γb表示)随取向差角的变化Fig. 2 Variation of energy of 〈011〉 STGBs with misorientation angle calculated using original sphere model (γ0),modified sphere model (γd) and block model (γb)

晶界能本质上由原子排列决定，上述原始模型和修正模型在16.10°、38.94°和109.47°所得晶界能的差别可从晶界处原子排布的不同来分析。以 θ=109.47°为例，图3所示为不合理原子被删除前后的晶界原子排布，分别对应于原始模型和修正模型所获得初始晶界结构。如图 3(a)所示，在原始模型构造的晶界处存在一些相距过近的原子(图3(a)中椭圆内的原子)。根据Lennard-Jones两体势[27]可知，在一次近邻范围内，原子势能随原子间距的减小而急剧升高。因此，晶界处间距过小的原子将具有较高的能量，直接导致晶界能升高。由图3(b)可看出，在修正模型构造的双晶体系中，存在于原始模型构造的双晶体系中的不合理原子(即具有较高能量的原子)被删除。因此，与原始模型相比，修正模型所构造的θ=109.47° 晶界具有更稳定的结构和更低的晶界能。

图2结果还表明，修正模型与块状模型所算晶界能在绝对值上仍然存在可见的差别。例如，修正球形模型所算结果相对于块状模型的平均误差和相对平均误差分别为14 mJ/m2和4%。这可能与两个方面有关。首先，与块状模型相比，修正球形模型在计算给定取向参数晶界的能量时，未对相邻两晶体在晶界面内进行刚体平移操作以获取尽可能多的初始晶界结构。与省略原子删除操作类似，省略刚体平移操作也可能导致球形模型计算的晶界能并非全局最小值，从而产生一定的误差。其次，球形模型考虑非周期性边界条件，双晶体系中的自由表面势必会影响晶界结构，导致晶界能的计算误差[23]。

图 3 基于球形模型所构造的 109.47°〈011〉对称倾斜晶界和初始晶界结构对比Fig. 3 Initial configurations generated before (a) and after (b)deletion of unreasonable atoms in sphere model for 109.47°〈011〉 STGB (Misorientation axis is perpendicular to plane of paper; circles and triangles stand for atoms in (011) plane in different layers in crystals; blue and red symbols stand for crystal A′ and B′, respectively)

2.2 晶界能与取向差角、晶界面倾斜角和重位因子的相关性

图4所示为〈001〉扭转晶界的晶界能随取向差角的变化情况。由图 4可知，当取向差角较小时(θ＜28.07°)，晶界能受取向差角的影响较大，随取向差角的增大呈现出升高的趋势，与 Read-Shockley公式所体现的趋势相符；当取向差角继续增大时，晶界能基本保持在一个较小的范围内波动。基于晶界的位错模型及相关理论可知[15]，当取向差角较小时，晶界由离散的位错构成，且位错密度(或原子错排程度)随着取向差角的增大而升高，从而导致晶界能随取向差角的增大而升高。当取向差角增大到一定值时，构成晶界的位错相互重叠，位错密度有所上升。然而，当取向差角进一步增大时，由位错重叠和位错间相互作用导致的位错湮灭将在一定程度上抵消由于取向差角增大所致的位错密度的升高，位错密度将不再发生明显变化，即晶界能不再升高而保持在一个较小的范围内波动。

图 5所示为 Σ3〈011〉非对称倾斜晶界的晶界能随晶界面倾斜角Ф的变化情况。结果显示，晶界能在Ф为0°～64.76°范围内随Ф的增大而逐渐升高，在74.21°附近达到最大值，随后缓慢降低。这种趋势与已有金属 Al的模拟结果[21]以及金属 Cu的相关实验结果[11]相符。Ф为0°和90°分别对应于对称共格孪晶界和对称非共格孪晶界，晶界面两侧的原子关于晶界面成镜面对称排布。与相邻晶界相比，这两个晶界的原子错排程度更低，对应的晶界能也更低。当晶界面偏离上述对称孪晶界时，两侧晶粒在晶界上的原子排布不同(例如 Ф=70.53°时，晶界面在两侧晶粒中的取向分别为{111}和{115})，即呈现非对称排布。此时，晶界处原子错排度和对应晶界能都将升高。

图 4 利用修正球形模型所算〈001〉扭转晶界的晶界能随取向差角的变化Fig. 4 Variation of energy of 〈001〉 TWGBs with misorientation angle, calculated using modified sphere model

图5 利用修正球形模型所算Σ3 (70.53°)〈011〉非对称倾斜晶界的晶界能随晶界面倾斜角的变化Fig. 5 Variation of energy of Σ3 (70.53°) 〈011〉 ATGBs with inclination angle calculated using modified sphere model

重位因子Σ是晶体点阵中重合点阵占体系总点阵比例的倒数，在一定程度上体现晶界结构的特点。早期研究[3,28]认为，Σ值越小，晶界所具有的能量就越低。为了考察晶界能与重位因子的相关性，图6所示为本研究中3类晶界的晶界能及对应的重位因子。由图6可以看出，高重位因子晶界的晶界能并不一定高于低重位因子晶界的，如Σ＜10晶界的能量并非都低于(甚至有时高于)Σ＞10晶界的能量。而且，同一 Σ值下(如Σ3)，不同晶界所具有的能量并不相同，甚至存在很大差异。由此看来，晶界能与Σ之间并不存在特定的相关性，无法通过Σ简单预测或评判晶界能的相对大小。OLMSTED等[22]采用块状模型计算并分析了388种重位点阵晶界的能量与重位因子的关系，得到了与本研究一致的结论。

图 6 利用修正球形模型计算的〈011〉对称倾斜晶界、〈001〉扭转晶界和〈011〉非对称倾斜晶界的晶界能随重位因子(Σ)的变化Fig. 6 Variation of energy of 〈011〉 STGBs, 〈001〉 TWGBs and 〈011〉 ATGBs with coincidence factor Σ according to predictions using modified sphere model

3 结论

1) 与原始球形模型相比，修正模型所构造的〈011〉对称倾斜晶界因删除了不合理的原子而更加稳定，所获得的晶界能及其随取向差角的变化趋势与已有实验和模拟结果相符更好。

2) 对于〈001〉扭转晶界，在取向差角较小时晶界能随取向差角的增大而升高，该趋势与 Read-Shockley晶界能计算公式所预测的结果一致；当取向差角继续增大时，晶界能保持在较小的范围内波动。对于 Σ3〈011〉非对称倾斜晶界，晶界能随晶界面倾斜角的增大而呈现出先升高后减小的趋势，在对称孪晶界处的晶界能低于相邻晶界的。

3) 〈001〉扭转晶界、〈011〉对称倾斜晶界和〈011〉非对称倾斜晶界的晶界能结果表明，晶界能与重位因子间不存在特定的相关性，不能基于重位因子预测晶界能的相对大小。

[1] HU Hui-e, YANG Li, ZHEN Liang, SHAO Wen-zhu, ZHAN Bao-you. Relationship between boundary misorientation angle and true strain during high temperature deformation of 7050 aluminum alloy[J]. Transactions of Nonferrous Metals Society of China, 2008, 18(4): 795-798.

[2] BRANDON D G. The structure of high-angle grain boundaries[J]. Acta Metallurgica, 1966, 14(11): 1479-1484.

[3] HASSON G C, GOUX C. Interfacial energies of tilt boundaries in aluminium: Experimental and theoretical determination[J].Scripta Metallurgica, 1971, 5(10): 889-894.

[4] CHEN Da, LU Min, LIN Dong-liang. Molecular static simulation of energy features of interaction between grain boundary and dislocations in Ni3Al alloy[J]. Transactions of Nonferrous Metals Society of China, 1994, 4(3): 67-70.

[5] 卢 山, 李子然, 昝 祥. 晶界对近片层 TiAl高温动态力学行为影响的数值模拟[J]. 中国有色金属学报, 2012, 22(2):379-387.LU Shan, LI Zi-ran, ZAN Xiang. Numerical simulation of dynamic mechanical behavior of near lamellar TiAl at elevated temperature with influence of grain boundary[J]. The Chinese Journal of Nonferrous Metals, 2012, 22(2): 379-387.

[6] 毛卫民, 陈 垒, 萨丽曼, 余永宁, 李云峰. 晶界对低压电解电容器铝箔腐蚀结构的影响[J]. 中国有色金属学报, 2004,14(1): 1-5.MAO Wei-min, CHEN Lei, SA Li-man, YU Yong-ning, LI Yun-feng. Influence of grain boundaries on corrosion structure of low voltage aluminum foil[J]. The Chinese Journal of Nonferrous Metals, 2004, 14(1): 1-5.

[7] YAN Wen, CHEN Jian, FAN Xin-hui. Effects of grain boundaries on electrical property of copper wires[J].Transactions of Nonferrous Metals Society of China, 2003, 13(5):1075-1079.

[8] 陈健美, 周卓夫, 唐建国, 张新明. 一种考虑晶界能各向异性模拟晶粒长大的Monte Carlo方法[J]. 材料导报, 2003, 17(8):77-79.CHEN Jan-mei, ZHOU Zhuo-fu, TANG Jian-guo, ZHANG Xin-ming. A monte carlo simulation of grain growth with grain boundary energy anisotropy considered[J]. Materials Review,2003, 17(8): 77-79.

[9] 李 旭, 周 庆, 陈明和, 王小芳. 基于晶界迁移率与晶界能各向异性的晶粒生长元胞自动机模拟[J]. 机械工程材料,2012, 36(7): 82-87.LI Xu, ZHOU Qing, CHEN Ming-he, WANG Xiao-fang.Cellular automata simulation for grain growth based on anisotropic grain boundary mobility and grain boundary energy[J]. Materials for Mechanical Engineering, 2012, 36(7):82-87.

[10] 魏承炀. 再结晶织构的相场模拟与机理研究[D]. 广州:华南理工大学, 2012: 49-58.WEI Cheng-yang. Phase field modelling of texture evolution and its mechanisms during recrystallization annealing[D].Guangzhou: South China University of Technology, 2012:49-58.

[11] WOLF U, ERNST F, MUSCHIK T, FINNIS M W. The influence of grain boundary inclination on the structure and energy of Σ3 grain boundaries in copper[J]. Philosophical Magazine A, 1992, 66(6): 991-1016.

[12] GJOSTEIN N A, RHINES F N. Absolute interfacial energies of[001] tilt and twist grain boundaries in copper[J]. Acta Metallurgica, 1959, 7(5): 319-330.

[13] MORI T, MIURA H, TOKITA T, HAJI J, KATO M.Determination of the energies of [001] twist boundaries in Cu with the shape of boundary SiO2particles[J]. Philosophical Magazine Letters, 1988, 58(1): 11-15.

[14] MORI T, ISHII T, KAJIHARA M, KATO M. Determination of the energies of [001] symmetric tilt boundaries in Cu from the shape of boundary SiO2[J]. Philosophical Magazine Letters,1997, 75(6): 367-370.

[15] READ W T, SHOCKLEY W. Dislocation models of crystal grain boundaries[J]. Physical Review, 1950, 78(3): 275-289.

[16] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅰ. Boundaries on the (111) and (100) planes[J].Acta Metallurgica, 1989, 37(7): 1983-1993.

[17] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅱ. Boundaries on the (110) and (113) planes[J].Acta Metallurgica, 1989, 37(10): 2823-2833.

[18] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅲ. Symmetrical tilt boundaries[J]. Acta Metallurgica et Materialia, 1990, 38(5): 781-790.

[19] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅳ. Asymmetrical twist (general) boundaries[J].Acta Metallurgica et Materialia, 1990, 38(5): 791-798.

[20] RITTNER J D, SEIDMAN D N. 〈110〉 symmetric tilt grain-boundary structures in fcc metals with low stacking-fault energies[J]. Physical Review B, 1996, 54(10): 6999-7015.

[21] TSCHOPP M A, MCDOWELL D L. Structures and energies of Σ3 asymmetric tilt grain boundaries in copper and aluminium[J].Philosophical Magazine, 2007, 87(22): 3147-3173.

[22] OLMSTED D L, FOILES S M, HOLM E A. Survey of computed grain boundary properties in face-centered cubic metals: I. Grain boundary energy[J]. Acta Materialia, 2009,57(13): 3694-3703.

[23] LEE B J, CHOI S H. Computation of grain boundary energies[J].Modelling and Simulation in Materials Science and Engineering,2004, 12(4): 621-632.

[24] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics, 1995,117(1): 1-19.

[25] LEE B J, BASKES M I. Second nearest-neighbor modified embedded-atom-method potential[J]. Physical Review B, 2000,62(13): 8564-8567.

[26] LEE B J, SHIM J H, BASKES M I. Semiempirical atomic potentials for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, Al, and Pb based on first and second nearest-neighbor modified embedded atom method[J]. Physical Review B, 2003, 68(14): 144112.

[27] LENNARD-JONES J E. Cohesion[J]. The Proceedings of the Physical Society, 1931, 43(5): 461-482.

[28] BALLUFFI R W, MEHL R F. Grain boundary diffusion mechanisms in metals[J]. Metallurgical Transactions A, 1982,13(12): 2069-2095.