Computational modelling of hydrogen assisted fracture in polycrystalline materials

A. Valverde-González E. Martínez-Pañeda e.martinez-paneda@imperial.ac.uk A. Quintanas-Corominas J. Reinoso M. Paggi IMT School for Advanced Studies, Piazza San Francesco 19, Lucca 55100, Italy Grupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain Department of Civil and Environmental Engineering, Imperial College London, London SW7 2AZ, UK AMADE, University of Girona, Polytechnic School, c/Universitat de Girona 4, 17003 Girona, Spain

Abstract

We present a combined phase field and cohesive zone formulation for hydrogen embrittlement that resolves the polycrystalline microstructure of metals. Unlike previous studies, our deformation-diffusion-fracture modelling framework accounts for hydrogen-microstructure interactions and explicitly captures the interplay between bulk (transgranular) fracture and intergranular fracture, with the latter being facilitated by hydrogen through mechanisms such as grain boundary decohesion. We demonstrate the potential of the theoretical and computational formulation presented by simulating inter- and trans-granular cracking in relevant case studies. Firstly, verification calculations are conducted to show how the framework predicts the expected qualitative trends. Secondly, the model is used to simulate recent experiments on pure Ni and a Ni-Cu superalloy that have attracted particular interest. We show that the model is able to provide a good quantitative agreement with testing data and yields a mechanistic rationale for the experimental observations.

keywords:

Phase field , Hydrogen embrittlement , Cohesive Zone Model , Elasto-plastic fracture , Finite Element Method

^†^†journal: International Journal of Hydrogen Energy

1 Introduction

Hydrogen is at the core of the most promising solutions to our energy crisis. Hydrogen isotopes fuel the nuclear fusion reaction, the most efficient potentially useable energy process. Moreover, hydrogen is widely seen as the energy carrier of the future and the most versatile means of energy storage. It can be produced via electrolysis from renewable sources, such as wind or solar power, and stored to be used as a fuel or as a raw material in the chemical industry. Hampering these opportunities, hydrogen is known for causing catastrophic failures in metallic structures. Through a phenomenon often termed hydrogen embrittlement, metals exposed to hydrogen containing environments experience a significant reduction in ductility, fracture toughness and fatigue resistance [1, 2, 3]. In the presence of hydrogen, otherwise ductile metals fail in a brittle manner, with cracking often nucleating and propagating along grain boundaries [4, 5]. This ductile-to-brittle shift of metallic alloys in hydrogeneous environments is arguably one of the biggest threats to the deployment of a hydrogen energy infrastructure.

There is a vast literature devoted to shed light into the physical mechanisms behind hydrogen embrittlement [6, 7, 8, 9, 10, 11], and to develop mechanistic predictive models that can prevent failures and map safe regimes of operation [12, 13, 14, 15, 16]. The vast majority of the computational models developed for predicting hydrogen assisted cracking fall into two categories: (i) discrete methods, such as cohesive zone models [17, 18, 19, 20], and (ii) diffuse approaches, such as phase field or other non-local damage models [21, 22, 23, 24, 25, 26]. Cohesive zone models and other discrete methodologies are suitable to describe the nucleation and propagation of sharp cracks through a predefined path. Phase field fracture methods have additional modelling capabilities and can also deliver predictions when the crack trajectory is unknown, when failure is triggered by defects of arbitrary shape, and when the fracture process is complex (e.g., involving the interaction between multiple defects). Both kinds of models can be readily coupled with the hydrogen transport equation and have been successful in qualitatively capturing the main experimental trends, such as the sensitivity to loading rate, hydrogen concentration, and material strength. However, these modelling studies treat materials as isotropic continuum solids, without resolving the underlying microstructure. A number of icrostructurally-sensitive works have been recently carried out (see, e.g., [27, 28, 29, 30, 31, 32, 33] and Refs. therein), but these are limited to capture the interplay between diffusion and deformation, and do not explicitly simulate fracture. The micromechanical fracture of polycrystalline materials in hydrogen-containing environments has been simulated in the works of Rimoli and Ortiz [34], Benedetti et al. [35], and De Francisco et al. [36]. In these works, a cohesive zone formulation was used to predict the failure of grain boundaries, neglecting transgranular cracks.

