Modelling Hydrogen Embrittlement using Density Functional Theory: A theoretical approach to understanding environmentally assisted cracking in 7xxx series aluminium alloys

The effects of H segregation to a Σ11 symmetric tilt Al grain boundary are investigated using atomistic simulations, as part of a wider study on cracking in 7xxx series alloys. Density functional theory based simulations of uniaxial straining of grain boundaries containing 11 different concentrations of H were performed under the cohesive zone fracture mechanics framework. The theoretical strength of grain boundaries is shown to be supressed by H segregation, and the cause of this is attributed to the prevention of the formation of Al ligaments across grain boundaries. Segregated concentrations of relevant alloying elements (Zn, Mg, and Cu) show minimal impact on the H embrittlement process investigated, namely H enhanced decohesion (HEDE). Further modelling, of H transport and grain boundary precipitates, is required to confirm the validity of the HEDE mechanism in the case of 7xxx alloys.


Introduction
7xxx series aluminium alloys, Al-Zn-Mg-(Cu), are used as high performance structural materials -AA7075 for instance carries an ultimate tensile strength 572MPa in the peak aged T6 temper [1]. However, material degradation arising from a collection of phenomena generally labelled as environmentally assisted cracking (EAC) has been a key challenge in the development of these alloy systems [2]; for alloys used by the aircraft industry, EAC is an intergranular phenomenon [3]. Considerable study within metallurgy has delivered two fundamental causes of EAC, namely anodic dissolution (AD) and hydrogen embrittlement (HE), and the problem has been recorded in a variety of environments where the relevant atomospheric conditions are reportedly humidity and chloride content [4].
Investigationg HE processes in Al alloys is a challenging endeavour because of the low equilibrium lattice solubility of H, and Al's inability to form stable hydrides [2]. Whilst evidence for H inducing cracks is limited, partitioning of H has been shown to degrade the material's resistance to crack growth. Detection of H's influence on localised effects has recently been achieved through advanced imaging techniques such as X-ray synchotron tomography [5] which confirms significant H enrichment at grain boundaries (GBs) and precipitate interfaces during deformation. A theoretical approach to understanding HE can circumvent experimental challenges by predicting material behaviour using fundamental physical and chemical principles [6] -instead of probing a material through successive experimental techniques, this work aims to model * e-mail: benjamin.wilson@manchester.ac.uk * * e-mail: joseph.robson@manchester.ac.uk * * * e-mail: christopher.race@manchester.ac.uk from first principles the relationship between intergranular fracture and H in 7xxx alloys by using electronic structure calculations of the cohesive strength of grain boundaries.

Hydrogen embrittlement
Proposed HE mechanisms are often speculative due to the difficulty in quantifying hydrogen levels accurately, although four main mechanisms are often discussed in the context of Al alloys. Hydriding has been detected in some Mg containing phases of aluminium and could potentially lead to brittle GBs [7]. Hydrogen-enhanced local plasticity (HELP) is a complex phenomenon involving high localised H concentrations increasing dislocation mobility by shielding dislocation interactions. Another related model of HE involves adsorbed H (i.e. external H) inducing the emission of dislocations [8]. The fourth mechanism is hydrogen-enhanced decohesion (HEDE) which posits that segregated H at the GB compromises cohesive forces between grains by disrupting bonds across the interface [9].
The HEDE mechanism itself requires several distinct steps [2]. H must be produced in large enough quantities to impact crack growth in any measurable way -a bare Al crack tip in the presence of water vapour will oxidise to provide a source of atomic H, but passivation will quickly inhibit this process. H transport along the GB has a low flux, meaning that under the HEDE model H is absorbed into the metal matrix and diffuses to the highly stressed region just ahead of the crack tip, Figure 1 shows a schematic of this process. If H in this region substantially embrittles the GB then the crack will advance, gen-  H is generated on fresh crack surfaces and is transported to the high stress region ahead of the crack tip, which leads to H embrittlement induced crack growth. The cohesive zone (CZ) is a small region immediately ahead of the crack tip and is simulated in this work using density functional theory. erating fresh surfaces for H production and the process is repeated.

