Computational predictions of hydrogen-assisted fatigue crack growth

The phase field description of fracture is based on the thermodynamics of fracture, as first postulated by Griffith [40]. Griffith’s theory describes fracture as a competition between the energy stored in the bulk and the energy required to create two new surfaces. In other words, crack propagation occurs when the energy stored in the material is sufficient to overcome the material toughness $G_{c}$. Denoting $\psi$ to be the density of elastic strain energy stored in the material and T a traction acting over the surface of the body, the total energy functional $\Pi$ of the cracked solid can be written as,

where $\mit\Gamma$ is the discrete crack surface and $\bm{\varepsilon}$ is the strain tensor, which is a function of the displacement field as $\bm{\varepsilon}=\left(\nabla\textbf{u}+\nabla\textbf{u}^{T}\right)/2$. Minimisation of Eq. (1) would enable capturing the evolution of cracks as an exchange between elastic and fracture energies. However, this problem is computationally intractable since the crack surface $\mit\Gamma$ is usually unknown. Phase field modelling enables overcoming this problem by introducing an auxiliary variable, the phase field order parameter $\phi$, to track the crack interface [41, 42]. As sketched in Fig. 2, the order parameter $\phi$ varies smoothly from 0, at intact material points, to 1, in fully cracked regions, smearing the sharp crack into a diffuse region. Subsequently, the spatial and temporal evolution of the crack, including its nucleation, propagation, branching, and merging behaviours, becomes a natural by-product of Griffith’s energy balance. Adopting the so-called conventional or AT2 phase field model, the total potential energy (1) can be approximated by

with $\ell$ being the length scale parameter that governs the size of the smooth fracture process zone. Through $\Gamma$-convergence, it can be rigorously proven that the regularised function Eq. (2) converges to the Griffith’s theory Eq. (1) when $\ell\rightarrow 0^{+}$. This smeared description of fracture has enjoyed great success, with phase field formulations being developed for a number of fracture problems [43, 44, 45, 46], from small-scale cracking in battery materials [47, 48] to crevasse growth in large ice sheets and glaciers [49, 50].

The balance equations can be readily derived by considering appropriate constitutive choices. First, following the usual assumption of a quadratic degradation function, the strain energy density in Eq. (2) can be given as,

where $\psi_{0}(\bm{\varepsilon})$ is the undamaged strain energy density, which for elastic solids is given by

with $\textit{{C}}_{0}$ being the (undamaged) linear elastic stiffness tensor, as defined by appropriate choices of Young’s modulus $E$ and Poisson’s ratio $\nu$. Since the focus is on fatigue crack growth under small-scale yielding conditions, plasticity effects can be neglected. Nevertheless, extending the present framework to incorporate plasticity is relatively straightforward (see Ref. [51]). Consistent with Eq. (4), the stress tensor is given by $\bm{\sigma}=g(\phi)\,\bm{\sigma}_{0}=\left(1-\phi\right)^{2}\partial\psi_{0}(\bm{\varepsilon})/\partial\bm{\varepsilon}$, with $\bm{\sigma}_{0}$ being the undamaged stress tensor. Considering these constitutive choices, one can derive the governing equations by taking the stationary of Eq. (2) with respect to the primary fields, $\mathbf{u}$ and $\phi$, and applying the divergence theorem; such that,

with the following associated boundary conditions,

Finally, it should be noted that the choice of $\ell$ defines the material strength $\sigma_{c}$ [52, 53]. This can be readily observed by considering the analytical solution for Eq. (5b) in the case of a homogeneous 1-D scenario ($\nabla\phi=0$), where $\psi_{0}=E\varepsilon^{2}/2$ and the stress and strain distributions reach a maximum at

Let us now consider the role of hydrogen. To this end, one must consider the diffusion of hydrogen within the metal, which governs the loading frequency effect. Hydrogen transport is governed by mass conservation, such that the rate of change in time ($t$) of the hydrogen concentration ($C$) must be equal to the concentration flux $\mathbf{J}$ through the boundary of the body,

Since mass conservation must hold for any volume, Eq. (8) can be re-written in a point-wise form as

The hydrogen flux $\mathbf{J}$ is given by