In this work, we present a new microstructurally-sensitive computational framework for predicting hydrogen assisted fractures. For the first time, the model combines a phase field description of bulk fracture with a cohesive zone model for intergranular cracking. This enables capturing both ductile transgranular fracture and brittle intergranular fracture, and the transition from one to the other. The mechanical and hydrogen transport problems are strongly coupled, with the hydrostastic stress driving hydrogen transport and the hydrogen content reducing the grain boundary strength. The fracture of polycrystalline solids is simulated, with the bulk deformation response being characterised by von Mises plasticity theory. Numerical experiments are conducted to gain insight into the mechanisms of hydrogen-assisted grain boundary decohesion. Focus is on Ni and Ni superalloys, where hydrogen assisted failures are known to be governed by grain boundary decohesion [4, 37]. Among other case studies, the model is used to provide a mechanistic rationale to two recent sets of experiments on Monel K-500 [38] and pure Ni [4] that have attracted particular interest in the hydrogen embrittlement community.

2 Modelling framework

2.1 Fundamentals

Our model deals with an arbitrary body $\Omega\in\mathbb{R}^{n_{dim}}$ with a delimiting external surface $\partial\Omega\in\mathbb{R}^{n_{dim}-1}$ of outward normal n. The body $\Omega$ contains a discrete internal discontinuity $\Gamma$ associated with fracture events in the bulk, and also pre-existing interfaces, arranged in the set $\Gamma_{i}$ . Through exposure to a hydrogen containing environment, hydrogen ingress takes place and thus hydrogen atoms can diffuse through the body, driven by gradients of chemical potential. As detailed in the forthcoming subsections, this leads to a three-field boundary value problem, where the displacement field, the fracture status, and the hydrogen concentration are the primal kinematic variables.

We build our framework upon the assumption of small strains. The vector $\mathbf{x}$ is used to denote the position of an arbitrary point in the global Cartesian system. From a mechanical perspective, the delimiting surface of the body is decomposed into two regions; one where the displacements u are prescribed through Dirichlet-type boundary conditions (BCs), $\partial\Omega_{u}$ , and one where tractions t are prescribed via Neumann-type BC, such that $\partial\Omega=\partial\Omega_{u}\cup\partial\Omega_{t}$ and $\partial\Omega_{u}\cap\partial\Omega_{t}=\emptyset$ . The deformation process is characterized by the small deformation tensor $\bm{\varepsilon}(\mathbf{x})$ , which is defined as the symmetric part of the displacement gradient: $\bm{\varepsilon}(\mathbf{x}):=\nabla^{s}\mathbf{u}(\mathbf{x})$ . Bulk fracture phenomena is here captured using the phase field fracture method [39, 40]. Hence, crack discontinuities are regularised within a diffuse region, whose thickness is characterised by a phase field length scale $\ell$ , and the evolution of the crack-solid interface is described by an auxiliary phase field variable $\phi$ . Regarding hydrogen transport, the external body surface is divided into two parts: the region $\partial\Omega_{q}$ , where the hydrogen flux $\mathbf{J}$ can be prescribed through a Neumann-type BC, and the region $\partial\Omega_{C}$ where the hydrogen concentration $C$ can be prescribed using a Dirichlet BC.

Refer to caption — Figure 1: Sketch summarising the modelling framework: the polycrystalline microstructure of metals exposed to hydrogen is explicitly simulated, with the phase field method being employed to describe the growth of ductile transgranular cracks while hydrogen-assisted grain boundary decohesion is captured by means of a cohesive zone model. Some of the key variables of the model are shown; namely, the phase field order parameter $\phi$ , the phase field length scale $\ell$ (governing the size of the interface region), the hydrogen coverage $\theta_{H}$ , and the cohesive tractions ( $\sigma,\tau$ ) and separations ( $\delta_{n},\delta_{t}$ ).

