Anales AFA Vol. 32 Nro. 4 (Enero 2022 – Abril 2022) 120-127
https://doi.org/10.31527/analesafa.2021.32.4.120
Materia Condensada
DIFUSIÓN DE HIDRÓGENO EN HCP Zr
HYDROGEN DIFFUSION IN HCP Zr
C. García1 y V. P. Ramunni*2,3
1 Instituto Sabato - UNSAM/CNEA, Av. Gral. Paz 1499, (1650) San Martín - Argentina.
2 CONICET Godoy Cruz 2390 (C1425FQD) - Argentina.
3 Gcia. Materiales-CNEA, Av. Gral. Paz 1499, (1650) San Martín - Argentina.
Autor para correspondencia: email: vpram@cnea.gov.ar
ISSN 1850-1168 (online)
Recibido: 02/06/2021 Aceptado: 27/07/2021
Resumen:
Estudiamos el comportamiento del hidrógeno (H) en las fases puras α, β Zr y β Nb, empleando los códigos de
primeros principios Siesta y Vasp, basados en la teoría del funcional de la densidad (DFT). Calculamos la energía
de solución del hidrógeno, la energía de unión de los complejos H-H y H-Nb, como también la energía de mezcla
del sistema Zr-Nb. Además, propusimos un modelo simple para el estudio de la difusión del H (DH) mediado por
un mecanismo intersticial en la fase hcp del Zr. Encontramos que el hidrógeno intersticial difunde isotrópicamente
según DH = 1.1×10−7 exp(−42KJ.mol−1/RT) (m2/s). Nuestros resultados están en buen acuerdo con los datos
experimentales y otras técnicas numéricas de Montecarlo cinético (KMC). Calculamos también, las energías de
formación de los tres hidruros más estables en la fase hcp Zr.
Palabras clave: Hidrógeno, Difusión, Hidruros de circonio, α Zr, cálculos ab-initio.
Abstract:
We study the behavior of hydrogen (H) in the α and β phases of Zr and in the β phase of Nb, using the Siesta and
Vasp first principles codes based on the density functional theory (DFT). We calculate the hydrogen solution
energy, the binding energy of the H-H and H-Nb complexes, as well as, the mixing energy of the Zr-Nb system.
Furthermore, we have proposed a simple model for the study of the hydrogen diffusion (DH) mediated by
interstitial mechanism in hcp Zr. We find that interstitial hydrogen diffuses isotropically according to DH =
1.1×10−7 exp(−42KJ.mol−1/RT) (m2/s). Our results are in good agreement with experimental data and other Kinetic
Monte Carlo (KMC) calculations. Finally, the formation energies of the three most stable hydrides in the hcp Zr
phase were also calculated.
Keywords: Hydrogen, Diffusion, Zirconium-Hydrides, Zr-Nb, ab-initio calculations.
I. INTRODUCCIÓN
In the context of the Hydride Delayed Cracking phenomena (DHC), Zirconium based alloys are currently used in
reactors due to their low neutron absorption and their good response to high temperature corrosion. Since hydrogen
can be present in the alloys as a remnant impurity of the manufacturing process, and can also enter during the
service life of the material, it is very important to evaluate the response of the alloys to hydrogen.
The microstructure, solubility and diffusion coefficient of H in the alloy are variables that allow us understand the
DHC phenomena. For example, the Zr-2.5Nb alloy is bi-phasic with grains of phase
(hcp, with 0.6 % Nb)
surrounded by plates of
Zr
phase (bcc, with 20 % Nb).
With thermal treatments, the original plates of
Zr lose their continuity and interrupts the fast paths for the
hydrogen diffusion in the axial direction of the tube and, at the same time, their percentage of Nb is altered
evolving to the stable phase
Nb (bcc, with 90 - 100%Nb).
From this experimental evidence, the following questions arise:
1. Since
Nb must have a hydrogen solubility lower than that of pure Nb, then: It is possible to calculate the
solubility of H in pure Nb and of the Zr-Nb alloys, with different percentages of Nb?
2. Is the loss of continuity of the original sheets effectively responsible for the observed decrease in the diffusion
coefficient of H as time passes at a given temperature?
3. Is the diffusion and / or solubility of hydrogen in the alloy affected by the Nb content?
This is a work prior to that carried out in Ref.,1 in which the authors have successfully answered questions 2 and 3
posed by our experiments. Here we focused on the study of the hydrogen’s diffusivity on
Zr.
Concerning atomic diffusion, several efficient methods have been developed in recent years for finding activated
states or, mathematically, saddle points from which diffusion can be studied. Here we employ the Monomer one,
developed previously in our laboratory2 and The Dimer method3 as implemented in Vasp.4 The Monomer, in
coupled to Siesta’s code5 to be interfaced as a subroutine.6
The Monomer computes the least local curvature of the potential energy surface using only the forces furnished by
Siesta. The force component along the corresponding eigenvector is then reversed (pointing ‘up hill’), thus defining
a pseudo force that drives the system towards saddles. Both, local curvature and configuration displacement stages
are performed within independent conjugate gradients loops. The method is akin to the Dimer one from the
literature,3 but roughly employs half the number of force evaluations.
These methods are very efficient at obtaining the migration barriers from which we study the atomic diffusion in
metals. In this way, Ishioka and Koiwa7 have studied the interstitial diffusion in hcp metals, involving octahedral
and tetrahedral sites. They have proposed an on-lattice random walk model to derive the diffusivities of impurity
atoms on a crystal lattice containing multiple sublattices, such as the octahedral and tetrahedral sites in hcp
crystals.
Here we focused on the study of the hydrogen’s behavior in Zr-Nb. In this sense we study some relevant diffusion
parameters, namely: (i) the hydrogen solubility, (ii) the hydrogen binding energies of the H-H and H-Nb
complexes, (iii) the mixing solution energy of Nb-Zr and (iv) the hydrogen diffusion coefficients in hcp Zr. With
this purpose, we perform ab-initio calculations with Siesta5 and Vasp,4 in order to obtain the hydrogen migration
barriers in the bulk of
Zr using the model developed in Ref..7 Classical calculations were also performed in
order to check some structure consistence.
The present paper is structured as follows: Section II describes our calculation methodology with the DFT based
codes as Siesta and Vasp, and from classical codes. Sections III show the way as the lattice parameter is calculated
from the fit of a third-order Birch-Murnagahn equation of state, EOS,.8 Section IV, is devoted to summarize the
energetic of Hydrogen in the Zr-Nb system. Section V describes the way as the H-diffusion coefficients were
calculated in
Zr together with our numerical results. In brief, Section VI, describes the structure of the three
stable zirconium-hydrides in
Zr. The last section presents some conclusions.
II. CALCULATION METHODOLOGY
Ab initio calculations were performed as implemented in Siesta and Vasp. Siesta, is a freely available code,5 which
is based on pseudopotentials and numerical, atomic-like, basis functions. The pseudopotential and basis for Zr are
according to Ref..6 Pseudopotential for H was downloaded from Siesta’s home page, and the corresponding base
was automatically generated using the code itself.
Present Siesta calculations were carried out within the generalized gradient approximation (GGA) for exchange
and correlation, according to the PBE parameterization described in Ref..9 A spatial mesh cutoff of 460.0 Ry were
used, with a smearing temperature of 0.15 eV, within a Fermi–Dirac scheme. Reciprocal space is partitioned in a
555
Monkhorst–Pack grid. All calculations were atomically relaxed, though the cell boundary remains fixed.
For hydrogen, we used a so called TZDP basis with a split norm parameter of 0.5, which is the accurate one that
most efficiently relaxes the bulk containing H.
Vasp uses a projector augmented-wave formalism to describe the interactions between atoms.10 An energy cutoff
of 460 eV and fine
-centered k-point meshes automatically generated with the Monkhorst-Pack scheme11 of
555
-k points mesh ensured electronic convergence of the calculations. For the term of exchange and
correlation we use the PBE-GGA approximation, with the Methfessel-Paxton broadening scheme with a
0.3
eV
width. Concerning optimum structures, relaxation was done until the residual forces were below 0.02 eV/Å.
In both, Siesta and Vasp, we have used two supercell sizes of 48 and 54 atoms, respectuvely for
Zr and
Nb.
The supercell of
Zr with 48 atoms is constructed by reproducing the conventional cell, of 4 atoms, along the
lattice vectors of the hexagonal primitive cell a, b and c, in
232
.
In addition, we have performed calculations with the molecular statics technique (MS), implemented in an in-house
code, provided with suitable EAM interatomic potentials for Zr.12 We use a simulation cell of
15 10 10
and
20 20 20
respectively containing 6000 Zr or Nb atoms, with periodic boundary conditions.
These methods are very efficient at obtaining the equilibrium positions of the atoms by relaxing the structure via
the conjugate gradients technique. Impurity and defect relaxation, includes interstitial H atoms in
Zr,
Nb, as
well as, in the Zr-Nb alloy structure. The needed migration barriers of several jumps were calculated coupling both
EAM and Siesta techniques to the Monomer method,2,6 while we used the DIMER as implemented in Vasp.3
Finally, our calculations were carried out at constant volume, and therefore the enthalpic barrier
H U p V
is equal to the internal energy barrier
U
.
The empirical formula (EF)
The EF is an expression that represents the simplest proportion in which the atoms are present in the alloy. In the
case of the Zr-2.5Nb alloy, we start from a hypothetical sample with 100 gr of Zr-2.5Nb. From the % of Nb we
defined that in 100 gr of the sample there are 2.5 gr of Nb while the rest, (100-2.5) = 97.5 gr, are of Zr. The moles
of Zr and Nb that are present in those 100 gr of Zr-2.5Nb, are obtained through the ratio between the grams of each
species in the sample by their respective molecular weights,
91.224
Zr
M
P
and
92.90637
Nb
M
P
, as,
97.5 /91.224 1.0057 ,gr Zr molesof Zr
(1)
2.5 /92.90637 0.0269 .gr Nb molesof Nb
(2)
In Eqs. (1) and (2) the ratio between species that contains the least number of moles gives the empirical formula
37 1
:Zr Nb
for the alloy. This assumes that Nb atoms are located at periodic positions and not randomly; Also, this
does not consider the different compositions of the alpha and beta phases.
From this result, we were considered two sizes of supercells, one of 36 atoms of Zr and another one with 48 atoms
of Zr.
III. STRUCTURE OPTIMIZATION
For
Nb and
Zr, in order to better describe bulk zirconium properties, we determine a cut-off for the needed
size of k-point mesh. On a
111
cell consisting of two Nb or Zr atoms, several calculations were done using k-
point meshes with sizes of
12 12 8
,
16 16 10
,
20 20 14
,
24 24 16
,
28 28 18
and
32 32 20
with the
use of 95, 180, 352, 549, 800 and 1122 irreducible k-points, respectively.
Using a general convergence cut-off of 5 meV/atom, it was shown that the choice of
20 20 14
k-point mesh
would be sufficient to represent the both, Nb and Zr
111
supercells while keeping the computational effort at
the minimum level. In subsequent calculations, choice of k- point mesh would be made based on the ratio between
the size of the new supercell and that of the
111
supercell.
With the chosen k-point mesh, equation of state calculations were done in determining the computational lattice
parameters of bulk
Zr and
Nb as follows: We perform an automated volume scan by means of a shell script
loop.sh, with lattice parameters above and below the experimental value. Then, the procedure was to get the energy
as a function of volume using by regression fitting DFT results to third-order Birch-Murnagahn equation of state,
EOS.8
3
2/3 '
0 0 0
00
2
2/3 2/3
00
9
( ) ( ) 1
16
( ) 1 6 4( )
V B V
E V E B
V
VV
VV
(3)
In Eq. (3),
0
V
was found from the regression fitting in Fig. 1, which corresponds to the cell volume that result in
the lowest possible total energy of the system. One last calculation to calculate the Zr lattice parameter was to
optimize the ionic structure of a zirconium cell with a fixed volume equal to the value of
0
V
. In addition, by getting
the energy as a function of volume using third-order Birch-Murnagahn EOS, information about bulk modulus was
also found.
The optimization procedure is also used with hcp structures. The fitting of the best
a
and
c
or
/ca
for hcp
crystals is more difficult than in the previous case of the cubic lattices. This is because for hcp alloys, one should
define a global minimum of the total energy of a crystal as a function of both, the volume of the unit cell and the
structural relation
/ca
.
We have obtained a lattice parameter of
3.304
Nb
a
Å and
3.301
Nb
a
Å for a supercell containing 16 and 54 Nb
atoms, respectively.
Also, the above procedure, is performed in order to obtain the optimal lattice parameter of the structure including
interstitial (H) and substitutional (Nb) defects. For
Zr and
Zr containing impurities, a test value of structural
relation
/ca
was fixed. The standard optimization procedure, like in the case of cubic lattices, was used to define
the optimal volume of the unit cell corresponding to the minimum of the total energy as a function of the unit cell
volume,
()( / )
i
min i
E c a
, where i=1,...,N, the number of test values for the structural relation
/ca
.
The obtained empiric data as the dependence
( / )
min
E c a
, were approximated by the functional dependence. The
minimum of this dependence corresponds to the optimal value of the structural relation
( / )opt
ca
. Next, for the
determined
( / )opt
ca
, the standard optimization procedure was made to determine a global minimum of the total
energy and the optimal unit cell volume, respectively.
Using the standard definition for the volume of the HCP unit cell
32
( / ) ( )
3
hcp opt hcp
V c a a sin
the optimal lattice
constant
hcp
a
was determined for each case. Table 1, summarizes present results of the optimized lattice parameters
from Zr and Nb.
Table 1 reports the calculated lattice parameters of
,
Zr and
Nb, namely: The lattice parameters,
a
and
/ca
, the Bulk modulus,
B
, the Cohsive energy,
c
E
, the vacancy formation energy,
V
f
E
, and the vacancy
migration energy,
V
m
E
, in comparison with experimental and theoretical values reported in the literature.
Fig. 1, show our results of the optimized lattice parameter from the fit of Eq. (3) of a supercell containing 16
atoms of
Nb. In this case the Bulk modulus was obtained, as well as, the resulting pressure of the supercell after
relaxation. The same procedure was followed from 2 up to 54 atoms of Nb and Zr atoms with bcc structure.
FIG. 1: Optimized lattice parameter of a
β
Nb supercell containing 16 atoms.
TABLE 1: Lattice parameter (in Å) and the vacancy formation and migration energy (in eV).
Zr
Vasp/Siesta
Vasp
Vasp
EAM
Exp.
N=48
N=36
Ref.13
N=6000
a
3.237/3.239
3.240
3.237
3.2171
3.23214
c
5.176/5.185
5.186
5.183
-
5.14115
/ca
1.599/1.601
1.601
1.601
1.6046
1.592914
f
v
E
1.89/1.96
2.202
1.79
1.75015
c
E
6.25
6.28
6.34
6.25
6.2516
v
f
E
1.79
2.00
2.02
1.7416
v
m
E
in
c
(eV)
0.61
0.61
0.61
0.6017
v
m
E
in
a
(eV)
0.66
0.65
0.66
0.6218
Nb
Vasp/Siesta
Vasp
Vasp
EAM
Exp.
N=108
N=54
Ref.19
Ref.20
a
(Å)
3.303/3.308
3.304
3.309
3.3000
3.30321
v
f
E
2.680/2.704
2.72
2.700
c
E
7.610/7.730
7.64
-
7.57
7.5716
v
f
E
2.76
2.81
2.8822
2.75
2.7523
v
m
E
0.65
0.70
0.5524
0.67
The comparison in Table 1 shows that the predicted geometries of the GGA (PBE) slightly overerestimated the
experimental lattice parameters. Although the regression fitting of the EOS in Eq. (3) overestimates
0
V
by about
1.0%, this discrepancy between the DFT results and the reported experimental data was acceptable for this study.
Then, the current set up of zirconium supercells with the chosen k-points mesh is approximately accurate in
representing the actual zirconium and niobium metal structure.
The dependence of the lattice parameter,
Zr
a
, for
Zr-Nb alloys versus niobium concentration was previously
evaluated in Ref..25 The authors have shown that with an increase in the niobium concentration the lattice constant
decreases according to the linear law. This result is consistent with the well known approximate empirical rule,
called in metallurgy Vegard’s law. Also, at 0% of Nb the lattice parameter of
Zr is
3.5634
bcc
Zr
a
Å while for
pure
Nb is
3.3004
bcc
Nb
a
Å. Present calculations give
3.5629
bcc
Zr
a
Å.
IV. THE ENERGETIC OF THE Zr-Nb SYSTEM
IV.1 Mixing solution energy
The mixing solution energy,
,
mix
, of a substitutional B atom in a matrix of A atoms with
or
structure, is
obtained as,19
00
[ ] [ ] ( ) [ ] [ ],
AB
N M M N M M M
E A B E A B N M E A E B
N
(4)
with
[]
N M M
E A B
the energy of the supercell with (N-M) A atoms and M B-atoms,
0[]
A
EA
the reference energy of
the matrix calculated with the same supercell and
0[]
B
EB
the energy per atom of B in its state of reference (bcc for
Nb
and hcp for
Zr
). This excess energy is the solution energy weighted by the solute concentration
BM
N
as in Ref.,19
,
[ ] ,
N M M B mix
E A B
(5)
High and positive values of Eq. (5), show a strong tendency to Zr and Nb atoms not to mix in the Zr-Nb alloy.
Here, for the diluted limit, we have used a supercell of 108 and 144 atoms respectively for bcc Nb and hcp Zr.
As
BM
N
and defining
1
Zr Nb
x
we obtain
( 0) 0.322
Nb x
eV for a substitutional Zr diluited in
Nb
, and
( 0) 0.632
Zr x
eV for a substitutional Nb atom in
Zr
. The last results are in agreement with
previous DFT ones of 0.61 eV26 and 0.68 eV.27
Our results show that Nb has a low solubility in
Zr
given the high value of 0.63 eV, while Zr has high solubility
in
Nb
. For the supercells with 54 and 48 atoms of
Nb and
Zr,
( 0) 0.34
Nb x
eV and
( 0) 0.641
Zr x
eV, respectively.
IV.2 Hydrogen solution energies in
Zr
and
Nb
The hydrogen solution energy,
S
E
, at stable sites were calculated according to the following expression,28
2
1
( , ) ( ) ( ) .
2
H
S ZPE
E E N H E N E H E
(6)
In Eq. (6)
( , )E N H
is the energy of the supercell containing N atoms plus an absorbed hydrogen at a determined
interstitial site,
()EN
the energy of supercell without defects,
2
1()
2EH
the energy of the hydrogen in the gaseous
phase, and
ZPE
E
, the zero point energy (ZPE) that corresponds to quantum corrections. The value of the solution
energy in a given site gives us information about the stability of such a site.
In Eq. (6), the third term corresponds to half of the total energy of the hydrogen molecule, is calculated using a cell
of sides
10 10 10
Å
3
containing only the dimer,
2
H
, in the center of the cell for a single k-point. While the
last term, correspond to the zero point energy due to quantum corrections.
Fig. 2 shows the obtained equilibrium bond length of 0.750 Å for the hydrogen molecule, and a dimer energy of
6.765 eV, from which we calculate a binding energy of
4.54
eV for the hydrogen molecule, in agreement with
experimental values of 0.741 Å, and 4.75 eV,29 respectively.
FIG. 2: The fitted hydrogen molecule,
2
H
, bond length.
The calculated
2
H
vibration frequency are 4401 cm
1
, that agrees with the experimental value of 4395 cm
1
.30
Finally, present result of the ZPE 0.27 eV, reproduces the experimental value 0.274 eV.29 Table 2, summarizes our
results of the H solution energy in the
Zr
and
Nb
respectively for a supercell of 48 and 54 atoms.
TABLE 2: Hydrogen solution energy,
H
S
E
(in eV) calculated with Siesta/Vasp.
System
H
S
E
H-Site
Refs.
48
Zr H
-0.481
T
-0.60930
48
Zr H
-0.407
O
-0.54930
54
Nb H
-0.323
T
-0.34031
54
Nb H
-0.052
O
-0.03031
48 2
Zr H
-0.963
T
-
48 2
Zr H
-0.637
O
-
High and negative values of the solution energies mean stable sites, while positive values can infer that the site will
be unoccupied or, the site is not favorable for hydrogen solubilization.
The results in Table 2 reveal that H is highly soluble on both the T and O sites of
Zr phase, with the following
energy relationship
( ) ( )
HH
SS
E T E O
. Concerning
Nb, H preferentially dissolves at the T site, whose solubility
is much higher than at the O site. In regarding the energy of solution for the 2H in the
Zr phase, we observe that
the T site is the preferential site for 2H whose energy of solution is more than double that of the isolated H in the
same phase, showing that there is a slight interaction between the 2H.
IV.3 Binding energies
The hydrogen binding energy in a crystal containing N atoms of Zr plus
k
atoms of Nb on substitutional sites, was
calculated with the following expression,
1 1 1
{ ( ) ( 1) ( )} { ( ) ( )}.
b N k k n N N N
E E Zr Nb H n k E Zr k E Zr Nb n E Zr H
(7)
In Eq. (7)
()
N k k n
E Zr Nb H
and
11
()
N
E Zr Nb
are respectively the energy of the supercell of Zr containing
k
substitutional atom of Nb, plus
n
-atoms of interstitial H, and the energy o the supercell with (
1N
) atoms of Zr
plus a substitutional Nb.
1
()
N
E Zr H
is the energy of the crystal containing
N
-atoms of Zr with a single interstitial
H at its minimum configuration, while
()
N
E Zr
is the energy of the perfect Zr lattice. Negative/positive values of
Eb
indicate repulsion/attraction between the H and the solute atom. Table 3 summarizes present results with Siesta
and Vasp.
TABLE 3: Binding energies (in eV) with Siesta/Vasp calculations for a supercell with N=48 atoms of Zr, n=1 or 2
H atoms and k=1 atom of substitutional Nb.
a
T
and
c
T
correspond to tetrahedral sites for H in the basal and axial
axes, respectively.
Zr
c
T
b
E
a
T
b
E
O
b
E
HH
0.47/0.31
0.01/-0.02
0.12/0.12
Nb H
0.11/0.29
0.08/0.06
0.05/0.04
2Nb H
-0.01/-0.02
0.31/0.44
0.01/-0.02
Our calculations predict T-sites as the preferential sites of hydrogen dissolved in
Zr, with
0.06
TO
EE
eV the
energy difference between T- and O- hydrogen sites. In presence of Nb, this energy difference was maintained.
Results in Table 3, show that
b
E
were mostly positive, indicating that the interactions are attractive, that is,
energetically favorable for the formation of the complexes, according to Domain et al..32
The exception is the binding energy of -0.01 / -0.02 eV, weak and negative, showing that the two T- hydrogen do
not interact with Nb in the axial plane. The other negative energy of -0.02 eV reveal a slight repulsion between the
two O-hydrogen first neighbors of Nb, in the basal plane.
From values in Table 3, we assume that the presence of Nb in
Zr, with a positive binding energy with H, can be
interpreted as a reduction in the hydrogen diffusivity by the hydrogen trapped locally by the Nb atom.
V. HYDROGEN DIFFUSION COEFFICIENT
Hydrogen diffusion in hcp Zr involving T and O sites driven by interstitial mechanism were studied with a model
proposed by Ishioka and Koiwa,33 using an on-lattice random walk model for hcp crystals. The diffusivity along
c
and
a
directions are expressed as,
2
(3 2 3 ) ,
4(2 3 )( 2 )
TO OO TO OO TT TT OT
c
TT TO TO OT
Dc
(8)
2.
2
TO OT
a
TO OT
Da
(9)
In Eqs. (8) and (9),
ij
are the mean jump frequencies
of the four possible hydrogen transitions, namely:
,TT
OO
,
TO
and
OT
.
In order to compute the mean jump frequencies
in Eqs. (8) and (9), we use the Vineyard formalism,34 that
corresponds to the classical limit, where the vibrational prefactors,
0
å
, do not depend on the temperature, that is
( / ).
mB
exp E k T
å
(10)
In (10),
m
E
are the vacancy migration energies at
0T
K, while
3 3 3 4
0
11
/,
m
NN
IS
ii
ii B
S
exp k
å
(11)
with
I
i
and
S
i
the frequencies of the normal vibrational modes at the initial and saddle points, respectively, and
0
the Deby’s frequency. Migration barriers were calculated with the Monomer2 coupled to Siesta,6 and the
DIMER Method3 as implemented in Vasp. Present results of the migration barriers (in eV) are summarized in
Table 4.
TABLE 4: Hydrogen migration energies,
m
E
(in eV), in
αZr
calculated with Siesta and Vasp for a supercell of
N=48 atoms of Zr and one interstitial H.
Jump
Siesta
Vasp
Ref.7
TT
0.14
0.12
0.129
TO
0.45
0.39
0.406
TT
0.36
0.32
0.346
OO
0.42
0.37
0.398
The anisotropy of hydrogen diffusion in
Zr has been widely studied from Monte Carlo Kinetic (KMC),7,35 using
individual H jumps between first neighbors (1NN) T and O sites. These authors have found an expression of the
diffusion coefficient
7
5.55 10 (0.41 / )
H
D exp eV RT
m
2
s
1
, which is higher than the experimental value of
7
(6.84 7.82) 10 ( 0.46 0.47 / )
H
D exp eV RT
m
2
s
1
36 for
1000T
K.
On the other hand, Sawatzky et al.35 have obtained the diffusion coefficient for hydrogen in annealed and extruded
specimens as
7
1.17 10
H
D
m
2
s
1
with
33.60
D
Q
KJ/mol over the temperature range 473 to 973 K from
KMC calculations. While, Zhang et al.,7 are who have obtained results closer to the experimental ones. Fig. 3
shows our results of
c
D
and
a
D
in
Zr
calculated from Eqs. (8) and (9) using both, Siesta and Vasp codes.
Experimental data from Ref.37 and previous numerical results from KMC calculations in Ref.,35 are also included.
Experimental results of hydrogen diffusion in
Zr-2.5Nb,38 are shown in blue line.
The mean jump frequencies were calculated from Eqs. (10) and (11) with the migration barriers listed in Table 4,
verify
1 3 4 2
.
Fig. 3(c) shows the behavior of the mean jump frequencies with the inverse of the temperature.
FIG. 3: Hydrogen diffusion coefficients vs 1/T (K−1). The hydrogen diffusion coefficients Dc and Da in black solid
and dashed lines respectively, from present calculations with: (a) Siesta, and (b) Vasp. In open triangles and stars
DH from experimental data 36 and KMC results 35. The experimental hydrogen diffusion coefficient in Zr-2.5Nb 38
(in blue solid line) is also shown. (c) Mean jump frequencies vs 1/T (K−1).
1 3 4 2.
Fig. 3(a)-(c), show our numerical results of
H
D
and the mean jump frequencies in
Zr in terms of 1/T (in K
1
),
compared with experimental data in pure metal and the alloy Zr-2.5 Nb. Our results with Vasp reproduce very well
the experimental data and the numerical results of KMC.
The anisotropy of the H diffusion in hcp Zr along the
c
and
a
axes, is measurable from the
/
ca
DD
ratio,
2
2
22
3( ).
8 3 2 3
c OO TT
a OT TT
Dc
D a TO
(12)
The dependence on the temperature of the ratio
/
ca
DD
, defined in Eq. (12) is through the mean jump frequencies
i
. Fig. 4 shows that the results of the ratio
/
ca
DD
in Ref.7 (symbols in orange) agree better with present results
with Vasp (in dashed line) than with Siesta (solid line). Possibly because Siesta overestimate the values of the
migration barriers.
FIG. 4:
/
ca
DD
vs
1/ T
(in
1
K
). Symbols in black show our results with Vasp (open stars) and Siesta (full stars). In
orange the results in Ref. 7
VI. ZIRCONIUM HYDRIDES
The three most stables zirconium-hydrides in Zr based alloys are
,
,
. The
hydride is the most frequent in
Zr-2.5Nb which appears at slow speed of cooling low than
2
C/ min. The composition of the
-hydride is
x
ZrH
with
1.53 1.66x
(56.7 to 66 % at Hydrogen, simulated as
1.5
ZrH
). This hydride forms plates located at the
grain boundaries or trans-granular interphases.
This particular form is associated with the following relation of orientation:
(0002) (111)
and
[1120] [110]
,
which generates a
(1017)
plane that is practically parallel to the basal plane of Zr39 and precipitates mostly in the
radial direction, as shown in Fig. 5.40
FIG. 5: Hydride’s Radial and circumferential orientation in tubes of Zr-2.5Nb.
We present our preliminary results in Table 5, of the three stable hydrides schematized in Fig. 6. Our results
reproduce well the ones obtained by Glasoff.41 We have followed the same methodology to familiarize ourselves
with this type of structures. The equilibrium lattice constants were determined by lattice relaxation, minimizing the
total energy of the system by varying atomic positions with the calculation parameters defined in Section II.
FIG. 6: Crystal structure of the 3 stable hydrides in zirconium: orthorrombic
γ ZrH
; cubic
1.5
δ ZrH
and
tetragonal
2
ε ZrH
.
In Table 5, the values of the lattice parameter of each cell vary little from one hydride to another. Our results agree
very well with the calculations made by Glazoff41 and with the available experimental data. The hydride formation
enthalpy,
x
f
E
, is calculated as Glasoff.41
TABLE 5: Lattice parameters (in Å) and the hydrogen solution energy (in eV) of the three stables hydrides in
αZr
with Vasp.
System
ZrH
1.5
ZrH
2
ZrH
a
3.232
4.807
3.538
a
13
3.223
4.482
3.537
exp
a
42
3.243
4.777
3.520
c
5.016
4.807
4.406
c
13
5.014
4.821
4.458
exp
c
42
4.948
5.777
4.449
/ca
1.55
1.0
1.25
/|
exp
ca
42
1.53
1.0
1.26
x
f
E
-0.408
-0.564
-0.573
x
f
E
13
-0.32
-0.548
-0.556
VII. CONCLUSIONS
In summary, we have studied the hydrogen diffusion by ab initio methods in
Zr. Also, the energetic of hydrogen
in the Zr-Nb system were also studied. We remark that,
• our DFT results of the lattice parameters are in well agreement with the experimental values.
• The calculated formation energies are within the range of values of previous references and experimental
values.
• The mixing energy in the dilute limit shows that Zr is energetically favorable to dissolve in the bcc phase of
Nb, and on the contrary, Nb has a low solubility in hcp Zr.
• The presence of Nb in hcp Zr moderately increases the solubility of H around the impurity, implying in a
slight reduction in diffusivity of the H by trapping H locally, forming the H-Nb complex.
• The hydrogen solution energy shows that tetrahedral sites are the most stable for H.
• The H-diffusion coefficient in hcp Zr, calculated with Vasp reproduce very well the experimental values and
previous results from KMC.
ACKNOWLEDGEMENTS
We want to thank Dr. Gladys Domizzi for her comments on the manuscript. This work was partially financed by
CONICET PIP-11220170100021CO.
REFERENCES
[1] G. Beltramo, V. Ramunni y R. Pasianot. Theoretical and numerical study of hydrogen diffusion in Zr-Nb
alloys. Anales AFA 32, 112-119 (2021).
[2] V. P. Ramunni, M. A. Alurralde y R. C. Pasianot. Search of point defect transition states in hcp twin
boundaries: The monomer method. Phys. Rev. B 74, 054113 (2006).
[3] G. Henkelman y H. Jonsson. A dimer method for finding saddle points on high dimensional potential surfaces
using only first derivatives. J. Chem. Phys. 111, 7010-7022 (1999).
[4] G. Kresse y J. Hafner. Ab initiomolecular-dynamics simulation of the liquid-metal–amorphous-semiconductor
transition in germanium. Phys. Rev. B 49, 14251-14269 (1994).
[5] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon y D. Sanchez-Portal. The SIESTA
method forab initioorder-Nmaterials simulation. J. Phys. Condens. Matter 14, 2745-2779 (2002).
[6] R. Pasianot, R. Perez, V. Ramunni y M. Weissmann. Ab initio approach to the effect of Fe on the diffusion in
hcp Zr II: The energy barriers. J. Nucl. Mater. 392, 100-104 (2009).
[7] Y. Zhang, C. Jiang y X. Bai. Anisotropic hydrogen diffusion in α-Zr and Zircaloy predicted by accelerated
kinetic Monte Carlo simulations. Sci. Rep. 7, 41033 (2017).
[8] F. Birch. Finite Elastic Strain of Cubic Crystals. Phys. Rev. 71, 809-824 (1947).
[9] J. P. Perdew, K. Burke y M. Ernzerhof. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett.
77, 3865-3868 (1996).
[10] P. E. Blochl, O. Jepsen y O. K. Andersen. Improved tetrahedron method for Brillouin-zone integrations. Phys.
Rev. B 49, 16223-16233 (1994).
[11] H. J. Monkhorst y J. D. Pack. Special points for Brillouinzone integrations. Phys. Rev. B 13, 5188-5192
(1976).
[12] R. Pasianot y A. Monti. A many body potential for α-Zr. Application to defect properties. J. Nucl. Mater. 264,
198-205 (1999).
[13] M. Christensen, W. Wolf, C. Freeman, E. Wimmer, R. B. Adamson, L. Hallstadius, P. E. Cantonwine y E. V.
Mader. H inα-Zr and in zirconium hydrides: solubility, effect on dimensional changes, and the role of defects. J.
Phys. Condens. Matter 27, 025402 (2015).
[14] D. R. Lide. CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data
(CRC press, Boca Raton, Florida, 1995).
[15] J. Goldak, L. T. Lloyd y C. S. Barrett. Lattice Parameters, Thermal Expansions, and Gruneisen Coefficients of
Zirconium, 4.2 to 1130°K. Phys. Rev. 144, 478-484 (1966).
[16] C. Kittel. Introduction to Solid State Physics (Wiley, New York, 2005).
[17] H. H. Neely. Damage rate and recovery measurements on zirconium after electron irradiation at low
temperatures. Radiation Effects 3, 189-201 (1970).
[18] G. Hood y R. Schultz. Defect recovery in electronirradiated α-Zr single crystals: A positron annihilation
study. J. Nucl. Mater. 151, 172-180 (1988).
[19] M. Cottura y E. Clouet. Solubility in Zr-Nb alloys from first-principles. Acta Mater. 144, 21-30 (2018).
[20] D.-Y. Lin, S. S.Wang, D. L. Peng, M. Li y X. D. Hui. Annbody potential for a Zr–Nb system based on the
embeddedatom method. J. Phys. Condens. Matter 25, 105404 (2013).
[21] R. Roberge. Lattice parameter of niobium between 4.2 and 300 K. Journal of the Less Common Metals 40,
161-164 (1975).
[22] P. Soderlind, L. H. Yang, J. A. Moriarty y J. M. Wills. First-principles formation energies of monovacancies
in bcc transition metals. Phys. Rev. B 61, 2579-2586 (2000).
[23] M. I. Baskes. Modified embedded-atom potentials for cubic materials and impurities. Phys. Rev. B 46, 2727-
2742 (1992).
[24] H. Ullmaier, P. Ehrhart, P. Jung y H. Schultz. Atomic defects in metals (Springer, Berlin, 1991).
[25] V. O. Kharchenko y D. O. Kharchenko. Ab-initio calculations for the structural properties of Zr-Nb alloys.
Condens. Matter Phys. 16, 13801 (2013).
[26] C. Domain. Ab initio modelling of defect properties with substitutional and interstitials elements in steels and
Zr alloys. J. Nucl. Mater. 351, 1-19 (2006).
[27] X. Xin, W. Lai y B. Liu. Point defect properties in hcp and bcc Zr with trace solute Nb revealed by ab initio
calculations. J. Nucl. Mater. 393, 197-202 (2009).
[28] D. Jiang y E. A. Carter. First principles assessment of ideal fracture energies of materials with mobile
impurities: implications for hydrogen embrittlement of metals. Acta Mater. 52, 4801-4807 (2004).
[29] K. Huber y G. Herzberg. Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules
(Van Norstrand Reinhold Co, 1979).
[30] C. Varvenne, O. Mackain, L. Proville y E. Clouet. Hydrogen and vacancy clustering in zirconium. Acta Mater.
102, 56-69 (2016).
[31] K. Lee, M. Yuan y J. Wilcox. Understanding Deviations in Hydrogen Solubility Predictions in Transition
Metals through First-Principles Calculations. J. Phys. Chem. C 119, 19642-19653 (2015).
[32] C. Domain, R. Besson y A. Legris. Atomic-scale Ab-initio study of the Zr-H system: I. Bulk properties. Acta
Mater. 50, 3513-3526 (2002).
[33] S. Ishioka y M. Koiwa. Diffusion coefficient in crystals with multiple jump frequencies. Philos. Mag. A 52,
267-277 (1985).
[34] G. H. Vineyard. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids 3,
121-127 (1957).
[35] A. Sawatzky, A. Ledoux, R. Tough y C. Cann. Miami International Symposium: April 1981 109-120 (Oxford:
Pergamon Press, Miami Beach, Florida, 1982).
[36] J. Kearns. Diffusion coefficient of hydrogen in alpha zirconium, Zircaloy-2 and Zircaloy-4. J. Nucl. Mater. 43,
330-338 (1972).
[37] C.-S. Zhang, B. Li y P. Norton. The study of hydrogen segregation on Zr(0001) and Zr(1010) surfaces by
static secondary ion mass spectroscopy, work function, Auger electron spectroscopy and nuclear reaction analysis.
J. Alloys Compd. 231, 354-363 (1995).
[38] G. McRae, C. Coleman, H. Nordin, B. Leitch y S. Hanlon. Diffusivity of hydrogen isotopes in the alpha phase
of zirconium alloys interpreted with the Einstein flux equation. J. Nucl. Mater. 510, 337-347 (2018).
[39] D. Northwood y R. Gilbert. Hydrides in zirconium-2.5 wt.% niobium alloy pressure tubing. J. Nucl. Mater.
78, 112-116 (1978).
[40] M. Leger, G. Moan, A. Wallace y N. Watson. Zirconium in the Nuclear Industry: 8 th Int. Symposium (ASTM
STP 1023, Philadelphia, 1979).
[41] M. V. Glazoff. Modeling of Some Physical Properties of Zirconium Alloys for Nuclear Applications in
Support of UFD Campaign inf. tec. (2013).
[42] S. S. Sidhu, N. S. S. Murthy, F. P. Campos y D. D. Zauberis. en Advances in Chemistry 87-98 (1963).
