Nuclear Inst. and Methods in Physics Research B Modelling the effect of hydrogen on crack growth in zirconium

Via molecular dynamics simulations, the effects of hydrogen on stress evolution of -zirconium and crack propagation in monocrystalline and multiple grained zirconium systems are investigated. Diffusion barriers are shown to reduce when strain is applied, which then causes hydrogen to accumulate at surfaces and grain boundaries. Crack growth is considered for a range of -zirconium systems, both with and without hydrogen, strained in multiple directions. The effects of crystal orientation are shown to be of high influence on the stress evolution of -zirconium irrespective of hydrogen content. Crack growth velocity is increased the most by hydrogen for -zirconium when uniaxial strain is applied in the [0001] direction. Simulations are conducted investigating the effects of single grain boundaries in normal and parallel orientations to crack growth showing a high importance on the location of interstitial hydrogen in crack growth behaviour. In addition, larger scale simulations show the effects of multiple grain boundaries and hydrogen content in the evolution of cracks.


Introduction
Zirconium (Zr) and its alloys are extensively used in nuclear reactors due to their low neutron absorption cross-section and good corrosion and fracture resistance. Their main application is for nuclear fuel cladding in boiling water reactors (BWRs) and pressurised water reactors (PWRs) [1] but in CANada Deuterium Uranium (CANDU) nuclear reactors a specific zirconium-niobium alloy (Zr-2.5%Nb) is extensively used in pressure tubes [2]. Each pressure tube contains a fuel bundle and a continuous flow of coolant when in operation. Zirconium alloys used in fuel cladding and pressure tubes are susceptible to hydrogen (H) embrittlement and hydride cracking during their service life in a nuclear reactor. While Zr-2.5%Nb is designed to be less susceptible to hydride embrittlement than other alloys, its use in pressure tubes, with an expected service life in the order of 30 years [3], means that hydride cracking is still a concern [4]. Hydride cracking can also be a mechanism of failure outside of operating conditions. Fuel cladding tubes can fail during long-term dry storage at room temperature [5]. The presence of hydrides in zirconium and zirconium alloy pressure tubes is associated with initial crack formation and subsequent crack propagation [6]. In operating conditions, hydrides can form in regions of sufficient tensile stress caused by cold spots, service induced damage and manufacturing flaws. Hydrides dramatically reduce the ultimate tensile strength of Zr and can cause cracks to form and quickly propagate. This crack mechanism is known as delayed hydride cracking (DHC). Hydrogen ingress occurs when the coolant (water or heavy water) reacts with the zirconium or zirconium alloy tubes to form an oxidised layer with the release of hydrogen: There are numerous models in the literature used to predict the delayed hydride crack velocity [7][8][9][10] but little is known about the mechanical properties of Zr and zirconium based alloys with H contaminates at the atomistic scale. Previous works have shown successful atomistic scale simulations investigating crack propagation in iron [11][12][13], aluminium [14], nickel [15,16] and silver [17]. A similar approach is considered in this work and is concentrated on the pure Zr material rather than individual zirconium based alloys. Structural properties such as average grain sizes, orientation of grains and local positioning of H within Zr may be important factors for DHC occurring and the velocity at which the crack propagates. An investigation into the best structural properties of Zr for DHC resistance may help dictate manufacturing processes for improved pressure tubes. This paper considers an atomistic approach of investigation into the effects of H concentrations and positioning along with the effects of various grain sizes and orientations of grain boundaries have on crack propagation in Zr. The diffusion of hydrogen in the Zr crystal and both in and close to grain boundaries is also considered. The method used is classical molecular dynamics (MD) which has previously been used to model radiation damage and crack propagation in strained Zr [18,19].