Cohesive zone fracture mechanics
For fracture events displaying limited plasticity, such as those found in EAC, cohesive zone (CZ) modelling is a useful tool which introduces a virtual region at the crack tip called a CZ [10]; Figure 1 shows a CZ in a HEDE model. The CZ captures separating atomic planes immediately ahead of the crack tip and crack growth is seen as the successive failures of CZs along the crack path. First-principles simulations can be used to model the separation of atomic planes, or decohesion, as the stretching and breaking of molecular bonds which exist across the crack path [11]. In the context of HEDE, a CZ will be weaker when decorated with H. The CZ model assumes linearly elastic behaviour in the regions surrounding a CZ. A macroscopic load on the material will translate into a microscopic traction τ(u) which, in the CZ approximation, is a unixial stress applied perpendicular to the GB plane; the separation of two grains can be quantifed by a displacement u. The energy of decohesion as a function of displacement χ(u) is evaluated in (1).
A CZ fails when the bonding across its interface is completely broken i.e. the two surfaces are displaced past their interaction range u f , with the work of separation (the energy required to complete the decohesion process) defined as W = χ(u f ). Modern electronic structure calculations, or ab-initio simulations, avoid the requirement for difficult to measure experimental parameters as an input, and instead rely on fundamental quantum mechanics based principles to predict the strength performance of materials, and have been effective in calculating the traction separation laws of CZ models [12]. These calculations are routinely performed using Density Functional Theory (DFT) which approximates electron interactions through the use of appropriate exchange-correlation potentials. Application of DFT reflects a trade-off between fundamental complexity of atomistic behaviour and computational feasibility -finite temperature effects are often ignored in these simulations [13]. DFT based ab initio tensile test simulations are used to obtain W, u f , the maximum traction τ * required for fracture, and the displacement u * at τ * , properties which are collectively referred to as a tractionseparation law for a CZ [14]. Once traction-separation laws have been evaluated they can be implemented in larger scale continuum models to constitute a full CZ based fracture mechanics model [15].

Electronic structure calculations
We use the Vienna Ab-initio simulation package (VASP) [16][17][18][19] with the projector augmented wave method [20,21] for DFT calculations. A Monkhorst-Pack k-point mesh [22] of 5 × 13 × 2 k-points and a plane wave cut off energy of 400 eV was found to be sufficient for an accuracy of 1 meV/atom. The PBE GGA exchangecorrelation functional [23,24] is used throughout, along with Methfessel-Paxton smearing of 200 meV [25]. Electronic structure is optimised using the Block-Davidson method [26] with a convergence criterion of 0.1 meV.
Our grain boundary was constructed for DFT calculations using the coincident site lattice model [27]  tilt boundaries with a misorientation angle of 129.5 • . The pure Al cell contains 70 atoms (35 per grain) and is volume optimised by fitting ten energy calculations of the cell for a range of lattice parameters to the Birch-Murnaghan equation of state [28]. Three more structures, each containing a cosegregant of Cu, Mg, or Zn in its most energetically favourable substitutional site on the GB, are further optimised along the GB normal axis through fitting the energies of 6 different cell lengths to a quadratic function. Figure 2 shows an example structure containing a Zn atom as a cosegregant. Ionic positions are optimised through a conjugate gradient algorithm [29,30] where ions are moved with respect to calculations of forces at each step, this has a convergence condition of 1meV per atom. N H = 0 to 10 hydrogen atoms are included in multiple configurations interstitial voids on the GB plane, leading to a total of 44 closed GB structures for DFT simulations. The arrangement of H atoms was determined through a systematic process of finding the most energetically favourable site for each additional H atom, excluding the option of H occupying a void already containing H. Given that there are 12 atoms on the GB in Hfree structures, the GB H concentration can be defined as

Ab-initio molecular dynamics
Most ab-initio tensile tests involve a predetermined failure mode, i.e. simulations have an artificially defined fracture path of bonds to break as grains are displaced [31]. For clean Al GBs the correct fracture path is easily identified as the GB plane, but with a high concentration of impurities the failure mode is unclear. To overcome this issue, we developed an ab-initio molecular dynamics tensile test. First, a small amount of vacuum (0.025 nm) is inserted on the GB plane. Second, an electronic structure calculation is peformed and ionic forces are calculated. Ions are moved by 0.025 nm in the directions of force, and ion positions are then optimised using the conjugate gradient algorithm. Finally, another slab of vacuum (0.025nm) is added and the process continues for 25 displacements u i corresponding to each vacuum slab added. To speed up calculations, the central region of the cell (containing the untested GB) is not included in the ionic position optimisation, see Figure 2. The decohesion energy is calculated using χ(u i ) = E(u i ) − E(0), where E(u i ) is the energy of a cell at displacement u i . Traction is obtained using central differencing of the decohesion energies with respect to displacement. It can be hard to decide when interactions across the displaced GB are negligible so a convention of u f = 2W/τ * is adopted to define the decohesion separation when the surfaces are deemed to be no longer interacting [32].
A total of 44 dynamic ab-initio tensile tests are performed, each with fixed Γ H . However small amounts of traction τ < τ * may occur at appreciable time scales to allow for changes in Γ H [33]. To understand the potential for small displacements to drive H transport to a GB, differences in solution enthalpy for H atoms at a clean unbroken GB and a GB separated by u containing Γ H are also calculated as Δμ H (u, Γ H ). Figure 3 shows the work of separation for all 44 tensile tests. For Γ H = 0% the Σ11 Al grain boundary has W = 15.0 eV/nm 2 , which is significantly higher than the result of 12.53 eV/nm 2 from [31]. This discrepancy is explained by [31] not accounting for ligamentation across the boundary through a dynamic tensile test procedure. Results in the literature often employ the rigid grain shift (RGS) method which takes a GB cell and peforms a DFT energy calculation at a range of separations. The RGS method models the interaction of two surfaces from their initial grain boundary equilibrium state to their non-interacting free surface state where χ(u) fits a universal binding energy curve. This method simulates an extremely high strain rate, and most research includes a surface relaxation at each of increment separation -however, these surface relaxations do not substitute our dynamic modelling used in this study. The decohesion energy χ(u) from the dynamic model can be seen in Figure 4 with Mg as a cosegregant. Energy curves do not follow the smooth path found in RGS simulations, which can be explained by discrete structural changes such as the formation (and dissolution) of ligament structures, and by sudden changes in the arrangement of surface H atoms.