where $D$ is the (apparent) hydrogen diffusion coefficient, $R_{g}$ is the gas constant, $T$ is the absolute temperature, $\bar{V}_{\mathrm{H}}$ is the partial molar volume of hydrogen in the solid solution, and $\sigma_{\mathrm{h}}$ is the hydrostatic stress. The last term in Eq. (10) account for the effects of mechanical straining on the increase in hydrogen solubility, a result of lattice distortion.

Inserting Eq. (10) into Eq. (9) renders the governing equation for hydrogen transport,

It remains to define the coupling between the hydrogen diffusion and fracture problems. Since hydrogen reduces the toughness of metals, this is commonly done through the definition of a hydrogen damage function $f_{\mathrm{H}}\left(C\right)$, such that the fracture energy becomes $f_{\mathrm{H}}\left(C\right)G_{c}$. Two approaches are generally taken to define $f_{\mathrm{H}}\left(C\right)$. A mechanistic one, whereby a degradation law is defined based on a multi-scale interpretation of the underlying mechanisms (e.g., quantitatively degrading $G_{c}$ via first-principles calculations [27]). Or a phenomenological one, whereby the toughness sensitivity to hydrogen content is characterised experimentally via fracture toughness tests at various hydrogen contents. Here, the latter is used since, as described below, significant fracture toughness versus hydrogen content data exist for pressure vessel steels. These data are best fitted through an exponential degradation law of the form [51],

where $\xi$, $\eta$, and $b$, the parameters that govern the extent of the fracture toughness degradation due to hydrogen, are to be calibrated with experiments.

The last ingredient of the model is the consideration of fatigue damage, and its couplings. The extension of phase field fracture to fatigue adopted here follows the framework established by Golahmar et al. [37], which is itself based on the work by Carrara et al. [54] and extends it to capture arbitrary S-N curve slopes, the fatigue endurance limit, and the mean stress effect. Accordingly, the material toughness is degraded by a fatigue degradation function $f_{\mathrm{F}}(\bar{\alpha})$, which evolves based on a cumulative history variable $\bar{\alpha}$, which is expressed as a function of a variable $\alpha$ to capture the physics of fatigue damage. Consistent with the fracture energy balance discussed above, fatigue damage is assumed to be driven by the energy stored in the solid and accordingly, the parameter $\alpha$ is defined as $\alpha=\psi$. The fatigue degradation function is given by [54]

where $\bar{\alpha}_{0}$ is a material parameter determining the rate at which hydrogen affects the degradation of fracture toughness, to be calibrated with experiments. The fatigue history variable $\bar{\alpha}$ at a time $t$ is given by

And the definition of $\Delta\bar{\alpha}$ is a key aspect of the model. Here, following Golahmar et al. [37], we use a generalised definition of $\Delta\bar{\alpha}$ that accounts for (i) the exponent term of the Paris law, (ii) the effect of the stress ratio, and (iii) the fatigue threshold. Specifically, $\Delta\bar{\alpha}$ is defined as

where $n$ is a material parameter to match the slope of the S-N curve, $\alpha_{\max}$ is the maximum value of $\alpha$ within a fatigue cycle, $\alpha_{n}=\sigma_{c}\varepsilon_{c}/2$ is a normalisation parameter to achieve dimensional consistency, $R$ is the stress (load) ratio, $\kappa$ is a material constant measuring the sensitivity of fatigue damage to mean stress effect, $H()$ is the Heaviside function and $\alpha_{e}$ is the fatigue threshold corresponding to the material endurance stress. As per Ref. [37], the exponent term $n$ can be related to the exponent $m$ in the Paris law equation, $\mathrm{d}a/\mathrm{d}N=C(\Delta K)^{m}$, as $n=0.49m-0.61$.

Considering both the fatigue and hydrogen degradation functions, the phase field evolution equation (5b) could be reformulated as⁸⁸8While this is the simplest form of the phase field evolution equation, it should be noted that if the toughness degradation is introduced in the weak form, gradient terms for the fatigue and hydrogen degradation functions naturally appear in the balance equation. Here, we follow Golahmar et al. [38] and consider only the former, neglecting the role of $\nabla f_{\mathrm{H}}(C)$.,