In addition to bulk fracture, our microstructurally-sensitive formulation employs a cohesive zone model to describe the failure of grain boundaries. As shown in Fig. 1, we explicitly model the polycrystalline microstructure of metals and predict intergranular and transgranular cracking events. The phase field fracture method is used to describe transgranular cracks, which are associated with a ductile fracture process, while a cohesive zone model is employed to predict the decohesion of grain boundaries. Since the latter correspond to a brittle failure and are triggered by the presence of hydrogen, the interfacial fracture energy $\gamma_{c}$ is defined as a function of the hydrogen concentration $C$ . Conversely, the bulk material toughness $G_{c}$ is considered to be a hydrogen-insensitive constant that characterises the resistance of the material to undergo microvoid cracking. Accordingly, for a solid with strain energy density $\psi(\bm{\varepsilon},\phi)$ , the internal functional of the mechanical system comprises the bulk ( $\Pi^{(b)}_{int}$ ) and interfacial ( $\Pi^{(i)}_{int}$ ) contributions, reading

\Pi_{int}(\mathbf{u},\phi,C)=\underbrace{\int_{\Omega/\Gamma_{l}}\psi(\bm{\varepsilon},\phi)\,\mathrm{d}V+\int_{\Gamma_{l}}G_{c}\,\mathrm{d}\Gamma}_{\Pi^{(b)}_{int}}+\underbrace{\int_{\Gamma}\gamma_{c}(C)\,\mathrm{d}\Gamma_{i}}_{\Pi^{(i)}_{int}}

(1)

In the subsequent subsections we proceed to describe each of the features and the constitutive choices of our micromechanical model.

2.2 Chemo-elastoplasticity

The deformation and diffusion problems are intrinsically coupled as hydrogen transport within the crystal lattice is driven by gradients of concentration and hydrostatic stress. We focus our attention on the transport of diffusible hydrogen and consider the influence of traps in the cracking process by using the Langmuir-Mclean isotherm to estimate the hydrogen coverage at grain boundaries. The role of microstructural traps upon slowing down diffusion can also be taken into account through an appropriate choice of the (apparent) diffusion coefficient $D$ . Mass conservation requirements relate the rate of change of the hydrogen concentration $C$ to the hydrogen flux through the external surface,

\int_{\Omega}\frac{\text{d}C}{\text{d}t}\,\mathrm{d}V+\int_{\partial\Omega}\mathbf{J}\cdot\mathbf{n}\,\mathrm{d}S=0

(2)

The strong form of the balance equation can be readily obtained by making use of the divergence theorem and noting that the expression must hold for any arbitrary volume,

\frac{\text{d}C}{\text{d}t}+\nabla\cdot\mathbf{J}=0

(3)

For an arbitrary, suitably continuous, scalar field, $\delta C$ , the variational statement (3) reads,

\int_{\Omega}\delta C\left(\frac{\text{d}C}{\text{d}t}+\nabla\cdot\mathbf{J}\right)\,\mathrm{d}V=0

(4)

Rearranging, and making use of the divergence theorem, the weak form renders,

\int_{\Omega}\left[\delta C\left(\frac{\text{d}C}{\text{d}t}\right)-\mathbf{J}\cdot\nabla\delta C\right]\,\mathrm{d}V+\int_{\partial\Omega_{q}}\delta Cq\,\mathrm{d}S=0

(5)

where $q=\mathbf{J}\cdot\mathbf{n}$ is the concentration flux exiting the body across $\partial\Omega_{q}$ . The diffusion is driven by the gradient of the chemical potential $\nabla\mu$ , with the chemical potential of hydrogen in lattice sites being given by,

\mu=\mu^{0}+RT\ln\frac{\theta_{L}}{1-\theta_{L}}-\bar{V}_{H}\sigma_{H}

(6)

Here, $\theta_{L}$ is the occupancy of lattice sites, which is related to the concentration and number of sites as $\theta_{L}=C/N$ . Also, $\mu^{0}$ denotes the chemical potential in the standard case, and $\sigma_{H}$ is the hydrostatic stress. The last term of (6) corresponds to the so-called stress-dependent part of the chemical potential $\mu_{\sigma}$ , with $\bar{V}_{H}$ being the partial molar volume of hydrogen in solid solution. The mass flux follows a linear Onsager relationship with $\mu_{\sigma}$ , which is often derived from (6) by adopting the assumptions of low occupancy ( $\theta_{L}\ll 1$ ) and constant interstitial sites concentration ( $\nabla N=0$ ) (see, e.g., [41, 42]), such that

\mathbf{J}=-\frac{DC}{RT}\nabla\mu=-D\nabla C+\frac{D}{RT}C\bar{V}_{H}\nabla\sigma_{H}