Methodology
The MD simulations are carried out using a second nearest neighbour (2NN) modified embedded atom method (MEAM) [20] to model the atomic interactions. To model all Zr-H interaction an existing MEAM potential for pure Zr [21], a slightly modified MEAM fitting for H interactions [22] and a fit for the binary Zr-H system [23] were used. This model has been shown to reproduce various structural and elastic properties for -Zr (hexagonal close packed, HCP), lattice constants, cohesive energy, dilute heat of solution of hydrogen, vacancy formation energy and planar defect properties in -ZrH 2 as well as diffusion coefficients of H in bulk Zr in good agreement with ab initio and experimental data by Kim et al. [21] and Lee et al. [23]. Whilst the diffusion coefficients and the preference for diffusion in the [0 0 0 1] direction of H in Zr agrees well with experimental results, this is at the expense of the potential overestimating the stability of the octahedral interstitial site compared to the ab initio predicted preferred tetrahedral site. Despite this, the capability of the MEAM potential to correctly predict mechanical properties means it is used in this work to describe crack growth. The model has previously been used to investigate stressstrain curves of -Zr with various concentrations of H when strained in the [0 0 0 1] direction and the effect of -ZrH 2 precipitates inside bulk Zr systems [19].
All molecular dynamics simulations were performed via the Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [24] on Loughborough University's own high performance computing (HPC) system. HPC systems are essential for large molecular dynamics simulations: the larger systems in the paper took 20,000 core hours to simulate less than 1 ns worth of crack growth. An isothermal-isobaric (NPT) ensemble is considered with Nosé-Hoover thermostat and barostat [25] for all MD simulations. For systems containing H, a time integration step of 0.1 femtoseconds (fs) is required for energy conservation due to the fast vibrations of H atoms whereas a larger timestep of 1 fs is adequate for pure Zr simulations. The system stress, , considered for stress-strain curves calculations is the virial definition (i.e. no kinetic energy contribution): In Eq. (2), f ij is the force between atoms i and j in the direction, r ij is the distance between i and j in the direction, V is the volume of the system cell, N is the total number of atoms and M is the total number of neighbouring atoms to atom i. Visualisations of the atomic configurations are produced via the Open Visualisation Tool (OVITO) [26] and the adaptive common neighbour analysis (ACNA) [27] is used for structural analysis. In crack growth simulations, the crack tip depths are plotted against simulation run time and the gradient of a linear line of best fit is used to find the average crack velocity.

