Effects of Nb concentration and temperature on generalized stacking fault energy of Zr–Nb alloys by molecular dynamics simulations

The effects of Nb concentration and temperature on the generalized stacking fault energy (GSFE) of basal, prismatic I, pyramidal I and II plane for Zr-Nb alloys are investigated by molecular dynamics simulations (MD). The stable and unstable SFEs of different slip systems show no significant change with the increasing Nb concentration (0, 0.5, 1.0, 1.5, 2.0, and 2.5 at.%) in Zr-Nb alloys at 0 K. Basal, pyramidal I and II planes slip of Zr-Nb alloys prefer to deform by full dislocation with the temperature increases. Additionally, plastic deformation anisotropy of Zr-Nb alloy is improved with the increasing temperature using both embedded atom method (EAM) and angular-dependent potentials (ADP). The present work provides a theoretical basis for understanding enhanced plasticity of Zr-Nb alloys under finite temperature.


Introduction
Zirconium (Zr) alloys are widely applied as nuclear structural materials because of their excellent anti-neutron irradiation, good high temperature mechanical properties, and outstanding corrosion resistance. However, according to Von Mises criterion, the ductility of hexagonal close-packed (hcp) Zr is severely limited because of the absence of enough independent easy slip systems [1]. The ductility, plastic deformation behavior, dislocation slip and twin of metals are closely related to Generalized stacking fault energy (GSFE) [2][3][4]. Therefore, investigating the GSFE of Zr alloys is beneficial to understand their plastic deformation mechanisms and further develop novel Zr alloys. GSFE was first proposed by Vítek in 1968 [5]. Since then, some significant researches have been done to study the relationship between the GSFE and the plastic deformation mechanisms of bcc, fcc, and hcp metals [6][7][8][9][10][11]. So far, experiments, first-principles calculations, and MD simulations are three main methods to study GSFE, but SFE values measured by experiments are difficult and usually associated with large error [12], whereas the first-principles method is the most widely used in GSFE calculations based on quantum mechanics because of its high accuracy. Recently, first-principles calculations were utilized to explore the solute atom effect on the GSFEs of hcp metals, especially in Mg alloys [13]. For example, Wang et al have proposed a design map of Mg alloys based on GSFEs by first-principles calculations, which provided a guideline for the development of novel Mg alloys [14]. Meanwhile, first-principles calculations were also used to study GSFEs in Zr alloys [15,16]. However, most GSFEs calculated by first-principles were at 0 K, and the temperature effects were ignored, yet temperature should be a crucial factor in material processing, and previous work indicated that the ductility and strength of solid solutions at elevated temperatures mainly rely on the SFEs [17][18][19][20]. What's more, solute atoms are usually located on stacking fault plane, this is not consistent with the random distribution of solute atoms in the metal matrix in experiments, and the number of atoms doesn't exceed 1000 in most cases, so it's difficult to estimate the effect of solute concentration and distribution on GSFEs by first-principles calculations. Molecular dynamics can simulate the properties of materials under finite temperature, and the number of simulated atoms is much greater than that by first-principles. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. contains C, N, Ti, Nb alloying elements at 300 K by MD simulations, the simulation results are in reasonable agreement with some published experimental/calculated data [21]. MD simulations of GSFEs were also carried out to predict dissociation paths of 〈c+a〉 dislocation in pure Mg [22]. The crystal structure of both Mg, Zr are hcp, thus they show similar stacking fault position and structure verified in previous study [10]. In addition, Nb is an important alloying element in Zr based alloys frequently-used as pressure tubes and fuel cladding in nuclear reactors, therefore, it is of great significance to study the Nb concentration and temperature effects on GSFEs of Zr alloys using MD simulations.
In the present work, MD simulations are employed to investigate Nb concentration and temperature effects on the basal, prismatic I, pyramidal I and II planes of GSFEs of Zr-Nb alloys. The stable and unstable SFEs, full and partial dislocation, pyramidal dislocation and deformation characteristics of different slip systems in Zr-Nb alloys are analyzed and discussed in detail.