(7)

where $T$ is the absolute temperature and $R$ is the ideal gas constant. Accordingly, the hydrogen transport equation becomes

\int_{\Omega}\left[\delta C\left(\frac{1}{D}\frac{dC}{dt}\right)+\nabla\delta C\nabla C-\nabla\delta C\left(\frac{\bar{V}_{H}C}{RT}\nabla\sigma_{H}\right)\right]\,\mathrm{d}V=-\frac{1}{D}\int_{\partial\Omega_{q}}\delta Cq\,\mathrm{d}S

(8)

As evident from (6)-(8), the presence of hydrostatic stresses (or volumetric strains) brings a reduction in chemical potential and an increase in hydrogen solubility. Thus, an appropriate description of the stress state in the solid is needed to quantitatively estimate the hydrogen distribution. Here, we choose to describe the deformation of the solid using conventional von Mises plasticity. Accordingly, the total strain energy density of the solid is given by the sum of the elastic and plastic parts,

\psi=\psi^{e}\left(\bm{\varepsilon}^{e}\right)+\psi^{p}\left(\bm{\varepsilon}^{p}\right)=\frac{1}{2}\lambda\left[\text{tr}\left(\bm{\varepsilon}^{e}\right)\right]^{2}+\mu\,\text{tr}\left[\left(\bm{\varepsilon}^{e}\right)^{2}\right]+\int_{0}^{t}\left(\bm{\sigma}:\dot{\bm{\varepsilon}}^{p}\right)\text{d}t\,.

(9)

where $\lambda$ is the first Lame parameter and $\mu$ is the shear modulus. Also, the Cauchy stress tensor is given by $\bm{\sigma}$ , and $\bm{\varepsilon}^{e}$ and $\bm{\varepsilon}^{p}$ respectively denote the elastic and plastic strain tensors. Isotropic power law hardening behaviour is assumed by adopting the following hardening law:

\sigma_{f}=\sigma_{y}\left(1+\frac{E\varepsilon^{p}}{\sigma_{y}}\right)^{(1/n)}

(10)

where $\sigma_{f}$ and $\sigma_{y}$ are the current and initial yield stresses, $E$ is Young’s modulus, $\varepsilon^{p}$ is the accumulated equivalent plastic strain and $n$ is the strain hardening exponent. One should note that, for simplicity, we have chosen to neglect the role of plastic strain gradients; however, large plastic strain gradients arise in the vicinity of cracks and other stress concentrators and lead to large crack tip tensile stresses and hydrogen concentrations [43, 44]. The extension of the present framework to account for the role of plastic strain gradients and geometrically necessary dislocations (GNDs) will be addressed in future work.

2.3 A phase field description of transgranular fractures

The phase field fracture method is used to regularise the internal discontinuity $\Gamma_{\ell}$ , representing the nucleation and growth of transgranular (ductile) cracks. An auxiliary phase field variable $\phi(\mathbf{x})$ is used to describe the evolution of the solid-crack interface, taking the value of $\phi=0$ for the pristine state and of $\phi=1$ for the fully damaged state. Phase field fracture methods have gained remarkable popularity in recent years and have been successfully used to predict cracking in a wide range of materials and applications, including functionally graded solids [45, 46], shape memory alloys [47, 48], and fibre-reinforced composites [49, 50]. The evolution of the phase field equation is dictated by the energy balance associated with the thermodynamics of fracture, as first presented by Griffith [51] and later extended to elastic-plastic solids by Orowan [52]. Thus, under prescribed displacements, the variation of the total energy $\Pi$ due to an incremental increase of the crack area d $A$ is given by,

\frac{\text{d}\Pi}{\text{d}A}=\frac{\text{d}\Psi}{\text{d}A}+\frac{\text{d}W_{c}}{\text{d}A}=0

(11)

where $W_{c}$ is the work required to create two new surfaces and $\Psi=\int\psi\,\text{d}V$ is the stored strain energy. The last term in Eq. (11) denotes the material toughness $G_{c}=\text{d}W_{c}/\text{d}A$ , which can be as low as some tens of J/m² for brittle solids or as high as thousands of kJ/m² for ductile metals where plastic dissipation enhances fracture resistance. The energy balance of Eq. (11) can be cast in a variational form as [53]:

\Pi=\int_{\Omega}\psi\left(\bm{\varepsilon}\right)\,\text{d}V+\int_{\Gamma}G_{c}\,\text{d}\Gamma

(12)

The energy balance is now global and fracture phenomena can be predicted by minimizing the total energy $\Pi$ . However, minimization of Eq. (12) is hindered by the unknown nature of the discontinuous crack surface $\Gamma$ . The phase field regularization can then be exploited to smear this sharp interface into a diffuse region, whose thickness is governed by a phase field length scale $\ell$ . Accordingly, the energy balance (12) can be approximated as [39]:

\Pi_{\ell}=\int_{\Omega}g\left(\phi\right)\psi_{0}\left(\bm{\varepsilon}\right)\,\text{d}V+\int_{V}G_{c}\gamma\left(\phi,\nabla\phi,\ell\right)\,\text{d}V

(13)

where $\gamma(\phi,\nabla\phi)$ is the so-called crack density functional and $g(\phi)$ is the degradation function. We choose to adopt the standard or AT2 phase field formulation [39], and accordingly make the following constitutive choices,

g(\phi)=(1-\phi)^{2}+\kappa\,,

(14)

\gamma(\phi,\nabla,\ell\phi)=\frac{1}{2\ell}\phi^{2}+\frac{\ell}{2}\lvert\nabla\phi\rvert^{2},

(15)

where $\kappa$ is a small numerical parameter to retain residual stiffness when $\phi=1$ , so as to avoid an ill-conditioned system of equations. Noting that $\bm{\sigma}=\partial_{\bm{\varepsilon}}\psi$ , the strong form of the coupled deformation-fracture problem can be readily obtained by inserting (14)-(15) into (13), taking the variation with respect to $\mathbf{u}$ and $\phi$ , and applying Gauss theorem. However, such a formulation would predict cracking also under compressive stress states. Hence, we adopt the so-called volumetric-deviatoric split [54] to decompose the elastic strain energy density into a tensile part,

\psi_{+}^{e}(\bm{\varepsilon}^{e})=\frac{K}{2}\langle\text{tr}[\bm{\varepsilon}^{e}]\rangle_{+}^{2}+\mu[(\bm{\varepsilon}^{e})^{\prime}:(\bm{\varepsilon}^{e})^{\prime}]

(16)

and a compressive part,

\psi_{-}^{e}(\bm{\varepsilon}^{e})=\frac{K}{2}\langle\text{tr}[\bm{\varepsilon}^{e}]\rangle_{-}^{2}

(17)

Here, $K$ is the bulk modulus, $\text{tr}[\bullet]$ is the trace operator, $\langle\bullet\rangle=(\bullet\pm|\bullet|)/2$ and $(\bm{\varepsilon}^{e})^{\prime}=\bm{\varepsilon}^{e}-\text{tr}[\bm{\varepsilon}^{e}]\mathbf{I}/3$ . Furthermore, a history field $\mathcal{H}$ is defined to enforce damage irreversibility. Among the various choices available (see, e.g., [55, 56]), we choose to assume that fracture is driven by the energy stored in the system, in consistency with the energy balance above [40]; accordingly: $\mathcal{H}=\text{max}_{\text{t}\in[0,\tau]}\,\psi_{+}^{e}(\bm{\varepsilon}^{e},\text{t})$ . The local governing equations then read,

\bm{\nabla}\cdot\left\{[(1-\phi)^{2}+\kappa]\bm{\sigma}_{0}\right\}=\mathbf{0}\quad\text{in}\ \Omega

(18)

{G}_{c}\bigg{(}\frac{1}{\ell}\phi-\ell\nabla^{2}\phi\bigg{)}-2(1-\phi)\mathcal{H}=0\quad\text{in}\ \Omega

(19)

where $\bm{\sigma}_{0}$ is the undamaged Cauchy stress tensor. As can be inferred by inspecting (18)-(19), a so-called hybrid approach is used, where the strain energy split is considered only in the phase field evolution equation [57].

2.4 Hydrogen-sensitive interface formulation for grain boundaries

In this modelling framework, the polycrystalline nature of the material is resolved and the decohesion of grain boundaries is explicitly captured by means of a cohesive zone model. Specifically, a traction-separation law is adopted that assumes a tension cut-off relation [58]. This interface formulation assumes a linear and reversible (elastic) evolution until the critical traction is reached. The damage criterion relates the normal $t_{n}$ and tangential $t_{t}$ tractions with their critical counterparts as follows,