Hence, both fatigue damage and hydrogen are assumed to degrade the material toughness synergistically. Importantly, we assume that there is no other coupling between the accumulation of hydrogen and the fatigue process. Fatigue crack growth rates will exhibit sensitivity to hydrogen content only through the hydrogen sensitivity of the material toughness. This reduces the number of parameters at play and is consistent with experimental observations - hydrogen significantly degrades the fracture toughness and accelerates fatigue crack growth rates but has a limited impact on the accumulation of fatigue damage, as evidenced by the relatively small effect on S-N curves. However, a more general version of our model is presented in 3.5, where the influence of hydrogen in accelerating the ability of a material to experience cyclic damage is accounted for. This leads to a slightly better agreement with experiments at the cost of having one additional parameter in the model.

Mimicking the laboratory experiments, simulations are conducted using Compact Tension (CT) specimens with the dimensions given in Fig. 3. Taking advantage of symmetry, only the upper half of the sample is modelled. The finite element mesh, shown in Fig. 3b, employs approximately 5,000 plane strain quadratic quadrilateral elements, with the mesh being particularly fine in the expected crack growth region. To ensure mesh-independent results, the characteristic element size is taken to be at least 6 times smaller than the phase field length scale $\ell$. A cyclic load $P$, of range $\Delta P$, is applied as a nodal force on the pin. Following the ASTM E647 standard [55], the stress intensity factor range $\Delta K$ is then obtained from the applied load using the following expression,

where $a$ is the current crack length, $W$ is the effective width, and $B$ is the sample thickness.

Consistent with the experimental data, we simulate the behaviour of a pressure vessel steel, referred to as ASME SA-723 Grade 1, Class 1. This material has a Young’s modulus of $E=210$ GPa and a Poisson’s ratio of $\nu=0.3$. In the absence of hydrogen, the fracture behaviour of this material is described by a toughness of $G_{c}=100$ $\mathrm{kJ/m^{2}}$ [56], and a critical strength that is four times the yield stress [51], which is $\sigma_{y}=715$ MPa [19]. Considering the relationship outlined in Eq. (7)a, this results in a phase field length scale of $\ell=0.27$ mm. In terms of its fatigue behaviour, the model requires defining the parameters $n$, $\kappa$ and $\bar{\alpha}_{0}$. The sensitivity to the stress ratio is characterised by a choice of $\kappa=0.78$, based on the experiments conducted by Dowling et al. [57] on a similar material. The parameter $n$, governing the slope of the Paris and S-N curves, is chosen to be $n=1.25$, to fit the slope observed in air. Finally, $\bar{\alpha}_{0}$ is taken to be equal to $\bar{\alpha}_{0}=8$ N/mm², so as to best capture the experimental fatigue crack growth curve in the absence of hydrogen. In this regard, we emphasise that the fatigue behaviour of the material in air ($p_{\mathrm{H_{2}}}=0$) is an input to the model (i.e., it requires parameter calibration), unlike the fatigue behaviour in a hydrogen-containing environment, which is obtained without any changes to the model parameters and is, therefore, a prediction. For completeness, we also define the fatigue threshold parameter, which is taken to be $\alpha_{e}=0.05$ (corresponding to an endurance stress of 150 MPa), although this parameter has no effect on the type of crack growth experiments conducted here [37].

Hydrogen is assumed to influence material behaviour only through the degradation of the material toughness $G_{c}$. To characterise the sensitivity of $G_{c}$ to hydrogen content, experimental fracture toughness versus hydrogen gas pressure data are fitted with a so-called hydrogen degradation function $f_{\mathrm{H}}(C)=G_{c}/G_{c}(C=0)=J_{Ic}/J_{Ic}(C=0)$, where $G_{c}(C=0)$ ($J_{Ic}(C=0)$) denotes the toughness in an inert environment. This function is an intrinsic material characteristic. To obtain $f_{\mathrm{H}}(C)$ from toughness versus $p_{\mathrm{H_{2}}}$ data, we first need to determine the hydrogen concentration corresponding to each of the hydrogen pressures considered in the experiments. To this end, Sievert’s law is used [58],