Simulation method
Open source code large-scale atomic/molecular massively parallel simulator (LAMMPS) is applied to MD simulations [23]. Two interatomic potentials, EAM and ADP are employed [24,25]. EAM potential provides an effective description for the interaction between Zr-Nb binary alloys. ADP is constructed based on a large amount of ab initio data, the potential is appropriate for the simulation of Zr-Nb binary alloys based on α-Zr.
The slip systems, atoms, coordinates, and box length in Zr-Nb single crystal are shown in table 1. The concentration of Nb (at. %) in this work is set to 0, 0.5%, 1.0%, 1.5%, 2.0%, and 2.5% considering that the concentration of Nb is generally less than 2.5% in Zr-Nb alloys used as industrial nuclear materials. Nb atoms are distributed in the Zr matrix randomly. The simulation for all Zr-Nb alloys with different Nb concentration is repeated 3 times to get an average value of stable and unstable SFEs. Zr-2.5 at.% Nb random solid solution model at 0 K is shown in figure 1 as a representative, VESTA software is used to build different slip systems of crystal models [26], and the atomic model visualization software is OVITO [27].
When calculating the GSFEs in Zr-Nb alloy at 0 K, the first step is relaxing the model by conjugated gradient minimization with an energy and force convergence of 10 -10 , and the maximum convergence step is 1000. The box is divided into two blocks as the upper block and the lower block. Pyramidal I, II slip systems allow atoms to relax in both Y and Z axes while other slip systems only allow atoms to relax in the Z axis. Except for {0001}〈10-10〉 slip system, the upper half of the model moves along the Y axis, the other models move along the X axis, and each step moves |b|/20, b is Burges's vector of each slip system, the conjugate gradient algorithm is used to minimize system energy for each step. Then the equation (1) is employed to calculate GSFE curves [5].
Where γ is GSFE, E 0 and E f are the potential energy of the system before and after creating the stacking fault, respectively, A is the area of the stacking fault. When calculating the GSFEs in Zr-Nb alloy at finite temperature, the first step is to equilibrate the model using a Nose-Hoover thermostat and barostat within an NPT ensemble at 300, 400, 500, 600 K (considering the service temperature of Zr-Nb alloy in the Pressure Water Reactor is around 600 K), respectively. The time step and relaxation steps are set to 0.001 ps and 100,000 steps respectively. A vacuum layer with 20 Å thickness is added at both ends of the Z axis after obtaining the finite temperature atomic model. One of the cross-sections of  figure 2 as a representative. Subsequently, the atomic model at finite temperature is divided into upper block and lower block, for each individual slip system, the upper half of the crystal slips a Burges's vector length with respect to the lower half of the crystal, each step slips |b|/20, the conjugate gradient algorithm is applied to minimize system energy for each step with an energy and force convergence of 10 -10 . Finally, GSFE curves are calculated by equation (1).

Results and discussion
3.1. GSFEs of Zr-Nb alloys at 0 K Figure 3 shows five slip systems of GSFEs in Zr-2.5 at.% Nb alloys at 0 K. GSFEs of {0001}〈11-20〉 slip system exist an unstable SFE (γ usf ) at 1/3〈11-20〉, representing the minimum energy barrier to overcome when slipping occurs in this direction [28], and also indicating that it is the lowest energy barrier for dislocation nucleation [29]. Two unstable SFEs (g usf I and g usf II ) and one stable SFE (γ sf or I 2 ) exist in the {0001}〈10-10〉 slip system, in which the stable SFE is the local minimum SFE [30]. As for {10-10}〈11-20〉 slip system, it is the midpoint of the  Burges's vector along the slip direction that is taken as unstable SFE (γ usf ). GSFEs of pyramidal plane {10-11} 〈11-23〉 and {11-22}〈11-23〉 slip systems also exist two unstable SFE (g usf I and g usf II ) and one stable SFE (γ sf ). In this work, we only focus on g usf I because the unstable SFE in the previous literature generally refers to the smaller one [14].
The stable and unstable SFEs in different slip systems for Zr-Nb alloys with different Nb concentration at 0 K are shown in figure 4. It is observed that in the composition range of Nb=0, 0.5, 1.0, 1.5, 2.0, 2.5 at.%, the stable and unstable SFEs in various slip systems for Zr-Nb alloys show no apparent change with the increase of Nb concentration at 0 K, which is different from the results in Ti alloys by first-principles method [31]. The results can be understood as follows: firstly, the concentration of alloying elements in Ti-X alloys is much higher than Nb concentration in Zr alloys; secondly, alloying elements located on stacking fault plane by first-principles calculations while Nb atoms distributed in Zr matrix randomly using MD simulations. Hence, the effects of alloying elements on the stable and unstable SFEs in Ti alloys are much larger than Nb concentration on Zr alloys.