Decohesion energy
The results of Figure 3 closely match those found in [34] which looked at a range of Γ H concentrations for Σ3,5, and 9 GBs and included higher Γ H than those explored in this work. In [34] it was reported that in some GBs, a high Γ H can lead to a negative W, although the tensile testing simulation was not performed and instead a simple initial    closed GB and completely separated GB calculation was performed to derive W. This simpler approach hides the possible effects of ligamentation on the final free surface structure and doesn't obtain traction-separation laws for a CZ model of fracture.
Each cosegregant was found to have a small embrittling effect. A 2011 study on Σ5 GBs [35] calculated the best possible substitutional site to place a Zn ion on the boundary in order to model the effects of Zn segregation on boundary cohesion -this study reported that with the Zn atom at different sites, the embrittlement could change from positive to negative; this result was attributed to changes in the lengths of bonds that naturally arise from the substitution of Zn atoms. Charge density difference plots from [35] demonstrate charge depletion which corresponds to preferred fracture paths where weaker bonds will break quickly. Also calculated in [35] is the electron density of states which show that Zn-Al bonds are metallic in character. However, results in [35] show that Zn has a lower bonding charge density than Al, therefore they are weaker than the Al-Al bond which is replaced after substitution. This supports the finding that Zn is an embrittler. On the other hand, in [35] they found that the reduction in W caused by Zn was very small, only 3%. With our dynamic tensile test of Σ11, Zn embrittlement is found to be significantly more potent with a reduction in W of 20%. Another more recent study [36] found strong segregation tendencies for Mg and Cu to the Σ5 GB, however Cu in fact preferentially occupies interstitial sites rather than substitutional. Mg acts as an embrittler and shows similiar charge depletion characteristics to Zn calculated in [35]. In contrast, the interstitial Cu segregants strengthen cohesion across the interface by making new Al-Cu bonds. Whilst the Cu atoms disrupt the mechanical integrity of the GB structure, their chemical bonding effect dominates leading to overall resistance from interfacial separation; in our work, Cu was limited to substitutional sites only, which negates this bonding effect.
The key observation from our atomistic tensile testing is the effect H has on the crystal structure around the GB under strain. For systems containing no H, separation of grains results in significant structural changes around the GB. Effectively, GB atoms transform into ligament structures to form bridges across the vacuum, see Figure 5. The presence of H quite clearly disrupts the formation of such structures, and crystals maintain their face-centered cubic structure throughout a test. This property is common to all cosegregants and provides a possible explanation for changes in W for different Γ H . There is a stark drop in W for increasing H concentrations. We know from RGS tests in the literature that W is reduced when ligaments are artificially prevented from forming, so it is presumed that the strength of low Γ H GBs is derived from these ligaments. There are no significant differences with respect to cosegregants, other than the alloying embrittlement already discussed at Γ H = 0.

Traction-separation laws
The corresponding τ(u) curves were obtained by taking the derivative of χ(u) using central differences with respect to displacement. Due to the discontinuous decohesion energy curves, the form of traction-separation calculated from finite differences is not a smooth curve such as those found using RGS methods. However, the properties of interest for traction-separation laws, maximum (or critical) traction τ * and the critical displacement u * and final displacement u f , can be extracted. Figure 6 shows the maximum traction τ * from each tensile test. A small embrittling effect from the cosegregants is observed at Γ H = 0 where there is a small decrease in τ * for Zn and Cu, and a slightly larger decrease for Mg. There is again a marked reduction in τ * with increasing H content. From Γ H = 0 to 45% there is a reduction in τ * of over 50%. A simple 2nd order polynomial fit of the data roughly follows the data's trend in the region of Γ H studied, where Cu cosegration shows improved strength over pure Al and Zn cosegregation, but Mg cosegregation tends to lead to weaker boundaries. At Γ H = 0, this result is in agreement with [36], and shows that this property is also true for Γ H > 0.
Combining the results for W and τ * allows us to estimate the decohesion separations u f , Figure 7. Whilst there is a lot of noise in the data, combining the quadratic fit from Figures 3 & 6 shows how the decohesion separation will likely change with Γ H . Furthermore Figure 8  shows the grain separation at τ * (= τ(u * )), also known as the critical displacement u * . Ignoring the three anamolous results showing values of approximately 0.4 nm, another decrease is observed with increasing Γ H . These distances play in important part in continuum mechanical modelling of fracture that use traction separation laws as input data. Additionally, in scenarios with variable Γ H , a strained GB could reach critical displacement without increasing traction by decohesion induced H segregation.

