International Journal of Nuclear Energy
Volume 2014 (2014), Article ID 912369, 6 pages
http://dx.doi.org/10.1155/2014/912369

Research Article

Molecular Dynamics Study of Hydrogen in α-Zirconium

Ravi Kiran Siripurapu,¹ Barbara Szpunar,² and Jerzy A. Szpunar¹

¹Department of Mechanical Engineering, University of Saskatchewan, 57 Campus Drive, Saskatoon, SK, Canada S7N 5A9
²Department of Physics and Engineering Physics, University of Saskatchewan, 116 Science Place, Saskatoon, SK, Canada S7N 5E2

Received 1 June 2014; Revised 20 September 2014; Accepted 12 October 2014; Published 9 November 2014

Academic Editor: Alexander B. Shick

Copyright © 2014 Ravi Kiran Siripurapu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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⁻⁷ cm²/sec and at 1200 K has a value 1.47 10⁻⁴ cm²/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  cm²/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 ZrO₂ and ZrH₂ using ab initio calculations [17]. The cut-off radius for atomic interactions for MEAM is chosen as 4.5 Å.

Table 1: Parameters used for elements Zr and H are the cohesive energy (eV), the equilibrium nearest-neighbor distance (Å), the exponential decay factor α for the universal energy function, the scaling factor for the embedding energy, the exponential decay factors for the atomic densities, the weighting factors for the partial electron densities, and the atomic density scaling .

Element

Zr
6.325
3.21
4.46
1.03
3.7
1.0
5.7
1
1
5
4.5
−3.5
0.63

H
2.235
0.74
2.96
2.5
2.96
3.0
3.0
3.0
1
0.2
−0.1
0
16

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.

Table 2: Calculated lattice parameters (Å) and cohesive energy (eV) for zirconium, in comparison with experimental data.

Parameters
MEAM
EAM
Experimental values

3.216
3.234
3.232 [16], 3.233 [15], 3.231 [12]

5.191
5.167
5.182 [16], 5.149 [15], 5.148 [12]

1.614
1.597
1.603 [16], 1.593 [12]

Cohesive energy
−4.647
−6.634
−6.32 [16], −6.36 [12]

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)

Figure 1: Increase of lattice parameters with increase in temperature for EAM and MEAM potentials of zirconium and comparison with experimental data for (a) lattice “” and (b) lattice “.”

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 .

Table 3: Calculated stiffness (GPa) matrix using EAM and MEAM potentials for zirconium, in comparison with published literature and experimental data.

Stiffness values
EAM
MEAM
Literature [13, 14]
Experimental values [21]

141.59
135.27
147
153

74.27
88.86
69
67

74.08
63.12
74
65

167.72
170.1
168
172

43.93
18.79
44
36

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.

Table 4: Calculated Young’s modulus (GPa), bulk modulus (GPa), and Poisson’s ratio for zirconium, in comparison with experimental data.

[Direction]

EAM
MEAM
Experimental values [22]

Young’s modulus (GPa)

    2    0, 0 1 0
93
74
98

    0 0 0 1
117
135
125

Poisson’s ratio

    2    0, 0 1 0
0.33
0.37
0.32

    0 0 0 1
0.34
0.28
0.30

Bulk modulus (GPa)
96.7
98.94
96.48

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.

Figure 2: Zirconium (brown) atoms distributed inside hydrogen (blue).

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)

Figure 3: Unstable zirconium (brown) atoms in the presence of hydrogen (blue) obtained using NPT at 1100 K using timestep of 0.00001 ps. (b) Stable structure observed using NVE at temperature of 1100 K using timestep of 0.0005 ps for (a) MEAM and (b) EAM.

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.

Figure 4: Logarithmic diffusion coefficients calculated using EAM and MEAM potentials for hydrogen in bulk zirconium at a temperature of 500–1200 K and compared with experimental data.

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.

Table 5: Calculated activation energies (kcal/mol) of hydrogen in zirconium using EAM and MEAM potentials, in comparison with experimental data.

Method
Activation energy (kcal/mol)

Makowitz [9]
Experiments
1.6 to 6.5

Johnston et al. [10]
Experiments
6.4

Kearns [7]
Experiments
10.83

Sawatzky [8]
Experiments
3.8 to 6.1

EAM
Present
11.302

Hydrogen in Nickel [16]
Ab-initio-(EAM)
11.618

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.

Figure 5: Calculated anisotropic hydrogen displacements are seen to be higher along the basal plane for zirconium.

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.

References

M. P. Puls, “The influence of hydride size and matrix strength on fracture initiation at hydrides in zirconium alloys,” Metallurgical Transactions A, vol. 19, no. 6, pp. 1507–1522, 1988.