where the solubility of steel is defined as $S=0.077$ wppm $\mathrm{MPa}^{-0.5}$ [59]. Under the sensible assumptions of small scale yielding and plane strain conditions, the toughness $G_{c}$ can be readily related to the critical $J$-integral or fracture toughness $K_{Ic}$ as,

Accordingly, data for various hydrogen pressures can be condensed in a $J_{Ic}/J_{Ic}(0)$ (or $f_{\mathrm{H}}$) versus $C$ plot, as shown in Fig. 4, which uses (static) fracture toughness experimental data for similar pressure vessel steels. The data can be suitably fitted with a degradation function of the type given in Eq. (12), with $\xi=0.12$, $b=2$ and $\eta=7$.

Finally, the transport of hydrogen requires the definition of the partial molar volume, $\bar{V}_{\mathrm{H}}=2000$ mm³/mol for steels, and the effective hydrogen diffusion coefficient, which is taken to be $D=2\times 10^{-4}$ $\mathrm{mm^{2}/s}$, based on the range of experimental data reported [61, 62]. The boundary conditions for the hydrogen transport problem mimic the conditions of the laboratory experiments. That is, starting from a hydrogen-free sample ($C(t=0)=0\,\forall\,\mathbf{x}$), the outer surfaces of the specimen are exposed to the hydrogen concentration associated with the H₂ pressure used in the test, as per Sievert’s law (18). This environmental hydrogen content is referred to as $C_{\text{env}}$. As in the experiments, no mechanical load is applied in the first 24 h of hydrogen gas exposure (the ‘soaking’ time associated with the stabilisation of the testing setup). Once the crack grows in a hydrogen-containing environment, we capture how the newly created crack surfaces are promptly exposed to the hydrogen environment by using a so-called penalty boundary condition [63, 64], which enforces $C=C_{\text{env}}$.

The aforementioned geometry, model, and parameters are used to run Virtual fatigue crack growth experiments, in both air and hydrogen gas environments, and for various choices of loading frequency, load ratio and hydrogen gas pressure, mimicking the conditions of the laboratory tests. The model can simulate every load cycle and therefore naturally accommodate different load ratios and loading frequencies, as in the experiments. A representative example of the qualitative outcome of these computational experiments is given in Fig. 5, as computed for the case of a H₂ pre-charged sample exposed to a loading frequency of $f=1$ Hz, a hydrogen pressure of $p_{\mathrm{H_{2}}}=106$ MPa and a loading ratio of $R=0.1$. Fig. 5a depicts the crack growth process, through the contours of the phase field variable $\phi$. The crack front is defined to be in the middle of the diffuse interface ($\phi=0.5$). Fig. 5b shows hydrogen concentration contours, showing how hydrogen accumulates in areas of high hydrostatic stress, such as the moving crack tip. From the crack extension and the applied load, see Eq. (17), d$a$/d$N$ vs $\Delta K$ curves are obtained.

The ability of the model to predict the behaviour of pressure vessel steels in the combined presence of hydrogen and cyclic loading is examined next. Recall that the model has been endowed with the material’s fatigue characteristics (in air) and with the sensitivity of $G_{c}$ to hydrogen, yet no criteria or parameter has been defined to predict hydrogen-assisted fatigue. As shown in Fig. 6, the predictive abilities of the model are benchmarked against a wide range of $\Delta K$ values and scenarios, spanning different hydrogen pressures, load ratios and loading frequencies. The results exhibit a remarkable degree of agreement across all the conditions examined.