\bigg{(}\frac{t_{n}}{t_{nc}}\bigg{)}^{2}+\bigg{(}\frac{t_{t}}{t_{tc}}\bigg{)}^{2}=1

(20)

Accordingly, a critical normal $\delta_{nc}$ and shear $\delta_{tc}$ separations can be defined, which leads to the following definitions of fracture energy for Mode I and Mode II fractures

\gamma_{IC}=\frac{1}{2}t_{nc}\delta_{nc},\quad\gamma_{IIC}=\frac{1}{2}t_{tc}\delta_{tc}

(21)

As sketched in Fig. 1, the role of hydrogen in weakening the grain boundaries is accounted for by degrading the interface fracture energy. Thus, the focus here is on materials that exhibit hydrogen assisted intergranular fracture, such a Ni and Ni alloys [4]. Atomistic calculations have shown a linear relationship between the surface (or fracture) energy and the hydrogen coverage (see, e.g., [59, 60]). Accordingly, and following Martínez-Pañeda et al. [21], we define the following relationship between the fracture energy and the hydrogen coverage $\theta_{H}$ ,

{\gamma_{c}}(\theta_{H})={\gamma_{C,0}}(1-\chi\theta_{H})

(22)

where $\gamma_{C,0}$ is the fracture energy in the absence of hydrogen and $\chi$ is the so-called hydrogen damage coefficient [21]. The same expression is employed for $\gamma_{IC}$ and $\gamma_{IIC}$ . Finally, the hydrogen coverage is estimated from the hydrogen concentration by means of the Langmuir-McLean isotherm:

\theta_{H}=\frac{C}{C+\exp\bigg(\frac{-\Delta g_{b}^{0}}{RT}\bigg{missing})}

(23)

where $\Delta g_{b}^{0}$ is the Gibbs free energy difference between the interface and the surrounding material, also referred to as the segregation energy. Unless indicated otherwise, a value of $\Delta g_{b}^{0}=\ 30\text{kJ}/\text{mol}$ is here employed, based on the spectrum of experimental data available for the trapping energy at grain boundaries [21, 17].

2.5 Numerical implementation

The theoretical model presented in Sections 2.1 to 2.4 is numerically implemented by means of the finite element method. Specifically, the model is implemented in the commercial finite element package ABAQUS/Standard via a user element (UEL) subroutine. In addition, the Abaqus2Matlab software [61] is used to generate the input files and the MATLAB supplementary codes given in Ref. [62] are used to generate the microstructure. Fig. 2 provides a flowchart of the steps followed in the definition and analysis of the microstructure sensitive boundary value problems investigated. Specifically, the microstructure is generated by using a Voronoi-based tessellation algorithm, programmed in MATLAB [62], and this step is followed by the introduction of the resulting microstructure into the ABAQUS input file using a Python script. The coupled deformation-diffusion-fracture system is solved in a staggered fashion, with every sub-problem being solved by means of a backward Euler solution scheme. Typical calculation times are of a few hours.

3 Results

The potential of the theoretical and computational framework presented in Section 2 shall be demonstrated through representative case studies. First, in Section 3.1, a benchmark example is analysed, whereby cracking is predicted in a single edged notched tension specimen containing 50 grains and exposed to various hydrogen-containing environments. Second, we simulate recent slow strain rate tests (SSRTs) on different Monel K-500 lots in Section 3.2, so as to assess the ability of the model in providing a quantitative agreement with experiments and shedding light into the suitability of SSRTs to measure hydrogen susceptibility. Finally, in Section 3.3, the model is used to simulate, for the first time, the seminal experiments by Harris et al. [4] on pure Ni under cryogenic and ambient temperature conditions.

3.1 Benchmark: fracture of a 50-grain SENT plate

