{'fraction_non_alphanumeric': 0.030216647662485746, 'fraction_numerical': 0.010547320410490307, 'mean_word_length': 5.0461509000383, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational modelling of hydrogen assisted fracture in polycrystalline materials

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.


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.

Fundamentals
Our model deals with an arbitrary body Ω ∈ R n dim with a delimiting external surface ∂Ω ∈ R n dim −1 of outward normal n. The body Ω contains a discrete internal discontinuity Γ associated with fracture events in the bulk, and also pre-existing interfaces, arranged in the set Γ 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 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), ∂Ω u , and one where tractions t are prescribed via Neumann-type BC, such that ∂Ω = ∂Ω u ∪ ∂Ω t and ∂Ω u ∩ ∂Ω t = ∅. The deformation process is characterized by the small deformation tensor ε(x), which is defined as the

Phase field
Cohesive zone model 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 φ, the phase field length scale (governing the size of the interface region), the hydrogen coverage θ H , and the cohesive tractions (σ, τ ) and separations (δ n , δ t ).
symmetric part of the displacement gradient: ε(x) := ∇ s u(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 , and the evolution of the crack-solid interface is described by an auxiliary phase field variable φ. Regarding hydrogen transport, the external body surface is divided into two parts: the region ∂Ω q , where the hydrogen flux J can be prescribed through a Neumann-type BC, and the region ∂Ω C where the hydrogen concentration C can be prescribed using a Dirichlet BC.
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 γ 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 ψ(ε, φ), the internal functional of the mechanical system comprises the bulk (Π In the subsequent subsections we proceed to describe each of the features and the constitutive choices of our micromechanical model.

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, 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, For an arbitrary, suitably continuous, scalar field, δC, the variational statement (3) reads, Rearranging, and making use of the divergence theorem, the weak form renders, where q = J · n is the concentration flux exiting the body across ∂Ω q . The diffusion is driven by the gradient of the chemical potential ∇µ, with the chemical potential of hydrogen in lattice sites being given by, Here, θ L is the occupancy of lattice sites, which is related to the concentration and number of sites as θ L = C/N . Also, µ 0 denotes the chemical potential in the standard case, and σ H is the hydrostatic stress. The last term of (6) corresponds to the so-called stress-dependent part of the chemical potential µ σ , withV H being the partial molar volume of hydrogen in solid solution. The mass flux follows a linear Onsager relationship with µ σ , which is often derived from (6) by adopting the assumptions of low occupancy (θ L 1) and constant interstitial sites concentration (∇N = 0) (see, e.g., [41,42]), such that where T is the absolute temperature and R is the ideal gas constant. Accordingly, the hydrogen transport equation becomes 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, where λ is the first Lame parameter and µ is the shear modulus. Also, the Cauchy stress tensor is given by σ, and ε e and ε p respectively denote the elastic and plastic strain tensors. Isotropic power law hardening behaviour is assumed by adopting the following hardening law: where σ f and σ y are the current and initial yield stresses, E is Young's modulus, ε 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.

A phase field description of transgranular fractures
The phase field fracture method is used to regularise the internal discontinuity Γ , representing the nucleation and growth of transgranular (ductile) cracks. An auxiliary phase field variable φ(x) is used to describe the evolution of the solid-crack interface, taking the value of φ = 0 for the pristine state and of φ = 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 Π due to an incremental increase of the crack area dA is given by, where W c is the work required to create two new surfaces and Ψ = ψ dV is the stored strain energy. The last term in Eq. (11) denotes the material toughness G c = dW c /dA, which can be as low as some tens of J/m 2 for brittle solids or as high as thousands of kJ/m 2 for ductile metals where plastic dissipation enhances fracture resistance. The energy balance of Eq. (11) can be cast in a variational form as [53]: The energy balance is now global and fracture phenomena can be predicted by minimizing the total energy Π. However, minimization of Eq. (12) is hindered by the unknown nature of the discontinuous crack surface Γ. 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 .
Accordingly, the energy balance (12) can be approximated as [39]: where γ(φ, ∇φ) is the so-called crack density functional and g(φ) is the degradation function.
We choose to adopt the standard or AT2 phase field formulation [39], and accordingly make the following constitutive choices, where κ is a small numerical parameter to retain residual stiffness when φ = 1, so as to avoid an ill-conditioned system of equations. Noting that σ = ∂ ε ψ, 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 u and φ, and applying Gauss theorem. However, such a formulation would predict cracking also under compressive stress states. Hence, we adopt the so-called volumetricdeviatoric split [54] to decompose the elastic strain energy density into a tensile part, and a compressive part, Here, K is the bulk modulus, tr[•] is the trace operator, • = (• ± | • |)/2 and (ε e ) = ε e − tr[ε e ]I/3. Furthermore, a history field 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: . The local governing equations then read, where σ 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].

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, Accordingly, a critical normal δ nc and shear δ tc separations can be defined, which leads to the following definitions of fracture energy for Mode I and Mode II fractures 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 θ H , where γ C,0 is the fracture energy in the absence of hydrogen and χ is the so-called hydrogen damage coefficient [21]. The same expression is employed for γ IC and γ IIC . Finally, the hydrogen coverage is estimated from the hydrogen concentration by means of the Langmuir-McLean isotherm: where ∆g 0 b 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 ∆g 0 b = 30kJ/mol is here employed, based on the spectrum of experimental data available for the trapping energy at grain boundaries [17,21].

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 Figure 2: Flowchart that summarises the methodology adopted in the numerical deformation-diffusion-fracture analysis of polycrystalline solids exposed to hydrogen.
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-diffusionfracture 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.

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  Crack nucleation and growth in a square plate microstructure of 50 grains is investigated. The plate, with dimensions given in Fig. 3a, 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 ofu 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 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.  The results obtained for the case of C env = 0 ppm wt are shown in Figs. (3b-3g), 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.  Contrarily, when the sample is exposed to hydrogen, cracking initiates along the grain boundaries and propagates in an intergranular manner. This is shown in Fig. 4, where the deformation and separation of the grains is shown as a function of the applied displacement for the case of

Virtual slow strain rate testing on Monel K-500 samples
Recent slow strain rate tests (SSRTs) on a Ni-Cu superalloy (Monel K-500) have revealed significant intergranular cracking depths, much larger than those expected from diffusion calculations for relevant exposure times [38]. Hence, these experiments suggest that crack initiation is likely to take place much before final failure, allowing for the hydrogen-containing solution to reach significant depths by following the crack. Early sub-critical crack growth would imply the need for a fracture analysis, compromising the suitability of the SSRT experiment and its metrics (e.g., time to failure) in assessing hydrogen embrittlement susceptibility. The micromechanical formulation presented can be used to gain complementary insight into these paradigmatic experiments and the early cracking hypothesis. Mimicking the testing conditions, notched cylindrical specimens with the dimensions shown in Fig. 6a are considered in our simulations. The samples are subjected to uniaxial loading and thus one can take advantage of axially symmetric conditions. Accordingly, the sample is discretised using axisymmetrical finite elements; a total of 64,835 quadratic elements are employed. A microscopic region is introduced near the notch, see Fig. 6b. This region spans a width of 0.5 mm and includes 280 grains, while the remaining part of the solid is modelled as an isotropic continuum without interfaces. In this way, the model can capture the two cracking modes observed in the experiments; in the absence of hydrogen, cracking occurred at the centre of the sample due to plastic instabilities, while in the presence of hydrogen, cracking took place near the notch tip and was of intergranular nature [38].   Each Monel K-500 heat was tested in four different environments: in air (i.e., in the absence of a hydrogen-containing solution) and while being exposed to solutions with the applied potentials The diffusible hydrogen concentration associated with each charging condition was experimentally measured and used as input in the model -a prescribed C env magnitude at the surface. The values measured are given in Table 3 (in ppm wt). Heat Simulations are first conducted in the absence of hydrogen. As shown in Fig. 7, damage initiates at the centre of the sample for sufficiently large loads, in agreement with experimental V SCE are shown in Fig. 10. In all cases the failure is intergranular, starting at the notch tip and triggering a significant drop in the load carrying capacity before any noticeable increase in the phase field variable is observed. Thus, the calibrated model is also able to capture the ductile-to-brittle transition observed in this case study. Interestingly, the cracking event appears to occur rather suddenly, with the intergranular crack propagating into regions with low hydrogen content. This would suggest that the crack resulting from the decohesion of the grain boundaries exposed to a high hydrogen content would be sufficiently large to propagate in an unstable fashion through grain boundaries that have only been negligibly weakened by hydrogen. If that were to be the case, then these numerical results would suggest that SSRT is not compromised by early cracking and thus remains a valid test for measuring hydrogen embrittlement susceptibility. Additional, albeit limited, insight can be gained by a simple estimate of the transition flaw size a t = K 2 Ic /(πσ 2 Y ). For Monel K-500 exposed to a relatively uniform hydrogen distribution resulting from an applied potential of E A = −1.1 V SCE , the transition flaw size would be on the order of 0.04 mm (see Table 2 and the K T H measurements of Ref. [7]). However, this quantity can increase to up to 4 mm in the absence of hydrogen. Thus, the magnitude of a t relevant to this scenario (a non-uniform distribution of hydrogen) falls between those two limiting cases, and could therefore be higher or lower than the crack extensions predicted (0.1-0.5 mm). Another source of uncertainty is the specific traction-separation law adopted, as assuming the existence of a damage dissipation region could add an additional source of fracture resistance. The interface parameters that provide a quantitative agreement with experiments are given in Table 4. The damage coefficient χ was estimated based on previous (microstructurally-insensitive) phase field simulations [38]. The values of χ used are higher than those estimated using atomistic calculations for most common types of Ni grain boundaries [60]. However, the choices of χ values are notably sensitive to the magnitude of the grain boundary binding energy considered in (23), and the estimation of this magnitude carries a degree of uncertainty [63]. The calibrated model is used to reproduce the entire experimental campaign, conducting virtual experiments on the four material lots over the four environments considered. We summarise As observed, the model is capable of predicting the occurrence of hydrogen-assisted brittle failures in all but one case -the experiment on the NRL LS heat at E A = −0.85 V SCE . As discussed in Ref. [38], this is a rare case as SEM images of the fracture region reveal intergranular features but the time to failure happens to be larger than that measured in the absence of hydrogen. It is thus concluded that the micromechanical model presented is capable of predicting hydrogen embrittlement upon appropriate calibration.  showed that embrittlement (hydrogen-assisted intergranular cracking) occurred even at cryogenic temperatures, where dislocation-hydrogen interactions are precluded. This suggests that hydrogenassisted decohesion of grain boundaries is a first-order mechanism in hydrogen embrittlement.
It was also found that intergranular microcrack evolution was enhanced at room temperature, relative to 77 K, but a mechanistic interpretation of this finding was deemed complicated due to the multiple factors at play. Here, we examine the ability of our micromechanical model to quantitatively reproduce the seminal experiments by Harris et al. [4] and use the numerical insight provided to gain further understanding on the role of temperature in enhancing embrittlement.
The material properties adopted correspond to those reported by Harris et al. [4], which are listed in Table 6 as a function of the temperature and the environment. Two environmental conditions are considered: (i) samples tested in air, without any hydrogen pre-charging, and (ii) samples exposed to a hydrogen content of 4000 appm (79.5 wppm). In the latter, gaseous hydrogen charging is used and the hydrogen is distributed uniformly within the samples. Mimicking the experimental conditions, the uniaxial load is prescribed by applying a remote vertical displacement with a rate ofu y = 0.0078 mm/s, while the bottom edge is completely constrained (u x = u y = 0).
The samples are cylindrical bars with the dimensions given in Fig. 11a. As in the previous case study, we take advantage of axial symmetry and model half of the 2D section using axisymmetric elements. A total of approximately 80,000 quadratic axisymmetric elements are used to discretise the model. As shown in Fig. 11b, we introduce a microstructural domain of 200 grains in the   Table 7 and, following the approach of Ref. [16], a phenomenological degradation law is adopted, such that with the Gibbs free energy being equal to 17 kJ/mol [64].
As shown in Figs. 11 and 12, the model is able to reproduce experimental measurements beyond the regimes of calibration, both qualitatively and quantitatively. Consider first the cracking patterns shown in Fig. 11. In the absence of hydrogen (Fig. 11c) parameter. However, when the sample is exposed to hydrogen, then cracking takes place in an intergranular fashion, as a result of the failure of the cohesive zone interfaces. The location of the grain boundary decohesion event that nucleates the brittle crack is random. For the microstructure and conditions of Fig. 11d, it occurs close to the edge of the sample, with the crack growing then both towards the outer surface and towards the centre of the sample. This change from ductile transgranular damage to brittle intergranular cracking due to hydrogen is observed at both 77 and 298 K, as in the experiments.
The quantitative results obtained for the four scenarios are shown in Fig. 12, in terms of the predicted and measured engineering stress-strain curves. A very good agreement with the experiments is observed. Interestingly, the good agreement observed for the case of the hydrogencharged sample suggests that the higher degree of embrittlement observed at room temperature can be rationalised by the additional accumulation of hydrogen at grain boundaries due to diffusion, without the need for additional contributions from mechanisms such as those associated with hydrogen-deformation interactions.

Summary and concluding remarks
We have presented a new microstructurally-sensitive deformation-diffusion-fracture formulation to predict hydrogen embrittlement in elastic-plastic solids. The modelling framework is able to explicitly resolve the microstructure and capture the transition from ductile transgranular fracture to hydrogen-assisted brittle intergranular cracking. This is achieved by combining a phase field description of transgranular cracks with a cohesive zone model that simulates decohesion at the grain boundary interfaces. Furthermore, mechanistic insight into the embrittlement of polycrystalline Ni is gained, showing that grain boundary decohesion is a first order effect, and that the differences between the responses observed at ambient and cryogenic temperatures can be rationalised by the additional content of hydrogen accumulated in grain boundaries due to diffusion. The numerical experiments conducted also showcase the computational robustness of the method, with significant cracking being predicted without convergence issues.