H segregation
During the tensile tests, traction gradually increases to τ * at u * and then decreases to u f . Before u * if Γ H were to increase, less traction will be required to displace the grains further. The enthalpy of solution for H was found to not be influenced by Γ H , i.e. the energy to dissolve a H atom into the GB is not affected by how much H is already present in the boundary for the range of concentrations used in this study. [34] found a similiar effect, although they also report that for high Γ H the segregation energy, and enthalpy of solution at the GB, will decrease due to an expansion of the GB caused by Γ H . In this study, these expansions were initially neglected and subsequently accounted for in the dynamic tensile test by shifting χ(u) so that its minimum is centered on u = 0. Figure 9 shows the change in enthalpy of solution for H from a closed GB with no segregants as a fuction of displacement, shown are mean values for each possible cosegregant averaged over all considered Γ H . Mg and Cu cosegregants act to increase hydrogen segregation by a small amount, but separation has a significant impact. For small displacements, due to tractions less than τ * , there is a strong trend in reducing the enthalpy of solution for H, which would drive further segregation, meaning less force is required for continued decohesion. To put these results in perspective, H has been shown using DFT calculations [37] to preferentially occupy vacancies over interstitial voids in bulk face-centered cubic Al with a difference of Δμ H = -0.4eV per vacancy.
By combining these results with models of H transport and crack growth rates, an assessment can be made as to the validity of the HEDE process in EAC. Experimental data obtained from desorption techniques in [38] suggests that there is no tendency for H segregation to GBs (with u = 0) that aren't associated with cracking. For displacements u < u * a question arises as to the effects of strain rate on both Γ H and the formation of ligament structures. In the fast fracture limit, the CZ will break according to RGS simulations with a constant Γ H and conversely, in the slow fracture limit Γ H will increase in accordance with Figure  9. Since DFT simulations suggest lattice stiffenening of elastic moduli from H in solid solution in bulk Al FCC structures, H undergoing transport to the GB but not quite on the GB plane could also enhance decohesion.

Conclusion
We have presented results which simulate decohesion of a CZ model of fracture in Al alloys, using a broad range of segregated H concentrations on a Σ11 symmetric tilt GB. The methods we used capture the complex phenomenon of MATEC Web of Conferences 326, 04006 (2020) ICAA17 https://doi.org/10.1051/matecconf/202032604006 ligament structures forming between grains during microscopic loading of the CZ by using a dynamic ab-initio tensile test. We found that the strength of the GB is derived from the formation of ligaments, and that dissolved H at the GB can disrupt their formation leading to H embrittlement. The work of separation W was found to decrease quickly when the H concentration at the GB, Γ H , increases -for a clean Al GB, W was calculated to be 15 eV/nm 2 but the same boudary containing Γ H =45% has a much lower W of only 4.2 eV/nm 2 . The effects of H on W was consistent with cosegregation of alloying elements Zn, Mg or Cu.
Traction-separation laws were derived for 10 concentrations of H and we found the same trend in the maximum traction τ * as for W. Much weaker correlations occur for critical and final displacements, although a definite trend exists with u * and u f both showing some decrease with respect to increasing Γ H . We also found a strong correlation between displacements of grains and H segregation, where the average enthalpy of solution for a H atom was found to drop by as much as 1eV with a displacement of 0.3nm; an extant Γ H did not seem to affectt the enthalpy of solution of additional H atoms.
There is one small noteable difference in segregation behaviour for a pure GB and each of the cosegregated GBs, namely that the presence of Mg lowers the enthalpy of solution of H -other than this, cosegregation does not influence HEDE. Our results encourage further work to demonstrate the practicality of these theoretical discoveries by use of a multiscale model to include H transport rates and crack growth rates. Our results show how small displacements of grains could potentially drive further H segregation and induce embrittlement, so a model predicting the rate of H transport required to significantly change Γ H during experimental crack growth timeframes is essential to confirm the validity of a HEDE mechanism. On top of this, H production and ingress into the material, along with the presence of GB precipitates, must also be considered as part of the overall HEDE process in EAC of 7xxx alloys.