GSFEs of Zr-Nb alloys at finite temperature
The γ sf /g usf I ratio has a significant influence on deformation behavior of metals. If the value of γ sf /g usf I deviates from 1.0, deformation by partial dislocation is preferred, while the γ sf /g usf I is close to 1.0, deformation is more likely to form full dislocations [6,[32][33][34]. To investigate the temperature effect on the deformation behavior of Zr-Nb alloy, we calculated the GSFEs of Zr-Nb alloy at finite temperature. Figure 5 presents the γ sf /g usf I ratio of  One can see that the γ sf /g usf I value of three types slip systems gradually increases as the temperature increases, indicating that basal {0001}〈10-10〉, pyramidal I {10-11}〈11-23〉 and II plane {11-22} 〈11-23〉 tend to slip by full dislocation with the increase of temperature from 0 K to 600 K. In addition, as for {0001}〈10-10〉 slip system, an increase of Nb concentration will decrease the value of γ sf /g usf I at the same temperature, in other words, the increase of Nb concentration will promote {0001}〈10-10〉 slip system deformation by partial dislocation if the temperature keeps constant. However, the increase of Nb concentration does not have much effect on γ sf /g usf I of {10-11}〈11-23〉 and {11-22}〈11-23〉 slip system at the same temperature.
In Mg alloys, ( / c usf P B ) x /( / c usf P B ) Mg parameter is defined to estimate the trend of promoting prismatic slip relative to basal slip when alloying elements are added to magnesium alloy [35]. Though basal slip is preferred in Mg alloys, Zr alloys favor to prismatic slip [16,36], however, both Mg and Zr crystal structure are hcp, similarly, Here we define two parameter χ 1 , χ 2 to estimate the Nb concentration and temperature effects on the trend of promoting pyramidal slip relative to prismatic slip. The χ 1 , χ 2 is defined as: Zr,0 is set to unity 1.0, represents the unstable SFE ratio of pyramidal I plane relative to prismatic I plane of pure Zr at 0 K, ( ) / g usf I P Nb,K represents the unstable SFE ratio of Zr-Nb alloys with different Nb concentration at finite temperature. Similarly, parameter χ 2 is defined to estimate the effect of pyramidal II plane. Figure 6 displays the values of χ 1 , χ 2 of Zr-Nb alloys with different Nb concentrations at 0 K and 300 K by two interatomic potentials. It can be observed from figure 6 that almost all the values of χ 1 , χ 2 calculated at 300 K are greater than those at 0 K, indicating that pyramidal slip are more inclined to activate at 300 K compared with 0 K, thus it may suggest temperature can improve plastic deformation anisotropy of Zr-Nb alloys.
Legrand et al have suggested an R parameter to judge whether basal or prismatic plane prefers to slip [37]. The criterion is R=γ basal /γ prism ×C 66 /C 44 . It is believed that Nb concentration has little effect on the C 66 /C 44 of Zr at finite temperature. We simplify the R=γ basal /γ prism . The value of R can be calculated and is shown in figure 7. One can see that R value increases gradually with the temperature increases from 300 K to 600 K, implying that the temperature increase can improve the activity of basal slip. Therefore, it could be concluded that the temperature increase can improve plastic deformation anisotropy of Zr-Nb alloys effectively.

Conclusions
In this work, MD simulations are applied to investigate the GSFE of basal, prismatic I, pyramidal I and II plane as a function of Nb concentration and temperature for Zr-Nb alloys by EAM and ADP. The following conclusions can be drawn: (1) In the composition range of Nb=0, 0.5, 1.0, 1.5, 2.0, 2.5 at.% at 0 K, the stable and unstable SFEs of different slip systems in Zr-Nb alloys show no apparent change with the increasing Nb concentration.
(2) Basal, pyramidal I and II planes of Zr-Nb alloys tend to slip by full dislocation with the increasing temperature.
(3) The plastic deformation anisotropy of Zr-Nb alloys are improved with the increasing temperature.