Effect of strain on -Zr
In a previous work by Chakraborty et al. [19], large scale MD simulations were used to investigate the effect of H on the stress-strain relationship in Zr when strained in the [0 0 0 1] direction only. Our research extends the work performed by Chakraborty et al. and considers the effect of strain in multiple directions. The system used for stress-strain curve calculations contained 12,000 atoms initially in a cubic cell ( × × 6.5 6. The results from our simulations show a significant difference in deformation behaviour and ultimate tensile strength (UTS) for pure -Zr strained in various directions. Fig. 1 shows the effect on Zr crystal structure and stress when uniaxial strain is applied. When a strain is applied in the [2 1 1 0] or < > a direction of the HCP structure, the bulk material spontaneously changes phase at 0.15 strain to the bodycentred cubic (BCC) -Zr phase and results in the lowest UTS of the 3 orientations considered. After further strain is applied, dislocations occur in the BCC structure in order to reduce the virial stress. For strain applied in the [0 0 0 1] direction, elastic deformation occurs until the strain reaches 0.16 and then dislocations initiate and propagate throughout the bulk. The last direction considered, [0 1 1 0], is perpendicular to the first 2. This direction of strain allows a large build up of stress and elastic deformation until a critical level of strain is reached. Once the strain exceeds 0.43, the whole of the originally HCP crystal spontaneously generates dislocations throughout the structure. The phase diagram of Zr has been investigated both theoretically [28] and experimentally [29,30]. At high temperature Zr adopts the BCC phase and at high compressive pressure, the (hexagonal) phase. In Fig. 1 there is a change in slope of the stress strain curve in the [0 1 1 0] direction at a strain of 3. We could not detect that this corresponded to a phase change such as to the structure which is known to occur under small compressive stress. For strain in the [2 1 1 0], where there is a change to BCC, we have not found experimental evidence to support this finding although Zr is known to adopt the BCC phase at higher temperature [30].
To investigate the effect of H in Zr, H atoms were randomly positioned in preferred octahedral interstitial sites within a bulk -Zr system and strained in the < > c and < > a directions (or [0 0 0 1] and [2 1 1 0] respectively) of the HCP crystal: the two directions with the lowest UTSs. Various concentrations of H were considered: pure Zr, 1% H, 2% H and 3% H in pure Zr. Note that all concentrations are given in atomic percent and 1at% H in Zr is equivalent to 110 wt ppm. Fig. 2 shows the effect of H concentrations on the stress-strain curves. The results in Fig. 2(a) show that the effect of H is to reduce the value of the strain at which the phase transformation occurs, shown in Fig. 1 and in Fig. 2(b) the generation of dislocations is enhanced. So the higher the concentration of H, the lower the UTS and strain required to permanently deform the material -a trend seen in the work by Chakraborty et al. The MD simulations show that the H remains isolated and does not cluster in the material and the effect H has on the strength of Zr appears to reduce once higher levels of H are introduced; the stressstrain relation is effected the most by the initial inclusion of 1 % H.
The values calculated for UTS for the various systems differ largely from experimental data [31,32] due to the small system sizes and extremely high strain rates required for the MD simulations. Smaller system simulations were also considered with a slower strain rate of 10 7 s −1 for pure Zr strained in various directions and a similar trend is seen albeit with a lower UTS for each case. By combining these results and the results in this work to those by Chakrabrty et al. it can be seen that by increasing the system size and reducing the strain rate, a reduction in calculated UTS can be seen. Despite this, MD simulations can predict the overall trend of mechanical behaviours well and give a good understanding of the effects of H on the strength of Zr.

H diffusion in Zr
Whilst stress-strain curves can vary with system size and strain rates, properties that are not time nor particularly size dependent include the energy barriers for H diffusion in crystalline Zr. Strain especially affects the H diffusion barriers in the preferred diffusion direction: [0 0 0 1]. Transition energy barriers are calculated via the climbing image nudged elastic band [33] technique. The minimum energy pathway band along with the initial and final states are relaxed with a force minimisation tolerance of 0.05 eV/Å via the conjugate gradient method [34].  (Table 1). The barrier for H to move onto a perfect -Zr (0 0 0 1) surface from a subsurface octahedral site is also reduced from 0.7 eV to 0.3 eV when 5% strain is applied. These calculated barriers are in broad agreement with experimental measurements of H diffusion in Zr [35,36]. The reverse barrier is also reduced from 2.5 eV to 2 eV but is still too large for H to penetrate the perfect Zr surface even at high temperature, suggesting that ingress to the material would be directly through grain boundaries or other defects.
Similar calculations near to different types of grain boundary show a range of transition energy barriers ranging from 0.2 to 0.5 eV in Zr without strain with a corresponding reduction when strain is applied although a full statistical analysis has not been carried out due the large number of potential configurations that would be required to be investigated.
To understand the grain boundary effect, MD simulations were run at 600 K by randomly distributing interstitial H near a tilt grain boundary in pure Zr without strain. PBCs were applied in directions parallel to the grain boundary and atoms at the surfaces were fixed. Fig. 3(a) shows an initial structure with a 65°tilt grain boundary and 1% H in a system consisting of 3200 atoms. Fig. 3(b) shows the diffusion of H within the Zr structure. The preferential diffusion direction   [13][14][15][16][17][18][19][20] above and below the grain boundary is clearly seen -corresponding to the crystalline < > c direction. There is also an occasional sideways hop. Also shown is an increased movement of H at the boundary itself and the attraction of H towards it in agreement with experiment [37,38]. A rough calculation estimated that the diffusion of H near to this grain boundary is 70% faster than in the bulk crystal.
These results therefore show (1) that H will accumulate at surfaces with diffusion there being enhanced by strain. (2) H can enter the material through grain boundaries and (3) even if H were to be present in the bulk, diffusion to the grain boundaries will occur.

Crack propagation in monocrystalline -Zr
The effect of interstitial H and orientation of applied strain on crack propagation in Zr was next investigated. The systems considered have a surface with an initial crack inserted perpendicular to the surface. To model crack propagation, a simplified effective '2D' system is considered. In one direction (x-direction) the system width of only 10 Å is considered. This width is larger than the potential cut-off used in the 2NN MEAM to prevent atoms self-interacting through PBCs. PBCs are applied in the thin, x, direction to simulate a 3D structure whilst the crack propagation, y, and strain, z, directions are perpendicular to this. An example system schematic is given in Fig. 4. The outside 2-3 Å of atoms in the z direction and bottom of the system are fixed relative to the simulation box. Before strain is applied, the system is relaxed at 1 K for 20 ps in an NPT ensemble to ensure a zero-stress condition in all directions. The simulation box is then extended in the strain direction and the relative atom positions are updated each time step according to the strain rate. A zero-stress condition is kept for the x-direction during the application of strain.
Three different orientations of monocrystaline -Zr were considered for crack growth simulations. The Zr system was positioned such that the strain directions for the three various orientations corresponded to two perpendicular directions, [ Overall trends of how H effects crack propagation have a dependence on the orientation of the crystal when being strained. In Fig. 5, the crack tip position is plotted against the simulation time for the three different orientations. In each case, systems with and without H are compared. For uniaxial strain applied in the [0 0 0 1] and the pyramidal slip direction, the crack tip position over time varies little depending on whether 1% H is introduced or not -see Fig. 5(b) and 5(c). However, a more marked effect is seen for the other strain direction. For the [2 1 1 0] case, including 1% H in the system reduces the crack growth rate by 30% -see Fig. 5(a). It is also interesting to note that a strain of 0.2 corresponds to phase changes and dislocation emission in bulk Zr and also the value at which the stress is suddenly reduced in the crack propagation case as shown in Fig. 5.
Also note the overall behaviour of the crack tip position as the system becomes strained; in 1 of the 5 [1 1 2 3] cases investigated, without H, a void appeared directly beneath the crack tip, causing a sudden increase in crack depth as the crack tip met the expanding void. This effect can be seen in the plot and explains the sudden increased length of the error bars in Fig. 5(c).
In all monocrystalline crack simulations considered, the fracture strength varied little depending on orientation (see Fig. 6); ranging between 5.3 and 5.8 GPa without H and between 4.8 and 5.5 GPa with H. In agreement with stress-strain analysis done on bulk -Zr, the strength was reduced slightly with the presence of interstitial H in each system explored.

Crack propagation in -Zr with grain boundaries
Materials used for pressure tubes and fuel cladding are not perfectly crystalline but instead contain many different grains and defects. Here we have considered two different sized systems: two smaller systems with only 2 grains and two larger systems with 6 grains. The smaller systems allowed us to focus on how the orientation of an isolated single grain boundary (GB) can affect the propagation of a crack when it comes into contact with it. The larger system then simulates a more realistic scenario with 6 grains of the same orientation generated by random rotation about the x-axis.
The multi-grained,effectively two-dimensional systems, whose sizes  are shown in the captions to Figs. 7 and 10 were produced by randomly placing N (number of grains) points within a system, creating Voronoi cells [39] for each point and filling the Voronoi cells with randomly rotated -Zr crystal. The atoms in the Voronoi cells are randomly rotated about the x-axis and if any atoms lies within 2.4 Å of another, one of them is removed. The whole system is then heated to 600 K and then relaxed to 1 K with zero-stress PBCs before a surface and initial crack is created.
Since our model suggests that H is likely to accumulate at grain boundaries we also consider crack propagation when H is accumulated near the grain boundaries and when H randomly distributed throughout the system. For both the 2 grain and 6 grain systems, 1% H is added to system in random octahedral interstitial sites. In the other case 1% H is concentrated within 15 Å of the grain boundaries.

System with 2 grains
For the 2 grain systems, two orientations of grain boundaries are considered. A vertical boundary that is positioned directly along the expected crack growth direction and a horizontal boundary normal to the expected crack growth direction positioned 3 nm below the initial crack tip (shown in Fig. 7(c)). The system is strained in the z axis at a constant rate of 10 9 s −1 at 1 K. Fig. 7(a) shows the resulting crack propagation, comparing 1% randomly positioned H and 1% H concentrated at the horizontal grain boundary whilst Fig. 7(b) compares 1% randomly positioned H and 1% H concentrated at the vertical grain boundary.
In the horizontal GB case, the high concentration of H near the grain boundary accelerates crack growth. However the crack propagation plateaus once the crack tip reaches the horizontal GB in both cases but is allowed to deform more with higher concentrations of H. For the vertical GB case, the evolution of the crack tip does not appear to vary depending on the positioning of interstitial H atoms.
Overall, there is a slight dependence on the positioning of interstitial H atoms in a multi-grained Zr system. However, the largest effect on crack propagation is in fact the orientation of the GBs themselves. GBs positioned perpendicular to crack growth may slow down crack propagation.

Systems with 6 grains
Large systems containing 340,000 atoms arranged into 6 grains, with the same orientation randomly rotated about the x-axis, with an average size of 1200 nm 2 , are then considered. The systems are strained in the same way as the 2 grain cases. By straining these large multigrained systems a better understanding of the effect of H on crack growth in a real system can be obtained. In Figs. 8 and 10 two systems are set up with 6 randomly orientated grains. Only two scenarios are considered here due to the large amount of computational power required to run these simulations. In both figures, the initially relaxed multi-grained system and crack (20 Å deep) is shown and compared to three different scenarios after 400 ps of simulation (equivalent to 0.4 strain). The three scenarios include (a): pure Zr, (b): Zr with 1% H randomly positioned and (c): Zr with 1% H positioned concentrated near the grain boundaries.
The crack depth is measured for all scenarios considered and plotted  against time in Figs. 9 and 11. The figures show that the crack propagates at approximately the same speed after 300 ps (0.3 strain) in all cases but that the crack depth is enhanced in the case of H near the grain boundary and in one case reduced when H is uniformly distributed. Although the statistics are not good enough to draw firm conclusions this is in agreement with the results of Fig. 5 where the effect of uniformly distributed H slightly reduces the crack depth, in one orientation but the ultimate speed at sufficiently large strain is about the same in all orientations. However in both multi-grained setups, the high concentration of H atoms around GBs has a marked effect on crack depth. Each case shows an increase in crack depth for a given strain, suggesting that Zr with hydrides at grain boundaries are the most susceptible to cracking.
In Fig. 10, multiple voids are shown to appear for the pure Zr case and a single large void arises where 3 grains meet when H is concentrated at the grain boundaries. However, in the 6 grain setup, no voids are seen when 1% H is randomly positioned; suggesting that if were to be possible to arrange small amounts of dissolved H in Zr a more ductile material would be created that can better withstand strain.

Discussion and conclusions
Results from all the various MD simulations show the importance of interstitial H, the position of the H interstitials and orientations of crystal structures in stress evolution and crack growth mechanisms. Hydrogen is shown to diffuse preferentially to grain boundaries and exposed Zr surfaces and can enhance crack propagation when it does (in agreement with experiment). Energy barrier calculations indicate that H enters the material through defects and grain boundaries.
The results show that in a perfect crystal the strain direction is important in the development of cracks and that in certain directions the presence of hydrogen can speed up the propagation of cracks. In a nanocrystalline system voids can develop around the triple junctions   where grains of different orientations meet but from the limited statistics available, no significant difference was observed between the cases when no H was present compared to when 1% H was randomly positioned. However, there is evidence to suggest that the location of H near grain boundaries in multi-grained systems does have a significant effect on the crack growth.
Although the strain rates are faster than those likely to be encountered in real pressure tubes and the effect of an oxidised surface layer was not considered, the results do give some insight into the role of H within the system. Further work could utilise long time scale techniques to access more realistic strain rates and develop models for the effect of the oxidised surface.

Data availability
The raw data required to reproduce these findings are available to download from DOI:https://doi.org/10.17028/rd.lboro.7151912 andhttps://doi.org/10.17028/rd.lboro.7151930. Fig. 10. Initial crack setup with 6 grains and crack growth propagation comparison of pure Zr, Zr with 1% H randomly positioned and Zr with 1% H concentrated at grain boundaries after 400 ps (0.4 strain). Initial system dimensions × × 1 100 70 nm 3 . Atoms are coloured by initial grain orientation with average grain size of 1200 nm 2 . Fig. 11. Plot comparing crack depth versus time for -Zr system with 6 grains shown in Fig. 10.