D. G. Hurst, Canada Enters the Nuclear Age: A Technical History of Atomic Energy of Canada Limited, McGill-Queen's University Press, Quebec, Canada, 1997.

C. Korn and S. D. Goren, “NMR study of hydrogen diffusion in zirconium hydride,” Physical Review B, vol. 33, no. 1, pp. 68–78, 1986.

J. A. Szpunar, W. Qin, H. Li, and N. A. P. K. Kumar, “Roles of texture in controlling oxidation, hydrogen ingress and hydride formation in Zr alloys,” Journal of Nuclear Materials, vol. 427, no. 1–3, pp. 343–349, 2012.

C. J. Simpson and C. E. Ells, “Delayed hydrogen embrittlement in Zr-2.5 wt% Nb,” Journal of Nuclear Materials, vol. 52, no. 2, pp. 289–295, 1974.

N. A. P. Kiran Kumar, J. A. Szpunar, and H. Zhang, “Preferential precipitation of hydrides in textured zircaloy-4 sheets,” Journal of Nuclear Materials, vol. 59, pp. 7010–7021, 2010.

J. J. Kearns, “Diffusion coefficient of hydrogen in alpha zirconium, Zircaloy-2 and Zircaloy-4,” Journal of Nuclear Materials, vol. 43, no. 3, pp. 330–338, 1972.

A. Sawatzky, “Hydrogen in zircaloy-2: its distribution and heat of transport,” Journal of Nuclear Materials, vol. 2, no. 4, pp. 321–328, 1960.

J. M. Makowitz, “Westinghouse Atomic power Report,” 1959

W. V. Johnston, W. R. Jacoby, J. S. Wallam, and A. H. Alberts, “Knolls atomic laboratory report,” Tech. Rep., 1959.

N. Singh, D. K. Avasthi, and A. Tripathi, “Electron structure and activation energy of hydrogen in α-Zr using nonlinear response theory,” Bulletin of Materials Science, vol. 20, no. 3, pp. 349–358, 1997.

S. M. Foiles, M. I. Baskes, and M. S. Daw, “Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys,” Physical Review B, vol. 33, no. 12, pp. 7983–7991, 1986.

M. I. Mendelev and G. J. Ackland, “Development of an interatomic potential for the simulation of phase transformations in zirconium,” Philosophical Magazine Letters, vol. 87, no. 5, pp. 349–359, 2007.

M. Szpunar, P. I Baskes, and T. Jochym, “Oxygen and hydrogen diffusion in Zirconia,” in Proceedings of the 8th International Conference on CANDU Fuel, pp. 551–557, Honey Harbour, Canada, 2003.

M. I. Baskes, “Modified embedded-atom potentials for cubic materials and impurities,” Physical Review B, vol. 46, no. 5, pp. 2727–2742, 1992.

B. Szpunar, L. J. Lewis, I. Swainson, and U. Erb, “Thermal expansion and hydrogen diffusion in nanocrystalline nickel,” Physical Review B: Condensed Matter and Materials Physics, vol. 60, no. 14, pp. 10107–10113, 1999.

B. Szpunar, M. I. Baskes, and P. T. Jochym, “Oxygen and Hydrogen diffusion in Zirconia,” in Proceedings of the 8th International Conference on CANDU Fuel, pp. 551–557, Honey Harbour, Canada, 2003.

B. Szpunar, M. I. Baskes, J. A. Szpunar, and P. T. Jochym, “Atomistic simulations Of Zr–O–H systems,” in Proceedings of the 9th International CNS Conference on CANDU Fuel, Belleville, Ontario, Canada, September 2005.

M. I. Baskes and R. A. Johnson, “Modified embedded atom potentials for HCP metals,” Modelling and Simulation in Materials Science and Engineering, vol. 2, no. 1, pp. 147–163, 1994.

J. Goldak, L. T. Lloyd, and C. S. Barrett, “Lattice parameters, thermal expansions, and grüneisen coefficients of zirconium, 4.2 to 1130°K,” Physical Review, vol. 144, no. 2, pp. 478–484, 1966.

D. Tromans, “Elastic anisotropy of HCP metal crystals and polycrystals,” International Journal of Research Reviews in Applied Sciences, vol. 6, no. 4, pp. 462–483, 2011.

Y.-M. Kim, B.-J. Lee, and M. I. Baskes, “Modified embedded-atom method interatomic potentials for Ti and Zr,” Physical Review B—Condensed Matter and Materials Physics, vol. 74, no. 1, Article ID 014101, 2006.