Consider first, Fig. 6a, where the influence of the H₂ pressure is investigated. As stated above, changes in H₂ pressure result in different environmental hydrogen concentrations, as per Sievert’s law (18) and this is introduced into the simulation through the boundary conditions $C=C_{\text{env}}$. The results show that the model, calibrated to fit the response in air, is able to adequately predict the experimental data for both $p_{\mathrm{H_{2}}}=55$ and $p_{\mathrm{H_{2}}}=106$ MPa. As in the experiments, these numerical calculations are obtained with a load ratio of $R=0.1$ and a loading frequency of $f=1$ Hz. Next, the load ratio is varied to assess the ability of the model to predict the interplay between stress ratio and hydrogen embrittlement susceptibility. The results, shown in Fig. 6b, reveal that the model quantitatively captures how fatigue crack growth rates become larger for higher $R$ values in the presence of hydrogen. The agreement is remarkable across the entire range, $R=0.1-0.7$. The numerical and experimental results shown in Fig. 6b were obtained for a loading frequency of $f=1$ Hz and a hydrogen pressure of $p_{\mathrm{H_{2}}}=55$ MPa, but experimental data for various $R$ values is also available for a pressure of $p_{\mathrm{H_{2}}}=106$ MPa. Thus, the ability of the model to predict fatigue crack growth rates across $p_{\mathrm{H_{2}}}$ and $R$ values can also be assessed - this is shown in Fig. 6c, where again a very good agreement with experiments is observed across the entire range of conditions considered. Finally, the role of the loading frequency is examined in Fig. 6d. Both the experiments and the numerical results show little sensitivity to the loading frequency. Overall, average trends in the experimental data suggest that fatigue crack growth rates are higher for lower loading frequencies, although the $f=0.1$ Hz and $f=1$ Hz curves overlap and cross each other for low $\Delta K$ values. The numerical results seem to reinforce this observation, with the $f=0.1$ Hz result being above the $f=1$ Hz one but with little differences being observed. While hydrogen-assisted fatigue crack growths are known to be sensitive to the loading frequency, the longer the loading cycle, the more time for hydrogen to accumulate in the process zone, a small sensitivity to the loading frequency is expected in the two extreme regimes (sufficiently fast or slow tests). Mapping the fatigue behaviour regimes that result from the interplay between material diffusivity and loading frequency is therefore key to optimising testing protocols and ensuring conservatism; this is the objective of the next subsection.

To unravel the interplay between hydrogen transport, loading frequency and fatigue crack growth rates in pressure vessel steels, numerical experiments are conducted over a wide range of loading frequencies, from $f=0.001$ Hz to $f=100$ Hz. There is a lack of experimental data at low frequencies due to the challenges associated with conducting tests at those conditions (large time scales, leaks). However, experimental data for a wide range of loading frequencies is available for SA-372 Gr. J steel, a similar material with $\sigma_{y}=760$ MPa (see Ref. [65]), and this is used to further benchmark model predictions. The results obtained are shown in Fig. 7, for three selected values of the stress intensity factor range: $\Delta K=8$ MPa$\sqrt{\text{m}}$, $\Delta K=15$ MPa$\sqrt{\text{m}}$ and $\Delta K=24$ MPa$\sqrt{\text{m}}$. The load ratio is $R=0.1$ and the H₂ pressure equals $p_{\mathrm{H_{2}}}=103$ MPa. Three observations should be particularly emphasised. Firstly, the model is shown to accurately capture the three expected regimes of behaviour: (i) a “fast test” regime where hydrogen does not have time to diffuse and d$a$/d$N$ values are lowest, (ii) a “slow test” regime, where hydrogen has sufficient time to diffuse to the fracture process zone region and d$a$/d$N$ values are highest, and (iii) an intermediate regime, where fatigue crack growth rates are susceptible to the loading frequency. Among others, the range of loading frequencies that this intermediate regime spans is dependent on the $\Delta K$ value and, mainly, the diffusivity of the material. Secondly, computational and experimental results appear to be in relatively good agreement, despite the significant scatter observed in the latter. Some differences are observed in the “slow test” regime (very low $f$) for the highest load $\Delta K$. These could be potentially related to the challenges intrinsic to low loading frequency tests (e.g., a higher propensity for leaks). Thirdly, model results show that the impact of loading frequency becomes more significant as $\Delta K$ increases. The sensitivity to $\Delta K$ is notable and can be explained by the dependency of the peak hydrostatic stress (and thus peak $C$) on $K$ and by the interplay between hydrogen transport and crack growth rate, with the latter increasing significantly with $\Delta K$. All virtual experiments conducted in the mid-low $\Delta K$ regime show a lower sensitivity to the loading frequency and a wider range of conservative test frequencies. A loading frequency of 1-10 Hz appears to be sufficient if $\Delta K$ is on the order of 8 MPa$\sqrt{\text{m}}$. On the other hand, for values of the stress intensity factor range as high as $\Delta K=24$ MPa$\sqrt{\text{m}}$, our results show that one has to reduce the loading frequency to roughly 0.1 Hz to be sufficiently close to the “slow test” regime and therefore conservative. Tests conducted at $f=1$ Hz (a common habit), appear to under-predict by a factor of 2 the fatigue crack growth representative of engineering applications (where $f\ll 1$ Hz). Nevertheless, the high $\Delta K$ regime is arguably of less technological significance as most pressure vessels for H₂ storage are currently designed to operate at mid to low $\Delta K$ values.