Crack nucleation and growth in a square plate microstructure of 50 grains is investigated. The plate, with dimensions given in Fig. 3.1a, contains a small notch and is subjected to uniaxial loading; a testing configuration often referred to as single edge notched tension (SENT) specimen. The load is applied by prescribing the vertical displacement at the upper edge, at a rate of $\dot{u}_{y}=10^{-10}$ mm/s, while the vertical displacement is constrained at the bottom edge. To prevent rigid body motion, the horizontal displacement is constrained at the bottom-right corner. The sample is exposed to hydrogen on its left side, where the notch is located. No pre-charging time in considered, with both hydrogen and mechanical charging starting simultaneously. The magnitude of the environmental hydrogen concentration $C_{\text{env}}$ is varied between 0 and 0.9 ppm wt with the aim of capturing a reduced critical load with increasing hydrogen content, and a transition from ductile (bulk) fracture to intergranular cracking. The material properties assumed for the bulk and the interface are given in Table 1. The sample is discretised using a total of 48,554 quadratic quadrilateral elements.

Table 1: Material properties of the bulk and the interface adopted in the SENT benchmark case study.

Bulk properties
$E\ (\text{GPa})$	$\nu$	$\ell\ (\text{mm})$	$G_{c}\ (\text{kJ}/\text{m}^{2})$	$D\ (\text{mm}^{2}/\text{s})$	$\sigma_{y}\ \mathbf{(\text{MPa})}$	n
185	0.3	0.025	0.05	$1.3\times 10^{-8}$	794.3	0.064
Interface properties
$k_{n},k_{t}\ (\text{MPa}/\text{mm})$		$t_{nc,0},t_{tc,0}\ (\text{MPa})$		$\gamma_{IC,0},\gamma_{IIC,0}\ (\text{kJ}/\text{m}^{2})$		$\chi$
$2\times 10^{8}$		$2.25\times 10^{3}$		$1.27\times 10^{-2}$		0.86

The results obtained for the case of $C_{\text{env}}=0$ ppm wt are shown in Figs. (3.1b-3.1g), in terms of the phase field contours. In the absence of hydrogen, the crack nucleates inside of the grain, in the vicinity of the notch tip, and propagates in a transgranular manner. The damage appears to accumulate within the grain region closer to its boundary, as a stress mismatch takes place between adjacent grains due to differences between the stiffness of the bulk and that of the interface.

Heat	$E\ (\text{GPa})$	$\nu$	$\sigma_{y}\ (\text{MPa})$	n
Allvac	180	0.3	794.3	0.058
NRL LS	198	0.3	715.7	0.054
TR2	202	0.3	795	0.055
NRL HS	191	0.3	910.1	0.050

Heat	$E_{A}=-0.85\ \text{V}_{\text{SCE}}$	$E_{A}=-0.95\ \text{V}_{\text{SCE}}$	$E_{A}=-1.1\ \text{V}_{\text{SCE}}$
Allvac	1.9	4.1	7.5
NRL LS	1.3	4.7	18.3
TR2	3.7	18.6	26.2
NRL HS	4.7	11.9	23.4

	Allvac	NRL LS	TR2	NRL HS
$k_{n},k_{t}$ ( $10^{8}\ \text{MPa}/\text{mm}$ )	$2.00$	$2.00$	$2.00$	$2.00$
$t_{nc,0},t_{tc,0}$ ( $10^{4}\ \text{MPa}$ )	$2.07$	$2.78$	$2.23$	$2.13$
$\gamma_{IC,0},\gamma_{IIC,0}$ ( $\text{kJ}/\text{m}^{2}$ )	$1.07$	$1.92$	$1.24$	$1.12$
$\chi$	0.85	0.82	0.86	0.79

$\mathbf{E_{A}}\ (\text{V}_{\text{SCE}})$	Allvac	NRL LS	TR2	NRL HS
0 (Air)	TG ✓	TG ✓	TG ✓	TG ✓
0.85	TG ✓	TG $\bm{\times}$	TG ✓	-
0.95	TG ✓	IG ✓	IG ✓	TG ✓
1.1	IG ✓	IG ✓	IG ✓	IG ✓

Temperature	77 K		RT
	no H	H	no H	H
$E\ (\text{GPa)}$	227	227	202	202
$\nu$	0.3	0.3	0.3	0.3
$\sigma_{y}\ (\text{MPa})$	222	233	182	192
$n$	0.159	0.159	0.140	0.140
$D\ (\text{mm}^{2}/\text{s})$	$10^{-15}$	$10^{-15}$	$10^{-9}$	$10^{-9}$