Abstract
Molecular dynamics approach is used to simulate hydrogen (H) diffusion in zirconium. Zirconium alloys are used in fuel channels of many nuclear reactors. Previously developed embedded atom method (EAM) and modified embedded atom method (MEAM) are tested and a good agreement with experimental data for lattice parameters, cohesive energy, and mechanical properties is obtained. Both EAM and MEAM are used to calculate hydrogen diffusion in zirconium. At higher temperatures and in the presence of hydrogen, MEAM calculation predicts an unstable zirconium structure and low diffusion coefficients. Mean square displacement (MSD) of hydrogen in bulk zirconium is calculated at a temperature range of 500–1200 K with diffusion coefficient at 500 K equals 1.92 10−7 cm2/sec and at 1200 K has a value 1.47 10−4 cm2/sec. Activation energy of hydrogen diffusion calculated using Arrhenius plot was found to be 11.3 kcal/mol which is in agreement with published experimental results. Hydrogen diffusion is the highest along basal planes of hexagonal close packed zirconium.
1. Introduction
The behavior of zirconium and its alloys under various operating conditions in nuclear reactors has been extensively studied. Zirconium (Zr) is a transition metal with strong anisotropic physical properties due to a hexagonal close packed (HCP) crystal structure. It is used in nuclear reactors because of low thermal neutron absorption and good corrosion resistance at high temperatures [1]. CANada Deuterium Uranium (CANDU) and almost all other nuclear reactors have structural elements made of zirconium alloys [2] making zirconium alloys the perfect choice for fuel cladding and pressure tubes in nuclear fuel channels.
The hydrogen generated during reactor operation has deleterious influence on mechanical properties of zirconium. The hexagonal close packed structure of zirconium consists of both tetrahedral and octahedral sites, and hydrogen is expected to diffuse and occupy these sites [3]. Zirconium can dissolve up to 450 wtppm of hydrogen at a temperature of 500°C, but the solubility drastically decreases to 65 wtppm at 300°C and it further decreases to 0.05 wtppm at room temperature. The low solubility of the hydrogen with temperature leads to formation of brittle hydride platelets [4]. These hydrides are brittle and crack on application of stress, which reduces life expectancy of nuclear components thereby increasing the cost of nuclear power.
Hydride precipitates play an important role in the hydrogen embrittlement of various zirconium alloys [1–10]. Fracture toughness tests under tension was performed on specimens of bulk δ-hydride obtained extremely low fracture toughness of ~1 MPa m1/2 at room temperature [5]. Because hydride in Zr alloy specimens may contain defects, such as voids or microcracks, and there is often misfit stress between the hydride and zirconium matrix, the results obtained for pure hydride may not be directly applied to the hydride in the reactor.
These hydrides are seen to preferentially precipitate in [6] crystallographic orientation. Under the influence of stress the hydrogen in hydride starts diffusing and is seen to reprecipitate in a direction perpendicular to applied hoop stress [4]. However the cause of diffusion of hydrogen into zirconium itself is not very well understood and is therefore of importance.
Hydrogen diffusion in zirconium was studied using diffusion sample method by Kearns [7]. This method involves measuring the hydrogen concentration along the length of samples to hydrogen free studs at a temperature range of 548–973 K. The diffusion coefficient of hydrogen in zirconium was found to be cm2/sec with activation energy of 10.7 kcal/mole. Activation energy obtained by Sawatzky [8] using temperature gradient technique for hydrogen diffusion was in the range of 3.8–6.1 kcal/mole and values ranging between 1.6–6.0 kcal/mole was obtained by Markowitz [9] using similar approach. Further Johnston et al. [10] have reported the value of 6.4 kcal/mole, and the scatter of these results clearly indicates a wide range of inconsistency in determining the activation energy of hydrogen in zirconium. Theoretical calculations [11] using spherical solid model potential (SSMP) in α-Zr predicted the activation energy of 0.285 eV or 6.57 kcal/mole. However the accuracy of these activation energies is still questionable due to material compositions and metallurgical structures.
In order to understand problem of hydrogen diffusion in zirconium we need to analyze different empirical methods for zirconium-zirconium interactions and zirconium-hydrogen interactions and put forward the best potential of interaction. The goal of this work is to understand hydrogen diffusion in bulk zirconium using molecular dynamics.
2. Simulation Methods
Calculations were carried out using molecular dynamics (MD) approach with the interactions between atoms represented by the embedded atom method (EAM) [12–14] and modified embedded atom method (MEAM) [15]. EAM hydrogen was previously used in calculation of diffusion coefficients in bulk and grain boundaries of Nickel [16]. In this research we introduce EAM potential to calculate diffusion coefficients and activation energies of hydrogen in bulk zirconium.
Modified embedded atom method (MEAM) functions have been developed [15], and this method describes very well the energetic, structural, and mechanical properties of solid metals and gases. The MEAM is an extension of embedded atom method with the inclusion of directional bonding for interaction of gases with metal atoms. MEAM functions for zirconia and was developed [17, 18] as shown in Table 1 and was successfully tested for ZrO2 and ZrH2 using ab initio calculations [17]. The cut-off radius for atomic interactions for MEAM is chosen as 4.5 Å.
However the MEAM potential needs to be tested for zirconium and zirconium-hydrogen system for studying diffusion of hydrogen in zirconium. Results of both EAM and MEAM are compared for zirconium as well as hydrogen diffusion in zirconium and will be used to propose the best potential for the present calculation.
3. Test of Methods Used for Calculation
3.1. Lattice Parameters
Both EAM and MEAM were used for the prediction of zirconium structure. It is well known that zirconium has HCP structure with different “” and “” lattice parameters. Lattice parameters are calculated using conjugate gradient (CG) algorithm at 0 K. An external pressure or stress tensor is applied to the electronic structure of zirconium and the energy is then minimized. At the lowest energy of structure with respect to the positions of atoms, the lattice parameters and cohesive energy which follows “universal of binding curve” were recorded [12].
The results presented in Table 2 indicate good agreement of lattice parameters of EAM with experimental data, and MEAM is seen to have a good agreement of with ideal ratio of zirconium which is 1.633 [19]. This result indicates that there is a better agreement of experimental lattice parameters with EAM potential than MEAM potential. However, MEAM results agree better than EAM with the experimental value of cohesive energy.
3.2. Thermal Expansion
We know that zirconium and its alloys in the nuclear reactor core are exposed to very high temperatures; it is therefore important to test the methods at various temperatures and evaluate changes of structure and properties.
Minimized structure of zirconium with lattice parameter () is equilibrated to the constant temperature using NPT with a timestep of 0.001 ps. The average value of lattice parameter () is calculated in a temperature range 300–1125 K.
As presented in Figure 1 the increase of temperature shows the increase of lattice parameters of HCP crystal of zirconium predicted by both EAM and MEAM. The thermal expansion coefficient for EAM (1/K) and (1/K) and for MEAM (1/K) and (1/K). Experimental calculations by Goldak et al [20] indicate (1/K) and (1/K). This means the thermal expansion is higher along “” direction when compared to “” direction of zirconium and this is in a good agreement with the experimental data. The potentials show good stability even at higher temperatures. EAM has a better agreement with experimental data than MEAM.
(a)
(b)
3.3. Mechanical Properties
Young’s modulus and Poisson’s ratio are the two most important properties considered in the design of engineering materials. The potentials must be verified by calculating these two properties in order to use it for other mechanical properties calculation. In order to obtain Young’s modulus and Poisson’s ratio the stiffness matrix has to be calculated.
The structure is subjected to positive (tensile) and negative (compressive) displacements and the resultant of the two is used for one row of the stiffness matrix [21]. This process is repeated for 3D structure in order to obtain stiffness constants () using Voigt [21] average in different directions as zirconium structure has anisotropic elastic properties as shown in Table 3. From the obtained stiffness matrix the compliance matrix () is calculated as the inverse of the matrix. Thus Young’s modulus and Poisson’s ratio in two directions can be obtained from (1) and (2). The bulk modulus is calculated using (3). Consider where .
As indicated in Table 3 MEAM shows lower stiffness along the shear directions () of zirconium indicating weaker interactions along that direction. Table 4 shows the calculated mechanical properties. Although both potentials show anisotropic Young’s modulus along the two calculated directions, lower Young’s modulus is observed along the basal plane direction as compared to the experimental data. Further, the calculated bulk modulus for both potentials is in good agreement with experimental data.
From the above tests for lattice constants, thermal expansion, and mechanical properties, we can conclude that EAM has a better agreement with literature than MEAM. Further, investigations are carried out to study diffusion of hydrogen in zirconium.
4. Results
4.1. Diffusion of Hydrogen in Zirconium
Hydrogen diffusion inside bulk zirconium is studied by modeling hexagonal crystal structure consisting of 6400 atoms of zirconium and 100 atoms of hydrogen as shown in Figure 2. Diffusion coefficients of hydrogen in zirconium are calculated using EAM and MEAM.
4.2. Mean Square Displacement (MSD)
The time-dependent diffusion coefficient was derived from mean square displacement. For the calculation of the diffusion coefficients of hydrogen, zirconium structure is minimized using energy minimization. NPT (constant temperature) equilibration is done until the constant temperature and pressure are reached. NVE (constant volume) conditions are then applied with timestep of 0.0001 ps for runs up to 2000 ps. Increase in pressure that follows the changes in the structure of zirconium in presence of hydrogen is observed for longer runs using MEAM potential as shown in Figure 3(a). Calculation using EAM potential zirconium structure is seen to be stable in the presence of hydrogen at higher temperatures of 1200 K as shown in Figure 3(b).
(a)
(b)
4.3. Diffusion Coefficients and Activation Energies
The diffusion coefficients () are calculated using Einstein relation as in the following equation: where ( is the average over all possible initial times) is the mean square displacement and is the dimensionality of space ( for bulk diffusion). Calculated diffusion coefficients of hydrogen increase with the increase in temperature as shown in Figure 4.
EAM potential of zirconium and hydrogen is seen to be in good agreement with experimental data. MEAM records some change of structure of zirconium and therefore lower hydrogen diffusion coefficients are obtained.
Activation energies () are calculated using diffusion coefficients for a temperature range 500–1200 K using the following equation: where “” is universal gas constant 8.31 J/mol-K and and are diffusion coefficients for temperatures and , respectively.
On comparison with various experimental results, as shown in Table 5, it has been seen that diffusion coefficients of hydrogen in bulk zirconium using MEAM potential are found to be in good agreement with experimental data. Table 5 shows calculated activation energy of hydrogen in zirconium using EAM that is in a good agreement with Kearns [7]. Similar activation energy and diffusion coefficients are observed for hydrogen diffusion in Nickel which has a face centered cubic structure.
4.4. Anisotropic Hydrogen Displacements
As noted above EAM potential for hydrogen in bulk zirconium has shown a good agreement of calculated diffusion coefficients and activation energies with the observed experimental data. The anisotropy of hydrogen movements during the diffusion is tested. Diffusion coefficients of hydrogen are calculated in 3 different directions of zirconium crystal , , and . MSD for hydrogen movements is calculated for a temperature of 500–1200 K and (4) is used in the calculation of the diffusion coefficients. As shown in Figure 5 hydrogen movements are higher along the two basal directions and which is more likely. These indicate preferential hydrides formation along these directions which are experimentally observed by Szpunar et al. [4] leading to preferential orientation of hydride along these directions.
5. Discussion
EAM used for zirconium structure and property calculation has shown a better agreement of lattice parameters and mechanical properties with available experimental data than MEAM method developed by Baskes [16]. By using EAM potential, stable zirconium structure is obtained when hydrogen is placed in the interstitial sites of zirconium. MEAM potential is seen to make zirconium structure unstable when hydrogen is present resulting in lower diffusion coefficients. MEAM potential could be possibly used for predicting hydride structure. Diffusion coefficients and activation energies obtained through EAM for hexagonal close packed structure of zirconium are similar to face center cubic structure of Nickel [19].
6. Conclusion
(1)Good agreement of lattice parameters calculated using EAM with experimental literature data is obtained.(2)MSD for hydrogen is tested for timesteps ranging from 0.00001 ps to 0.0001 ps using both EAM and MEAM. Structure of zirconium in the presence of hydrogen at higher temperatures as calculated with MEAM potential was unstable.(3)Diffusion and activation energies for hydrogen diffusion obtained through EAM potential are in good agreement with the experimental data.(4)Anisotropy of hydrogen diffusion is predicted; diffusion is higher along the basal directions of zirconium crystal and this observation may be used to justify formation of hydrides along these directions.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the Natural Sciences and Engineering Research Council of Canada (NSERC) for the financial support. Also, they are grateful to WestGrid (Compute/Calcul Canada) for computing resources. The authors would also like to thank S. M. Folies, M. I. Baskes for providing potentials.