Finally, the model is used to gain insight into the influence of the so-called ‘soak time’, the period when the sample is exposed to hydrogen before mechanical loading starts. In the experiments considered so far, this time was approximately 24 hours, as this is the time that it takes to stabilise the load cell. With state-of-the-art Linear Variable Differential Transformer (LVDT) load cells this can be reduced to minutes. In general, samples are not typically exposed to a hydrogen-containing environment prior to charging to reduce the testing time and to emulate the scenario where the material is protected from hydrogen uptake until a fatigue crack develops (e.g., because of the presence of an oxide or an engineering coating). However, in components experiencing long-term exposure to hydrogen, one can readily envisage a situation where the material has absorbed hydrogen before the onset of fatigue crack growth. Consequently, we use our model to quantify the deviation between the fast, no pre-charging condition and the potentially more technologically relevant scenario of hydrogen saturation before mechanical testing. To this end, calculations are conducted with the following initial conditions: $C(t=0)=C_{\mathrm{env}}\,\forall\,\textbf{x}$. That is, a uniform hydrogen concentration (the one corresponding to the relevant applied pressure) throughout the sample. The results obtained are shown in Fig. 8. The numerical experiments were conducted for a load ratio of $R=0.1$, a loading frequency of $f=1$ Hz and two H₂ pressures: $p_{\mathrm{H_{2}}}=55$ MPa and $p_{\mathrm{H_{2}}}=106$ MPa. The obtained fatigue crack growth curves show a very small effect at low $\Delta K$ values but a larger deviation at high $\Delta K$ values. Also, differences appear to be slightly smaller with increasing $p_{\mathrm{H_{2}}}$, as the toughness sensitivity to changes in hydrogen content is smaller for high $C$ values (see Fig. 4). The qualitative trends observed are not surprising, as in the high $\Delta K$ regime the crack is growing into a region that was free of hydrogen in the conventional test condition. This non-uniform distribution of the material toughness is much less pronounced in the fully charged numerical experiments, leading to higher d$a$/d$N$ values. Less sensitivity to the loading frequency is also expected in the pre-charged condition, although there will still be a dependency due to the interplay between hydrostatic stresses and hydrogen transport.

The model proposed is capable of accurately predicting fatigue crack growth rates based solely on one element of coupling between the damage process and hydrogen: the degradation of fracture toughness with increasing hydrogen content. This assumes that hydrogen only reduces the fracture resistance of the material and has no effect on the ability of a material to develop damage due to cyclic loading. This is arguably a sensible assumption, as the impact of hydrogen in fatigue is mostly through an enhancement of fatigue crack growth rates, and the results obtained reveal that this assumption enables establishing a model capable of delivering accurate predictions with minimal calibration. However, changes in S-N curve behaviour with increasing hydrogen content have been reported [66], which suggests that hydrogen can also accelerate the nucleation of fatigue cracks. While the present model can naturally capture this through the relationship between strength and (hydrogen-degraded) toughness, the extent to which this agrees quantitatively with experiments has not been thoroughly investigated yet. Thus, an argument could be made for the development of a more flexible model, capable of accounting for a hydrogen-fatigue interplay through the addition of an additional parameter. To showcase the potential improvement in accuracy that such a model would bring, fatigue crack growth curves are re-calculated considering a different $n$ value for hydrogen-containing environments. Specifically, Fig. 9 shows the results obtained using $n=1.9$ for $p_{\mathrm{H_{2}}}>0$ (with $n=1$ for $p_{\mathrm{H_{2}}}=0$). A slightly better agreement is attained, relative to the results presented in Fig. 6.