A modelling framework for coupled hydrogen diffusion and mechanical behaviour of engineering components

In this paper, we propose a finite element formulation for solving coupled mechanical/diffusion problems. In particular, we study hydrogen diffusion in metals and its impact on their mechanical behaviour (i.e. hydrogen embrittlement). The formulation can be used to model hydrogen diffusion through a material and its accumulation within different microstructural features of the material (dislocations, precipitates, interfaces, etc.). Further, the effect of hydrogen on the plastic response and cohesive strength of different interfaces can be incorporated. The formulation adopts a standard Galerkin method in the discretisation of both the diffusion and mechanical equilibrium equations. Thus, a displacement-based finite element formulation with chemical potential as an additional degree of freedom, rather than the concentration, is employed. Consequently, the diffusion equation can be expressed fundamentally in terms of the gradient in chemical potential, which reduces the continuity requirements on the shape functions to zero degree, C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}_{0}$$\end{document}, i.e. linear functions, compared to the C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {C}}_{1}$$\end{document} continuity condition required when concentration is adopted. Additionally, a consistent interface element formulation can be achieved due to the continuity of the chemical potential across the interface—concentration can be discontinuous at an interface which can lead to numerical problems. As a result, the coding of the FE equations is more straightforward. The details of the physical problem, the finite element formulation and constitutive models are initially discussed. Numerical results for various example problems are then presented, in which the efficiency and accuracy of the proposed formulation are explored and a comparison with the concentration-based formulations is presented.


Introduction
In the presence of hydrogen, many materials experience significant degradation of their mechanical properties, which is commonly known as hydrogen embrittlement (HE) [1][2][3]. These detrimental effects are reported in the literature in terms of reductions in ductility, strength and toughness. As a result, HE may result in premature catastrophic failure at load levels much lower than in the hydrogen free case. A number of mechanisms of embrittlement have been proposed in the literature [2,3]. The two mechanisms that have received the most attention are hydrogen induced (or enhanced) decohesion (HID or HEDE) [4,5] and hydrogen enhanced local plasticity (HELP) [6,7]. In the HEDE mechanism [8][9][10], hydrogen accumulates at an interface and reduces the cohesive strength and hence the energy required to fracture the material. In the HELP mechanism [11,12], hydrogen is assumed to enhance the dislocation mobility, change the core energy and also reduce the elastic interaction between dislocations, generally termed elastic shielding, each of which can result in material softening and localisation of plastic flow. In recent years, there has been significant interest in developing computational models for hydrogen embrittlement. Although there is still significant controversy and debate in the scientific literature concerning the detailed mechanisms responsible for embrittlement, it is generally agreed that any model should contain the following ingredients: -Diffusion of hydrogen through the material and hydrogen accumulation within different microstructural features. -The effect of hydrogen on the plastic response of the material. -The effect of hydrogen on the cohesive properties of any interfaces in the material.
Most published finite element models of the embrittlement process follow the approach originally proposed by Sofronis and McMeeking [13] in which the governing equation for diffusion is expressed in terms of the gradient in concentration and gradient in the mean (or hydrostatic) stress, and the nodal degrees of freedom for the coupled mechanical/diffusion problem are taken as the displacements and concentration. As noted by Barrera et al. [14], the dependence of the diffusion equation on the stress gradient requires the use of finite elements with shape functions that are at least quadratic in terms of displacement. Also, for situations where interfaces are modelled, either between similar or dissimilar materials, there can be a discontinuity in concentration at an interface. This can lead to numerical problems. In this work, we return to fundamentals and express the diffusion equation in terms of the gradient in chemical potential and we take the chemical potential as a degree of freedom, rather than the concentration. A consequence of this is that the finite element formalism no longer requires stress gradients to be determined, allowing the use of elements with linear shape functions and simplifying the finite element implementation (particularly when implemented in finite element packages such as Abaqus [15]). Also, the chemical potential is a continuous function within the body, and we do not need to confront problems arising from discontinuities in our solution dependent variables, particularly at interfaces.
In the following section, we provide a more detailed background to the class of physical problems considered here and a general statement of the combined mechanicaldiffusion analysis of multi-component and/or multi-grained boundary value problems. This is followed by a statement of the fundamental diffusion and mechanical equations for the bulk material and interfaces. The governing finite element equations are developed in Sect. 3, using a Galerkin approach based on statements of conservation of mass and mechanical equilibrium. Implementation of the finite element approach further requires constitutive models for the mechanical response of the matrix and interfaces. Suitable models are presented in Sect. 4 and these are used in Sect. 5 to analyse some simple representative problems, where we compare the efficiency and accuracy of the approach with that described by Barrera et al. [14].

Background
In this study, we investigate a class of problems in which a material or multiple materials are subjected to mechanical loading in the presence of hydrogen, i.e. combined chemicalmechanical loading. In particular, we are interested in the response of polycrystalline metallic engineering alloys. Generally, hydrogen diffuses in metals through interstitial lattice sites (NILS). Dislocations and interfaces can provide additional paths for diffusion, but these sites tend to trap hydrogen (as discussed further below) and diffusion along them can be much more sluggish than through the matrix [16,17]. We therefore ignore these minor contribution here. Hydrogen resides in NILS and the trapping sites, i.e. dislocations, grain boundaries, carbide/matrix interfaces, microvoids and other defects. The NILS are considered to be the dominant sites. However, the effect of shallow traps on diffusivity is very significant and deep traps can play an important role acting as a gutter for hydrogen and therefore cannot be ignored. Under mechanical loading, the hydrostatic stress gradient also plays an important role in determining the flux of hydro-gen within a body. Moreover, in elastic-plastic materials, the number of trapping sites increases as the dislocation density increases with increasing plastic deformation. Hence, the chemical and mechanical fields are strongly coupled and the material behaviour is controlled by the characteristic time of the diffusion and deformation processes. Further, depending on the length scale and density of traps, they can either be treated as a continuous distribution in the material, such that their presence modifies the effective diffusivity in the matrix, or as discrete features, such as grain boundaries, carbide/matrix interfaces, precipitate, etc., wherein a direct interaction between the trap and mobile hydrogen is considered. We consider each of these here. In the following sub-section, we introduce a model problem that can be used to represent different material systems such as grains and their boundaries, precipitates in a matrix or any other situation where failure can occur at an interface.

The model problem
We consider a body which is comprised of two different materials A and B that can have different mechanical and diffusional (and thermodynamic) properties, which are separated by an interface, see Fig. 1. A common Cartesian coordinate system for the reference and deformed configurations X i and x i , i = 1, 2, 3, respectively, is assumed. The body occupies a volume Ω with boundary Γ -materials A and B occupy Ω A and Ω B respectively, i.e. Ω = Ω A ∪ Ω B , in the current configuration. The interface is defined by Γ int . Failure can occur along the interface as the material deforms, leading to complete separation of the interfacial planes either side of the interface. To allow for this, a material point along the interface is defined by the two normal vectors n − i and n + i either side of the interface, as shown in Fig. 1ii, where n − i = −n + i , i.e. the initially intact material point splits into two points with unit normals acting opposite to each other and into the material on either side of the interface (on the upper and lower interface surfaces Γ ± int of Fig. 1ii, respectively). The body is assumed to be subjected to mechanical loading and contains hydrogen that may diffuse through the body and across interfaces and boundaries. The total concentration of hydrogen in the bulk, i.e. materials A and B, at a material point is defined by the total number of H atoms in NILS and traps per unit volume (H atoms/m 3 ). The traps in the bulk are assumed to be continuous and locally in equilibrium with mobile hydrogen following Oriani's theory [18]. Hence, the total concentration in materials A and B, C A T and C B T , respectively, can be written as where C L is the hydrogen concentration in the lattice, C Tr is the concentration associated with the traps and the superscript Fig. 1 The schematic of the model problem. The figures illustrate: i a continuum body that is subjected to macroscopic deformation E i j , Σ i j in the presence of hydrogen that is defined by the chemical potential and concentration at the lattice and discrete traps (μ, C L , C Tr ); and ii the microstructure is comprised of a bi-material system of materials A and B that are separated by an interface where the microscopic deformation is ε i j , σ i j and the hydrogen chemical potential and concentration in the lattice, distributed and discrete traps indicates the material. The molar concentration (mol/m 3 ) are where N A = 6.0232 × 10 23 atoms/mol is Avogadro's number and a bar over the symbol indicates a molar quantity. The interface is considered as a discrete trap with a concentration-now defined as the total number of H atoms per unit area of the interface (similarly, the molar concentration in the trap isC int = C int /N A ). It should be noted that the hydrogen concentration is not continuous across the interface, i.e. the concentration in the bulk materials is different to the concentration at an adjacent point in the interface. Furthermore, the hydrogen concentration in the interface is assumed to be in equilibrium with the mobile hydrogen in the bulk adjacent to the interface. We limit our consideration to the situation where hydrogen can cross the interface, i.e. in the normal direction, n, but not along the interface direction, t. This is consistent with the recent Monte Carlo simulations of Du et al. [17]. The chemical potential μ (J/mol) is assumed to be continuous at the interface, i.e. it has the same value in the matrix and the interface. Diffusion across the interface is allowed and this is assumed to be driven by a small difference in chemical potential of hydrogen in materials A and B either side of the interface. Hence, we effectively consider the interface to have a small thickness h, such that the chemical potential μ int is assumed to be given by a continuous function where μ ± are the chemical potentials at the upper and lower interface surfaces Γ ± int in the bulk materials, respectively, see 1ii. In practice, when we assign a chemical potential to the interface, this is taken as the mean of the values either side of the interface. Further, the diffusion distance for transport across the interface is small (of the order of atomic dimensions) and we assign a large kinetic constant for the diffusion constant, so that in the simulations μ + ≈ μ − . This is discussed further below. Materials A and B are assumed to exhibit elastic or elastic-plastic behaviour. Thus, under mechanical loading, the chemical potential may change due to the change of hydrostatic stress, allowing the hydrogen to diffuse. Also, the density of traps can increase as dislocation density increases during plastic flow, which results in an increase in concentration associated with the traps, C Tr , and a concomitant reduction of the concentration in the lattice. Further, swelling of the bulk material (i.e. volumetric deformation) may result from the hydrogen in the NILS and/or traps. The interface is assumed to be defined by a traction-separation law that undergoes full separation after reaching a critical value. The chemical potential at the interface changes under deformation allowing the concentration to increase or decrease. Moreover, the presence of hydrogen may reduce the cohesive strength of the interface and induce a separation (jacking) that may influence the deformation in the bulk adjacent to the interface and lead to a faster separation. Our representation of this process is guided by the recent atomistic simulations of Katzarov and Paxton [19] and is consistent with the thermodynamic framework described by Mishin et al. [20].

Hydrogen diffusion
The mass conservation of hydrogen states that the rate of change of the total hydrogen concentrationC T in a volume Ω is equal to the fluxJ i through its boundaries Γ . The swelling and elastic deformations are assumed to be small and only plastic deformation is taken to be large which is essentially isochoric and therefore the total rate of change of the hydrogen concentration is equal to its local rate of change (for further details see "Appendix A"). Thus, it can be written as where dV and dS are the infinitesimal volume and surface elements in the current configuration, respectively, and n i is the unit outward normal on Γ . Applying the divergence theorem, we obtain The molar hydrogen fluxJ i (mol/(m 2 · s)) in the bulk is related to the gradient of the chemical potential μ (J/mol) as where D L is the lattice diffusivity (m 2 /s). Hydrogen is allowed to diffuse across the interface and not along the interface, i.e. hydrogen diffusion within the interface is only in the n-direction. The hydrogen flux across the interface is written as whereJ int (mol/(m 2 · s)) is the flux in the direction normal to the interface, h int is the hydrogen transfer coefficient across the interface (a kinetic constant for the interface) (mol 2 /(m 2 · s · J)) and Δμ int = μ + − μ − is the chemical potential jump across the interface. The boundary conditions of the diffusion equation for the model problem are defined bŷ μ = μ on Γ μ , whereμ andĴ are the applied chemical potential and hydrogen flux on the boundaries Γ μ and Γ J , respectively, and n i is the outward normal to Γ . The boundaries Γ μ and Γ J should satisfy The conditions on the interface Γ int are defined as: where μ int is the chemical potential in the interface Γ int and the hydrogen flux is given byJ int =J i n int i where n int i is the outward normal vector to Γ int . It should be mentioned that the interface is assumed to undergo complete separation when it fails such that the newly created free surfaces have flux boundary conditions specified. Hence, we introduce a scalar damage variable d that is assumed to evolve monotonically from 0 to 1, i.e. from an undamaged to a fully damaged state. It follows that the no flux boundary condition is defined as

Mechanical equilibrium
The balance of linear momentum, i.e. the mechanical equilibrium equations, in a quasi-static state and the absence of body forces can be written as where σ i j is the Cauchy stress. The boundary conditions for the model problem is defined bŷ whereû i andt i are the applied displacement and traction values acting on Γ u and Γ t , respectively, that should satisfy The conditions on the interface Γ int are defined as: where Δu int i is the displacement jump across the interface Γ int , u ± i are the displacement vectors in the upper and lower interface surfaces Γ ± int , respectively, and the traction acting on the interface is given by t int i = σ i j n int j . In the case of complete separation, the initially intact interface becomes two traction free surfaces that are defined by

Finite Element implementation
In this section, we present the implementation of the model problem in Sect. 2. We consider the coupled diffusionmechanical system given by the mass and linear momentum balance equations in Eqs. (4) and (11), respectively. These equations are usually expressed in terms of two primary variables that are the displacement field u i and the total hydrogen concentrationC T . Alternatively, we propose the use of the chemical potential μ as the diffusion degree of freedom instead of the total hydrogen concentrationC T . In most coupled diffusion-mechanical analysis, researchers use the concentration as a degree of freedom, e.g. Sofronis and McMeeking [13], Barrera et al. [14] and Oh and Kim [21]. Di Leo and Anand [22] have proposed a theory for coupled diffusion and large elastic-plastic deformations using chemical potential as a degree of freedom which they later used to study fracture of ferritic steels in the presence of hydrogen using a continuum damage approach (Anand et al. [23]). Further, the chemical potential has been used in similar problems such as coupled diffusion and mechanical deformation in elastomeric gels (Chester and Anand [24]) and in Li-ion batteries (Miehe et al. [25]). In comparison with these formulations, a broader description of material systems is consider here by including the behaviour of interfaces. We propose a thermodynamically consistent framework based on the work of Mishin et al. [20], see Sect. 4. The proposed formulation has many advantages in comparison with the standard formulation and a detailed comparison between these formulations will be discussed later in this section. In the following, we present the details of the proposed Finite Element formulation of the coupled diffusion-mechanical analysis in the presence of interfaces.

The variational formulation
The variational forms of the governing equations are required for the development of the Finite Element equations. Consider the mass balance equation in Eq. (4), and a virtual variation of the chemical potential δμ, which is continuous and is assumed to vanish where the chemical potential is defined on a boundary, i.e.
Multiplying Eq. (4) by δμ and integrating over the volume gives the variational form of the mass balance equation Thus, integrating by parts and applying the divergence theorem yields Similarly, an admissible virtual displacement field, δu i , is considered which vanishes on boundaries where the displacement is prescribed, i.e.
Now, the variational form of the equilibrium equation of (11) takes the form and integration by parts and application of the divergence theorem gives

Time integration
A fully implicit backward Euler scheme is employed here for time integration of the diffusion equations. Thus, the total concentrationC T at time t p+1 is approximated bȳ where indices p + 1 and p denote variable values at instants t p and t p+1 , respectively, and Δt = t p+1 − t p is the time increment.

The spatial discretisation
In the discretisation of the body in Fig. 1, two regions are considered: the bulk (materials A and B); and the interface regions. We use continuum and cohesive elements to discretise the bulk and interface, respectively. We adopt the standard Galerkin method where the primary variable and their variations are discretised using the same shape functions, see for example Belytschko et al. [26] and Wriggers [27]. It should be noted that different interpolation may be used for the continuum and cohesive elements. Therefore, fields for the chemical potential and its variation are discretised in the bulk and interface regions as and respectively, where the superscript (•) int is used to distinguished the variables associated with the interface from those in the bulk (no superscript is used for the bulk), N I are the standard finite element shape functions, N node are the total number of nodes in the mesh, μ I and δμ I are the values of the chemical potential and its variation at node I , respectively, and a node is considered to be in the bulk if I ∈ Ω and at the interface if I ∈ Γ int . The displacement field and its variation are interpolated in the same manner as and where u i I and δu i I are the values of the displacement fields and its variation at node I , respectively.

The Finite Element equations
In order to derive the Finite Element equations, we need to determine the discrete forms of the diffusion and mechanical equilibrium equations. This can be done by substituting the interpolations of the chemical potential, displacement and their variations, and the time integration into the variational forms of the diffusion and mechanical equilibrium equations. Thus, substitution of the time derivatives of the total concentration in Eq. (22) and chemical potential and its variation in Eqs. (23) and (24), respectively, into Eq. (18) gives the discrete form of the diffusion equation as the variation of a functional Π μ : where the first sum is over the nodes of Ω, the second sum is over the nodes on Γ J and the third and fourth sums are over the nodes on Γ int . This variation can be written in terms of the nodal chemical potential variation and diffusion forces as where the nodal chemical force in the continuum elements is where A int I = H I N int I and H I is a multiplier that relates the interface chemical potential to the bulk nodal chemical potential as in Eq. (2), i.e. μ int I = H I μ I . Since the chemical potential variation δμ I can be chosen arbitrarily, the residual of the diffusion equation at node I is then determined as Similarly, substitution of the test and displacement fields in Eqs. (25) and (26) into Eq. (21) gives the discrete form of the mechanical equilibrium equation as the variation of a functional Π u : where B j I = ∂ N I /∂ x j . This variation can be written in terms of the nodal displacement variation and mechanical forces as where the nodal mechanical force in the continuum elements is B int I = M I N int I and M I is a linear multiplier that relates the displacement jumps to the nodal displacements, i.e. Δu int i = M I u i I . Since δu i I can be chosen arbitrarily, the components of the residual of the mechanical equilibrium equation at node I are determined as To this end, we need to solve the coupled nonlinear diffusion-mechanical discrete equations, i.e. Eqs. (30) and (34), for the nodal chemical potential μ I and displacements u i I . In the next sub-section, the solution procedure is discussed.

Solution procedure
In order to solve the nonlinear set of discrete equations represented by (30) and (34) with a Newton-Raphson algorithm, a suitable linearisation is required. Thus, an iterative scheme can be determined using a first-order Taylor expansion of the residuals, i.e. the chemical and mechanical nodal residuals at node I are R μ I and R u I , respectively. Assume that the nodal chemical potential and displacement components at node J at the start and end of an iteration r are μ J , u j J and μ J + Δμ J , u j J + Δu j J , respectively. Hence, the residuals and unknown values can be written in vector form as respectively. Therefore, the linearised form of the residuals becomes where Δd J is the nodal increment vector of the unknown quantities, O Δμ 2 J , Δu 2 j J , . . . represents the higher order terms in vector form and the superscript denotes the iteration. Hence, by neglecting the higher order terms, the incremental form of these equations can be written as where K I J is the tangent stiffness matrix that is given by and the elements of the tangent stiffness matrix, i.e. the tangent operators, are The details of these operators are provided in "Appendix A". It is worth mentioning that the residual and tangent matrix can be split into the bulk and interface contributions. In this study we use continuum and cohesive elements to implement the diffusion and mechanical models in the commercial FE code Abaqus. In particular, we use the user-defined subroutines UMAT or UHARD and UMATH for the bulk and UEL for the interface element. The implementations are provided in "Appendix B".

The concentration versus chemical potential based formulations
To compare the two formulations, we recall the discrete form of the diffusion equation in Eq. (27). In the concentration based formulation, the total concentration in the bulk,C T , and the trap concentration in the interface,C int , are given as whereC T I andC int I are the values of the total concentration and trap concentration at node I , respectively. Hence, the linearisation of the discrete diffusion equation yields a system of equations similar to Eq. (35) where the chemical potential increment Δμ I is replaced by the total concentration and interface concentration increments ΔC T I and ΔC int I , respectively. As a result, the derivative ∂J /∂μ in the tangent operator K μ μ I J is replaced by ∂J /∂C T , see "Appendix A". These derivatives are determined from the definition in Eq. (5) as and Examining the continuity requirements on the shape function, in the chemical potential based formulation, Eq. (40) implies that the shape functions in Eq. (23) should exist and be continuous, i.e. they satisfy C 0 continuity. On the other hand, in the concentration based formulation, the gradient of the chemical potential contains the gradient of hydrostatic stress ∂σ h /∂ x i , see Sect. 4. Hence, the strain gradient should now exist and be continuous. Consequently, the shape functions in Eq. (25) should be continuously differentiable, i.e. satisfy C 1 continuity, which means that shape functions that are at least quadratic should be used. Therefore, the chemical potential based formulation has less restrictive continuity requirements which allow linear approximations to be used. This may result in a lower computational cost. The spatial discretisation of the model problem requires that the top and bottom faces of the cohesive elements to be attached to the continuum elements, see Fig. 2. Hence, in the case of a concentration based formulation, the concentration is defined at the nodes of these faces. Thus, the concentration will be discontinuous at each node bearing two different values from the bulk and interface. This inconsistency can be overcome using the chemical potential based formulation. The chemical potential is essentially continuous across the interface which leads to a consistent definition of its values at these surfaces. Further, the concentration at the interface is defined in the mid-surface of the cohesive element (the constitutive surface) which allows the concentration to differ from that in the bulk.

Constitutive descriptions
In this section, we present the hydrogen transport and mechanical constitutive models for the bulk and interface.

The Bulk model
Consider a material in which the hydrogen concentration in the lattice is C L and that from the distributed traps is C Tr . Thus, the total hydrogen concentration is C T = C L + C Tr . The hydrogen concentration in the lattice is given by Sofronis and McMeeking [13] where β L is the number of interstitial sites per solvent atom, θ L ∈ [0, 1] is the fraction of lattice sites occupied by hydrogen atoms and N L = N A /V M is the number of atoms of solvent per unit volume (atoms/m 3 ). The molar volume of the solvent lattice is V M = M/ρ, where ρ and M are the density and relative atomic mass of the lattice element, respectively. We consider n Tr types of traps such that the total concentration in the traps is given by where C i is the concentration in the ith trap that is defined by where α i = 1 denotes the number of trapping sites, N i (atoms/m 3 ) denotes the number of atomic trapping sites per unit volume and θ i ∈ [0, 1] is the fraction of trapping sites occupied by H atoms. We consider two types of traps-those in which the number of trap sites are fixed (fixed traps i ≡ F) and those that increase with increasing plastic deformation (dislocation traps, i ≡ D). The molar concentration of the ith trap is denoted asC i = C i /N A . Hence, N F is taken to be constant and N D is assumed to be proportional to the dislocation density ρ D , which is assumed to be a function of the accumulated equivalent plastic strain ε pl e : where a is the lattice constant and the prefactor is 1/b, where b is the magnitude of the a/2 <1 1 0> Burgers vector in FCC materials, which is a reasonable approximation for the atomic spacing along a dislocation line. The traps are modelled using Oriani's theory, in which the hydrogen concentration in the lattice is assumed to be in equilibrium with the concentration in the traps. This means that the concentration in the traps can directly be estimated from the concentration in the lattice. Hence, for finite populations of H in the lattice and trap sites, the equilibrium condition for the ith trap is given by where K i = exp (−W i /RT ) = const., represents equilibrium between the lattice and the ith type of trap site and W i < 0 is the binding energy for the ith trap, R is the universal gas constant and T is the absolute temperature. It should be noted that Eq. (46) reduces to θ i / (1 − θ i ) = K i θ L in the case of small lattice concentration, i.e. θ L 1, which is generally the situation in practice (note, however that the corresponding occupancy in the traps can approach 1, particularly for high trap binding energies [28]). The Gibbs free energy per unit volume of the bulk material, G, can be written as where ψ is the free energy density per unit volume of the bulk material, σ i j is the applied macroscopic stress, ε i j is the corresponding strain, μ 0 represents the chemical potential at a suitably defined standard condition, T is the absolute temperature and s is the entropy per volume in the bulk material, which we assume is dominated by the configurational entropy (i.e. entropy of mixing): whereC max L is the hydrogen concentration corresponding to the saturated state in the lattice. The free energy density for elastic and elastic-plastic materials can be expressed as where ε pl i j is the plastic strain, ε s i j = ε s δ i j = 1 3 V MCL δ i j is the swelling strain (here we have assumed that swelling is only associated with hydrogen in the lattice and that trapped hydrogen does not contribute to swelling observed macroscopically-this assumption can be readily relaxed if necessary, i.e. see "Appendix A)" and C i jkl is the 4th order elasticity tensor. It is worth noting that both the swelling and elastic deformation are assumed to be small.
The chemical potential in the bulk is determined from the Gibbs free energy in Eq. (47) as It should be noted that at low concentrations (C L C max L ), the chemical potential becomes The mechanical response is determined by assuming that the bulk is in mechanical equilibrium, i.e. by solving ∂G/∂ε i j = 0 to give σ i j . Note, when considering two different materials, e.g. A and B of Fig. 1b, μ 0 can be different for each material. This needs to be taken into account when solving any boundary value problem -particularly when invoking continuity of chemical potential across an interface. In the following, for consistency of notation, we take H in one of the materials, say A, as the reference, for which we designate Oriani's relationship, Eq. (46), can then be used to determine the occupancy in B in terms of that in A at interfaces where the chemical potential is required to be continuous.

The interface model
The interface is treated as a discrete trap. The description here parallels that for the bulk described above in Sect. 4.1 and is consistent with the thermodynamic description of Mishin et al. [20]. We assume that the hydrogen concentration, C int (atoms/m 2 ), in the interface can differ from that in the bulk (i.e. the interface binding energy can be different from the bulk) and is defined by where α int = 1 denotes the number of H trapping sites per interface site, N int (atoms/m 2 ) denotes the number of atomic trapping sites per unit area of interface and θ int ∈ [0, 1] is the fraction of trapping sites occupied by H atoms. In this study, we limit our consideration to situations where failure at an interface is driven by the component of the traction normal to the interface, allowing the tangential response to be given by a stiff linear mechanical model. We can then simply concentrate on the response under the action of a normal traction to derive the necessary thermodynamic equations. The Gibbs free energy per unit area of the interface, G int , can then be written as where φ is a free energy density per unit area, T n is the traction in the normal direction, δ n is the separation in the normal direction. W int is the excess cohesive/interface zone binding energy compared to that in the reference matrix material (e.g. material A of Fig. 1ii) and s int is the configurational entropy per unit area of interface, and is given by whereC max int is the hydrogen concentration corresponding to the saturated state in the cohesive zone. The constitutive response for the cohesive interface is modelled in terms of the relationship between the traction and the displacement jump across the interface. The presence of hydrogen in the cohesive zone may introduce a volumetric change that can be expressed in terms of a jacking or swelling separation δ s n [19] that can be defined as where V int M is the molar volume of hydrogen in the interface zone (m 3 /mol) which can have a different value to V M in the bulk. Hence, the total separation in the normal direction δ n can be decomposed into mechanical and jacking components: where δ m n is the separation in the normal direction due to mechanical loading. Note that φ in Eq. (53) is a function of δ m n , i.e.
The chemical potential in the cohesive zone is determined from the Gibbs free energy in Eq. (53) as The traction-separation law (TSL) is determined from the Gibbs free energy assuming that the interface is in mechanical equilibrium, i.e.
Making use of Eqs. (56), (57) and (59), Eq. (58) becomes The free energy density φ can be physically based or a phenomenological model that describes the relation between the traction and the displacement jump across the cohesive zone. The initial behaviour is assumed to be linear elastic until the onset of damage. When damage is initiated, the material stiffness decreases and the degradation is defined by the scalar damage variable d that evolves monotonically from 0 to 1, i.e. from an undamaged to a fully damaged state, and can take different forms, i.e. linear, exponential, trapezoidal etc. The irreversible traction-separation law is then written as where K n is the initial elastic or penalty stiffness and δ c n is the critical separation in the normal direction. The linear elastic unloading behaviour is determined by the degraded stiffness K n = (1 − d) K n . A bilinear traction-separation law is used in this study. Hence, the damage variable, d, for linear softening is determined as where δ max n is the maximum value of the displacement attained during the loading history and δ f n is the separation at failure.
The initial cohesive strength in the normal direction is defined by σ c0 = K n δ c0 n where δ c0 n is the hydrogen free critical separation. The cohesive strength is assumed to decrease with increasing hydrogen concentration. Hence, we assume that the softening takes the general form where Φ is a softening function. For example, a linear form is given by Barrera and Cocks [29] as where ξ int ∈ [0, 1] is a hydrogen softening parameter,C min int is the concentration at the onset of softening andC max int is the maximum concentrations associated with softening. (The form of the softening function can be determined by fitting of experimental results (e.g. Raykar et al. [30] and Martínez-Pañeda et al. [31]) or using atomistic simulations (e.g. [32] and Katzarov and Paxton [19].) It is worth mentioning that this study is limited to pure mode-I loading; the damage evolution is only considered in the normal direction and the response of the tangential direction is taken to be purely elastic as where T t and δ t are the traction and total separation in the tangential direction.

Numerical examples
In this section, we validate and demonstrate the effectiveness of the proposed formulation through studying two problems. Firstly, we consider a 2D fully coupled elastoplastic diffusion problem of a deep notched specimen in the absence of damage. This problem is analysed using both concentration and chemical potential based formulations and a detailed comparison showing the similarities in the solutions is provided. Additionally, we investigate the failure of the notched specimen using the interface model described in Sect. 4.2. The second problem is a micromechanical investigation of hydrogen-enhanced decohesion (HEDE) in a dissimilar metal weld. More specifically, we consider the failure in the presence of hydrogen of the carbide-matrix interfaces around fine M 7 C 3 precipitates generated adjacent to the interface of an 8630 steel/IN625 nickel alloy dissimilar weld, which has been analysed extensively by Barrera et al. [33].

Analysis of deformation, diffusion and failure in a deeply notched specimen
A plane strain plate with deep double-edged notches is considered as shown in Fig. 3. The width, height and thickness of the specimen are denoted by W , 2 H and B, respectively, and the notch radius, depth and angle are denoted by r , d and ϕ, respectively. The plate is made of a material that is assumed to exhibit elastic-plastic behaviour. The presence of hydrogen in the lattice and dislocation traps are considered. An interface model is used to simulate failure such that damage may accumulate at the interface leading to loss of load carrying capacity. Further, hydrogen can be trapped at the interface which is treated as a discrete trap. We assume a full coupling between the mechanical and diffusion responses and the constitutive descriptions are taken to be suitable for a high strength steel. The purpose of this analysis is to: (1) compare the concentration and chemical potential formulations; and (2) investigate the failure of the specimen using the chemical potential formulation. The deep double-edged notch geometry provides a high degree of plastic constraint and high level of hydrostatic stress across the minimum section which allows us to compare the two formulations for the fully coupled diffusion and elastic-plastic deformation processes. Furthermore, notched specimens are commonly used to study hydrogen embrittlement in steel [34][35][36]. The constitutive behaviour of the bulk material is assumed to be described by an isotropic von Mises plasticity model and the flow stress is assumed to be independent of the hydrogen content. Following Sofronis et al. [37], the hardening function is taken to be a power law, i.e. σ y = σ 0 1 + ε pl e /ε 0 where σ 0 is the initial yield strength, ε pl e is the equivalent plastic strain, ε 0 = σ 0 /E, E is Young's modulus and n is the hardening exponent. Further, the dislocation density ρ D is assumed to be linearly related to the equivalent plastic strain according to ρ D = ρ D0 + γ ε pl e for ε pl e ≤ 0.5 and ρ D = ρ max D for ε pl e > 0.5, where ρ D0 is the initial dislocation density and γ is a material parameter [21,29]. The cohesive model in Sect. 4.2 is used to model failure. The hydrogen in the bulk material is assumed to reside in the lattice and dislocation sites, i.e.C T =C L +C D , and it is also trapped in the interface. The hydrogen concentration in the dislocation trapsC D increases with increasing plastic strain, i.e. due to the increase of the number of trap sites associated with the dislocation density Eq. (45). Further, the hydrogen trapped in the dislocations and the discrete trap represented by the interface are in equilibrium with the concentration in the lattice following Oriani's theory, according to Eqs. (46) and (9), respectively. Hence, the displacement and concentration or chemical potential can be determined by solving a coupled mechanical-diffusion problem using the Finite Element Method based on the concentration or chemical potential formulation. The concentration formulation is adopted from Barrera et al. [14] and the chemical potential formulation described above is implemented in the FE code Abaqus [15]. The isotropic von Mises plasticity model is implemented using a UMAT subroutine and the diffusion equations are implemented using the UMATHT subroutine (the details of the implementations are provided in Appendices A and B). The material parameters used in the analysis are summarised in Table 1.
The computational domain of the deep double-edged notched specimen shown in Fig. 3 is analysed under plane strain conditions. (Only a quarter the specimen is analysed due to symmetry.) The dimensions are taken as W = 10 mm, H = 12.5 mm, B = 1 mm, r = 5 mm, D = 5 mm and ϕ = 30 • . The 4-node bilinear and 8-node quadratic coupled temperature-displacement plane strain elements (CPE4RT and CPE8RT, respectively) are used in the different dis-cretisations. In particular, the CPE4RT element is used in the failure simulations and the CPE8RT element is used to compare the two formulations. It should be mentioned that quadratic elements are necessary to compute the stress gradient in the concentration formulation (a detailed discussion is presented in Barrera et al. [14]). A 4-node two-dimensional linear interface element was implemented in Abaqus using the user-defined subroutine UEL. In the comparison between the formulations, no damage is considered and the Finite Element model comprised of the bulk elements only.
In the failure analysis, the Finite Element model is divided into two regions, in which the bulk and interface elements are defined. The interface elements are inserted along the predefined failure path, g W ≥ x ≥ D + r and y = 0, and bulk elements are defined elsewhere. The interface elements are modelled with zero initial thickness (i.e. the top and bottom face nodes coincide) and their top faces are attached to the bulk elements. In failure simulations, the mesh has 5682 elements, of which 5469 are bulk elements and 213 are interface elements. The total number of elements is reduced to 5469 elements in the simulation that are used to compare the formulations. This number of elements is found to be necessary to obtain converged solutions in both cases.

Comparison between the concentration and chemical potential formulations
Firstly, the concentration and chemical potential based formulations are compared. We assume that the hydrogen concentration in the lattice is initially uniform throughout the plate, i.e. the initial lattice concentration isC L0 = 4.3 mol/m 3 which corresponds to an occupancy of θ L0 = 5×10 −6 . All boundaries are insulated. The mechanical loading is applied in the form of a prescribed displacement in the y-direction, u y , at the upper surface (i.e. y = H ) as shown in Fig. 3. The displacement is ramped up monotonically from 0 to 0.025 mm which corresponds to an overall engineering strain of ε yy ∈ [0 − 0.002]. Under these loading conditions only elastic deformation is expected. The process time t p is taken to be equal to the characteristic diffusion time, τ = l 2 c /D L where l c is the characteristic length scale, taken to be equal to the notch radius (1 mm). Hence, transient diffusion is considered.  Figure 4a, b illustrate the distribution of the hydrogen concentration in the latticeC L normalised by the initial con-centrationC L0 and the chemical potential μ normalised by the initial chemical potential μ ini = 30.26 J/mol at strain level ε yy = 0.002. The Fig. 4i, ii show that the two formulations agree. In the elastic regime, the concentration of hydrogen in the lattice is controlled by the hydrostatic stress distribution. Thus, the maximum concentration is at the notch root where the hydrostatic stress is a maximum. Far from the notch, the hydrostatic stress is much lower, which leads to a decrease in the hydrogen concentration locally to satisfy conservation of volume of hydrogen within the specimen.

Failure of the double-edged notched specimen
In this section, we investigate the use of the interface model and the chemical potential formulation to simulate failure of the deep double-edged notched specimen. Similarly, we assume initially uniform hydrogen concentration in the lattice throughout the plate and the mechanical loading is applied in the form of a prescribed displacement in the y-direction at the upper surface. The initial hydrogen concentration in the lattice is taken to beC L0 = 4.3 mol/m 3 and zero flux at the boundaries is prescribed. The displacement at the upper surface is ramped up monotonically from zero to complete failure. The process time is taken to be t p = τ which yields a loading rate of 0.25 mm/s. The ratio between the interface cohesive strength and yield strength is chosen to be σ c /σ 0 = 3.0 (i.e. ξ int = 0) and the separation at failure is taken to be δ f n = 0.01 mm. It should be mentioned that these parameters are suitable to describe steel at the given concentration level [21]. Further, the softening of the interface (i.e. the reduction of cohesive strength) due to hydrogen is not investigated in this example. The interface stiffness is taken to be K n = 10 6 MPa/mm which is found to be necessary to prevent artificial compliance [38,39]. The interface binding  Figure 5 shows the macroscopic stress-strain responses in the y-direction for damaged and undamaged cases, i.e. the relationship between ε yy and σ yy as indicated in Fig. 3. The behaviour illustrates a typical elastic-plastic response for a notched specimen in which the onset of macroscopic yielding usually occurs at a lower stress level than the yield strength of the material due to the stress concentration, i.e. σ yy /σ 0 ≈ 0.6. Damage initiation takes place at ε yy ≈ 0.012 and final failure occurs at ε yy ≈ 0.0135. In what follows, we investigate four points in the loading history, these points are associated with: A, elastic deformation; B, damage initiation; C, damage progression and D, complete failure of the interface, see Fig. 5. Figure 6i-iv illustrate the distribution of the hydrogen concentration in the latticeC L normalised by the initial con-centrationC L0 at points A, B, C and D of Fig. 5, respectively. At point A, prior to any plastic deformation, hydrogen is concentrated at the notch root in the region of maximum hydrostatic stress. After the onset of plastic deformation and before damage initiation, the maximum hydrostatic stress moves to a position below the notch root. Consequently, the maximum concentration shifts to this location, i.e. the maximum value ofC L is at x/W ≈ 0.65 at point B. After damage initiation, the stress relaxes along the interface which results in diffusion of hydrogen from the interface. At complete failure, the specimen is only loaded by the residual stresses that are caused by the plastic deformation. Hence, the hydrogen is redistributed according to this residual stress field.
In order to investigate the interface, we consider the interface and the bulk material adjacent to the interface, i.e. along W ≥ x ≥ D+r and y = 0. Figure 7i, ii show the distribution of the hydrogen concentration in the latticeC Tr normalised by the initial concentrationC Tr0 and the damage variable d, respectively, along the interface at different instants indicated by points A, B, C and D as shown in Fig. 5. At point A, prior to damage initiation, the damage variable is zero along the interface. The increase of deformation causes damage initiation at a distance from the notch root, i.e. x/W ≈ 0.65, where the stress normal to the interface is maximum. Consequently, the damage spreads towards the notch root and the middle of the specimen, as indicated by the damage distribution at points B and C. At failure, the interface is completely Fig. 6 The lattice hydrogen concentration distribution at: i the onset of the plastic deformation, ii damage initiation, iii damage progression and iv complete failure. These instants are associated with points A, B, C and D in Fig. 5, respectively damaged such that the damage parameter is 1.0 along the entire interface. Initially, hydrogen is trapped at the interface and its concentration isC Tr /C Tr0 which is solely determined by the interface's binding energy as the material is initially stress free. Since V int M has been taken to be zero and there is no jacking due to the hydrogen that diffuses to the interface, the concentration remains uniform along the interface as the stress is increased. It's magnitude gradually decreases, however, to supply some of the hydrogen that flows into the regions of high hydrostatic stress as the applied load is increased. Thus the concentration at A is less than that under zero stress. As damage develops the stress in the matrix decreases, hydrogen diffuses away from regions where the hydrostatic stress was high and some of this flows into the interface, leading to an increase inC Tr . After more significant damage development, the traction decreases, which causes additional trapping of hydrogen in the interface, i.e. points C and D. Figure 8i-iv show the distribution of the hydrostatic stress σ h normalised by the yield strength σ 0 , equivalent plastic strain ε pl e , hydrogen concentration in the latticeC L and hydrogen trapped in the dislocation sitesC D normalised by the initial concentration in the latticeC L0 in the bulk material adjacent to the interface at different instants that are indicated by points A, B, C and D as shown in Fig. 5. The distributions of hydrostatic and equivalent stress imply that prior to the onset of plastic deformation, at point A, the equivalent plastic strain is zero and the hydrostatic stress is maximum at the notch root. The plastic deformation develops at the notch root and spreads toward the centre of the specimen.
Consequently, the maximum hydrostatic stress changes location from the notch root. As damage develops, the stresses relax which halts the development of plastic deformation and causes the maximum hydrostatic stress location to move to the partially damaged ligament at the centre of the specimen. After final failure, the hydrostatic stress is determined by the distribution of the residual stresses. The distribution of the hydrogen concentration in the lattice reflects the hydrostatic stress distribution. Further, the hydrogen trapped in dislocations is in equilibrium with the hydrogen in the lattice which is demonstrated by the similarities of the distributions. The number of dislocation sites in the vicinity of the notch increases as suggested by the plastic deformation distribution. However, a small change in the trapped hydrogen concentration is observed.

Analysis of hydrogen enhanced decohesion (HEDE) in a dissimilar metal weld system
In this example, we investigate hydrogen enhanced decohesion (HEDE) in an AISI8630/IN625 dissimilar weld system. (A comprehensive investigations of this system is provided in Barrera et al. [33].) In particular, we are interested in studying the decohesion process that takes place around fine M 7 C 3 carbide particles in a 'featureless' zone located on the Nickel side of the interface in the presence of hydrogen. Barrera et al. [33] report analyses of this system using TEM imagery and used a mesh generation scheme that converted these images into a finite element mesh [40]. Their simulations have shown that during deformation microcracks are initially formed at the carbide-matrix interface which then propagate along the interface. Thereafter, several microcracks connect together and form macrocracks through the localisation of plastic flow in a region adjacent to the region where (1) hydrogen content is high and (2) the carbide/matrix interface has debonded.
The featureless zone is located next to the fusion line on the Nickel side of the weld and is typically about 20 μm wide. It is rich in M 7 C 3 carbides that occupy a volume fraction of about 15% and the average length of their major axis is of the order of 40 nm. In this analysis, we adopt the representative microstructure of the featureless zone from Barrera et al. [33] (see Fig. 9i) in which several precipitates are distributed in the Nickel matrix. Further, we focus on analysing a single precipitate and consider the volume element in Fig. 9ii assuming plane strain conditions. The objective of this analysis is to initially reproduce the results obtained by Barrera et al. [33] and then explore some of the features of this problem. Bar- Fig. 9 The representative volume of the 'featureless' region: i cluster of M 7 C 3 carbides in a 623 Nickel alloy; and ii single precipitate RVE. The single precipitate vertices are numbered from 1 to 10 and a local coordinate ξ along the interface with origin at vertex 1 is introduced rera et al. [33] did not undertake a fully coupled mechanical diffusion analysis. They assumed a simple steady state distribution of hydrogen, and expressed the material parameters in terms of this distribution. The current analysis consists of a fuller physical description of the problem, and captures the detailed evolution of the hydrogen concentration in the material as plastic deformation and damage develop. In particular, we are interested in investigating: (1) hydrogen enhanced decohesion at the interface and how this influences the hydrogen distribution; (2) the effect of swelling of the matrix and precipitate on the decohesion process and (3) the effect of hydrogen charging on deformation in the matrix. To analyse this problem, as noted above, we assume a full coupling between the mechanical and diffusion responses. The matrix material (i.e. 625 Nickel) is assumed to be described by an isotropic von Mises plasticity model with power law hardening as in Sect. 5.1. The flow stress in the bulk material is assumed to be independent of the hydrogen content. The relevant matrix material parameters are: Young's modulus E m , Poisson's ratio ν m , initial yield strength σ 0 , hardening exponent n, initial dislocation density ρ D0 and dislocation density parameter γ . The precipitate ( M 7 C 3 -carbide) is taken to be elastic with Young's modulus E p and Poisson's ratio ν p . The bilinear interface/cohesive model described in Sect. 4.2 is used to model the interface. The hydrogen in the bulk and precipitate materials is assume to reside in the lattice, dislocations and fixed traps, i.e.C T =C L +C D +C F . Hydrogen is also trapped in the interface, with concentrationC Tr . The Nickel alloy parameters are given in Table 1 and the precipitate parameters are taken to be E p = 2 E m and ν p = ν m . The reference chemical potential for the matrix and the precipitate are taken to be μ m 0 = 0 and μ p 0 = μ m 0 + W p = W p , respectively, where W p is the binding energy of the precipitate. The molar volume of the matrix and precipitate are taken to be V p M = V M = 2 × 10 −6 m 3 /mol. The effects of the trapping energies will be discussed later. The interface parameters are W int , V int M = 2 × 10 −6 m 3 /mol and h int = 8 × 10 −6 mol 2 /(m 2 · s · J). The chemical potential formulation and its implementations in the FE code Abaqus [15], described in Sect. 3 and Appendices A and B, is used. The Finite Element model of the single precipitate problem is shown in Fig. 10i, ii. The model is divided into three regions, namely; the matrix, precipitate and interface. Interface elements are inserted along the interface and the bulk elements are defined in the precipitate and the matrix regions. The interface elements are modelled with zero initial thickness (i.e. the top and bottom face nodes coincide) and their top and bottom faces are attached to the bulk elements in the matrix and precipitate regions, respectively, as shown in Fig. 10i. Further, a uniform refined element region is created adjacent to the interface for controlling the interface element length. The 4-node bilinear coupled temperaturedisplacement plane strain elements (CPE4R) are used in the discretisation for plane strain conditions, respectively. The 4-node two-dimensional linear cohesive user element in Appendices A and B is used to discretise the interface. The mesh has 19408 elements, of which 10662 and 8095 are bulk elements in the matrix and precipitate, respectively, and 651 are interface elements. The number of elements is found to be necessary to obtain converged solutions.
The hydrogen concentration in the lattice is initially assumed to be uniform throughout the body, i.e. the initial lattice concentration isC L0 = 4.3 mol/m 3 which corresponds to an occupancy θ 0 = 5.0 × 10 −6 . All boundaries are insulated. The mechanical loading is applied in the form of a prescribed displacement in the y-direction, u y , at the upper surface as shown in Fig. 10i. The displacement is ramped up monotonically from zero to complete failure. The loading time t p is taken to be equal to the characteristic time for diffusion τ = l 2 c /D L = 31.1 ns where the characteristic length scale is taken to be equal to the precipitate spacing l c ≈ 50 nm. In the following sections we will present the results of the different investigations.

Analysis of hydrogen enhanced decohesion at the interface
In order to study the hydrogen enhanced decohesion (HEDE) process along the interface, the interface strength is explicitly reduced. Thus, the ratio between the interface strength and yield strength is chosen to be in the range σ c /σ 0 ∈ [1.0 − 3.0], i.e. ξ int = 0. It should be mentioned that, in the absence of hydrogen, this ratio is typically greater than 1 which allows plastic deformation to spread in the matrix material prior to the failure of the interface. The separation at failure and interface stiffness are taken to be δ f n = 0.5 μm and K n = 10 6 MPa/mm, respectively. The reference chemical potentials for the matrix and precipitates are taken to be μ m 0 = μ p 0 = 0, i.e . W p = 0. Figure 11 shows the macroscopic stress-strain responses in the y-direction for different values of σ c /σ 0 , i.e. the relation between Σ yy and E yy as indicated in Fig. 10i. The macroscopic stress-strain response shows that the material is initially elastic and yields at a stress level that is approximately equal to the macroscopic yield strength. After yielding, the material shows strain hardening behaviour which is followed by a significant drop in the stress due to damage development at the interface. Further, reducing the cohesive strength causes a significant decrease in the ductility. More specifically, at higher values of σ c /σ 0 , the strain to failure drops by 50% as the result of reducing the cohesive strength by 17%. Further, the reduction in ductility is smaller at lower values of σ c /σ 0 , i.e. ≈ 20%. Additionally, stress softening is observed which is attributed to damage development along the interface.
To investigate the details of the deformation and diffusion processes, the distribution of lattice occupancyC L normalised by the initial concentrationC L0 and equivalent plastic strain ε pl e for the case of σ c /σ 0 = 2.0 at different instants are plotted in Fig. 12a, b, respectively. In particular, the results are explored at four points A, B, C and D in the loading history as illustrated in Fig. 11, represent- . At point A, prior to any plastic deformation, hydrogen is concentrated at some of the vertices of the precipitate, i.e. 1, 5, 6, 7, 9 and 10, where the hydrostatic stress attains its maximum values. As deformation proceeds, plastic deformation develops at precipitate vertices and radiates toward the upper and lower surfaces. Further, the maximum hydrogen concentration remains at the vertices, where the hydrostatic stress is greatest. The damage at the interface initiates at the uppermost facet of the precipitate, i.e. facet 6 − 7. This facet makes the largest angle with the loading direction, i.e. y-direction, which results in a maximum normal traction distribution along this segment. After damage initiation, the traction decreases along this part of the interface which causes stress relaxation in the matrix and in the precipitate adjacent to the interface. As a result, hydrogen diffuses away from the interface to regions where the hydrostatic stress is larger. Further, the damage continues to increase along this segment leading to a complete separation. The newly created crack propagates along interface segments 5 − 6 and 7 − 8. The intensive plastic deformation at the crack tips stops further growth of the crack. Further, hydrogen diffuses to the crack tips where now the hydrostatic stress is greatest. Simultaneously, the traction at the lowermost facet of the interface, i.e. facet 10 − 1, increases and approaches the cohesive strength causing damage to initiate. It should be mentioned that this segment makes the second largest angle with the loading direction after the uppermost facet. Eventually, plastic deformation spreads widely in the element forming several bands wherein the equivalent plastic strain reaches a maximum value of 50%. At this strain level, plastic localisation is expected to take place between adjacent precipitates leading to complete failure.
To analyse the interface, the distribution of the damage variable d along the interface at the instants A, B, C and D are plotted in Fig. 13. The result shows that at point A, prior to damage initiation, the damage variable is zero along the interface. The damage initiates at segment 6 − 7 and propagate to segments 5 − 6 and 7 − 8. Further, damage initiates along segments 1 − 2 and 10 − 1 at a later stage of deformation. It should be mentioned that the hydrogen trapped in the interface remains constant during the deformation. After complete failure of the interface, i.e. point D, a small increase is observed along segments 5 − 6 and 6 − 7.
The result of this simulation is consistent with the findings reported in Barrera et al. [33], although the current study provides more details of the way in which the diffusional, deformation and interfacial failure processes are coupled, which provides more insight into the sequence of events leading to failure.

The effect of the swelling in the precipitate and matrix
To analyse the effect of the differential swelling of the matrix and precipitate on the failure of this material, the ratio between the precipitate and matrix reference chemical potentials, μ p 0 /μ m 0 , is varied over the range 0.5 − 2.0. In particular, the reference energy of the matrix is taken to be μ m 0 = −10 kJ/mol and the reference energy of the precipitate, μ p 0 , is then selected to give the required ratio. Hence, the simulations cover the cases where the matrix swells more than the precipitate and the vice versa. The ratio between the interface and yield strength is chosen to be σ c /σ 0 = 2.0 and the separation at failure and interface stiffness are kept as δ f n = 0.5 μm and K n = 10 6 MPa/mm, respectively. Figure 14 shows the macroscopic stress-strain responses in the y-direction for different values of μ p 0 /μ m 0 . The result implies that both critical stress (i.e. the stress at which damage is initiated) and ductility depend significantly on the swelling. More specifically, at lower values (μ p 0 /μ m 0 ≤ 1.0), smaller changes are observed. At larger values (μ p 0 /μ m 0 > 1.0), the swelling causes an increase in both the critical stress and the ductility. To better understand this increase in the critical stress and ductility, the initial distribution of the normal traction T n normalised by the cohesive strength σ c along the interface is plotted in Fig. 15 for the different values of the ratio μ Hence, when the element is loaded, the traction along the interface increases from the initially negative values which results in an increase in both the critical stress and ductility. This behaviour can be explained by investigating the relative volumetric expansion between the matrix and precipitate. Thus, at the lower values of μ p 0 /μ m 0 , the hydrogen concentration in the matrix is larger than the precipitate which leads to a larger volumetric expansion in the matrix and a positive relative volumetric expansion, see Fig. 15. On the other hand, at the higher values μ p 0 /μ m 0 , the hydrogen concentration in the precipitate is larger than the matrix which leads to a larger volumetric expansion in the precipitate and negative relative volumetric expansion, see Fig. 16. Hence, the interface is subjected to a tensile loading at lower values of μ p 0 /μ m 0 due to the positive relative volumetric expansion and compressive loading at higher values due to the negative relative volumetric expansion.

The effect of hydrogen charging on deformation
Another aspect of the problem is to investigate the nature of the deformation induced by hydrogen charging. The deformation in the RVE is mainly induced by volumetric expansion of the matrix and the precipitate. Hence, the ref-   Figure 17 shows the distribution of the equivalent plastic strain ε pl e for different values of the initial lattice concen-trationC L0 . The result shows that the plastic deformation initiates at the precipitate vertices where the stress is a maximum and at a concentration ofC L0 = 4.30 mol/m 3 . Further, increasing the initial concentration increases the plastic deformation. Thus, for the given values of μ m 0 , μ p 0 , σ 0 and E, there is a hydrogen concentration level at which plastic deformation may initiate. Moreover, the size of the plastic zone and the extent of plastic deformation additionally depends on the hardening exponent n.

Concluding remarks
In this paper, a finite element formulation for solving coupled mechanical/diffusion problems is proposed. The purpose is to investigate hydrogen diffusion in metals and its impact on their mechanical behaviour (i.e. hydrogen embrittlement). In particular, the hydrogen enhanced decohesion and hydrogen enhanced local plasticity mechanisms are incorporated. The hydrogen atoms are assumed to reside either at NILS or at trapping sites such as dislocations and/or fixed traps (e.g. precipitates, interfaces, etc.) that are taken to be in equilibrium following Oriani's theory. The HELP mechanism is modelled using a continuum plasticity framework, e.g. isotropic von Mises model, in terms of a decrease in the flow stress   [13] and Barrera et al. [14]. The HEDE mechanism is modelled using a cohesive zone modelling approach such that an interface is treated as a discrete trap and its cohesive strength decreases with the increase of hydrogen concentration at the interface. The proposed formulation adopts a standard Galerkin method in the discretisation of both the diffusion mass conservation and mechanical equilibrium equations. In particular, a displacement based finite element formulation with chemical potential as an additional degree of freedom is employed. Hence, in comparison with the commonly used concentration based formulation, we have substituted the concentration by the chemical potential in the discretisation of the diffusion equation. It follows that the diffusion equation can be expressed fundamentally in terms of the gradient in chemical potential, which reduces the continuity requirements on the shape functions to zero degree (i.e. C 0 ). Therefore, a linear shape function can be used to interpolate the chemical potential. In the concentration based formulation, the continuity requirements is of the first degree C 1 , i.e. at least quadratic polynomial. Thus, using the chemical potential offers a lower continuity requirements in comparison with the concentration based formulation. We introduced a consistent interface element formulation that can be achieved due to the continuity of the chemical potential across the interface -concentration can be discontinuous at an interface and this can lead to numerical problems in formulations where concentration is employed as a state variable. Moreover, the coding of the FE equations is more straightforward in the proposed model.
In order to investigate the efficiency and accuracy of the proposed formulation, various example problems are studied. Initially, we have investigated a 2D fully coupled elastoplastic diffusion problem of a deep notched specimen in the absence of damage. The solutions using both concentration and chemical potential based formulations are compared and found to be identical. We have then studied the failure of the specimen using the chemical potential formulation. We have explored the effect of deformation and damage initiation and propagation on the hydrogen of the bulk and the cohesive zone. The results show that the concentration of hydrogen in the lattice is controlled by the hydrostatic stress distribution which takes its maximum value at the notch root in the elastic deformation regime and at a distance from the notch root in the elastic-plastic deformation regime. After the onset of damage, the stress relaxes along the interface which results in diffusion of hydrogen away from the interface and the maximum concentration of hydrogen in the lattice takes place at the tip of the region of extensive damage. The second problem is a micromechanical investigation of hydrogen enhanced decohesion (HEDE) in a dissimilar metal weld. More specifically, we consider the failure in the presence of hydrogen of carbide-matrix interfaces around fine precipitates generated adjacent to the interface of an 8630 steel/IN625 nickel alloy dissimilar weld, which has been analysed extensively by Barrera et al [26]. We have investigated the hydrogen induced decohesion at the interface and have shown that the reduction in cohesive strength significantly reduces the ductility and influences the hydrogen distribution. Further, we have studied the effect of differential swelling of the precipitate (arising from differences in the solubility of hydrogen in the matrix and precipitate) on the decohesion process and have illustrated that this differential may delay the onset of damage and consequently increase the ductility. Finally, we have shown that hydrogen charging may introduce localised plastic deformation in the matrix in the absence of any applied mechanical loading, which results in the generation of a residual state of stress in the vicinity of a particle. unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
A Finite element implementations of the coupled elastic-plastic mechanical, diffusion and interface models.
The Finite Element implementations of the continuum plasticity, diffusion and interface models in Sect. 3 are discussed in this Appendix. The purpose is to provide explicit expressions for the residual force vectors and consistent tangents for the Newton-Raphson solution procedure. Hence, we consider the system of equations that are given by the residual vectors and the solution procedure in Sects. 3.4 and 3.5, respectively. The model problem domain is divided into bulk and interface regions. The bulk material is assumed to exhibit elastic-plastic behaviour and the interface behaviour is determined by the cohesive zone model in Sect. 4.2. The hydrogen transport model for the bulk and interface are described in Sects. 4.1 and 4.2, respectively.

A.1 The bulk material
The contributions of the bulk material to the global residual vectors in Eqs. (30) and (34) are defined as where node I belongs to Ω and Γ (i.e. I ∈ Ω, Γ ) and the volume Ω and the boundary Γ define the bulk region. The lattice concentrationC L is determined as described in Sect. 4.1 and the Cauchy stress σ i j is obtained using an elastic-plastic constitutive model that will be discussed later in this section.
In the solution procedure in Sect. 3.5, the linearisation of the residuals yields the tangent stiffness matrix in Eqs. (37) and (38). Hence, the tangent operators for the residuals in Eqs. (A.1) are obtained as where nodes I and J belongs to Ω and Γ , and D μμ , D μu j , D uμ i j , and D uu ik jl are the tangent moduli that can be determined from the variation of the total hydrogen concentration and the Cauchy stress. Therefore, the derivatives and tangent moduli in Eqs. (A.2) will be determined in the subsequent subsections.

A.1.1 The elastic-plastic model
The elastic-plastic constitutive relation is written in terms of the rate of deformation tensor D i j and Zaremba-Jaumann stress rate of Cauchy stress σ i j . The Zaremba-Jaumann stress rate is used to meet the objectivity requirements, i.e. material invariance for rigid body rotation, and it is given by whereσ i j is the material time derivative of the Cauchy stress, and D i j = sym L i j and W i j = asym L i j are the rate of deformation and spin tensors that are equal to symmetric and antisymmetric parts of the velocity gradient L i j =Ḟ ik F −1 k j , respectively, and F i j = ∂ x i /∂ X j is the deformation gradient. Following Sofronis et al. [37], the total rate of deformation is taken to be the sum of the elastic D el i j , plastic D pl i j and swelling D s i j parts: The elastic strains are assumed to be small and given by where S i jkl = C −1 i jkl is the elastic compliance tensor. The presence of hydrogen may cause volumetric change that is given by the swelling rate of deformation .6) and the material time derivative of ε s isε s = 1 3 V MĊT . It should be mentioned that, for the sake of generality, the swelling is taken to be associated with hydrogen in the lattice and traps, i.e. the total concentration, which is different from the definition used in Sect. 4.1. Consequently, the stress contribution to the chemical potential in Eq. (51) becomes μ σ = −σ h V M ∂C T /∂C L . Further, the swelling deformation is assumed to be small. The plastic part of the rate of deformation tensor is given by the von Mises flow rule as where λ is the a positive scalar function and f is the yield function, i.e. f < 0 implies elastic unloading and f = 0 represents the yielding, that is given by where the yield function f is a function of ε pl e andC T . Therefore, using Eqs. (A.6-A.9), the plastic part of the rate of deformation tensor becomes 10) and the rate tangent expression for the Zaremba-Jaumann derivative of Cauchy stress is where K is the bulk modulus, ∂ f /∂ε pl e = ∂σ y /∂ε pl e = H and ∂ f /∂C T = ∂σ y /∂C T can be determined from a hardening law and the associated hydrogen softening function (see Sect. B.2.1), respectively. The different derivatives in Eq. (A.11) determined from Eqs. (42)(43)(44)(45)(46) and (51) are The change of Cauchy's stress with respect to the chemical potential is determined as (A.14) It follows that the tangent moduli in Eq. (A.2) can be determined as (A.15)

A.1.2 The diffusion model
The variation of the chemical potential is expressed as: where the diffusion moduli are defined by (A.17)

A.2 The interface/cohesive model
The contributions of the interface material to the global residual vectors in Eqs. (30) and (34) are where node I belongs to Γ int (i.e. I ∈ Γ int ), the surface Γ int defines the interface region and T int i = t int i . The total con-centrationC int and the traction T int i are obtained using the cohesive model in Sect. 4.2 and the details of the processes will be discussed later in the subsequent subsections. Similarly, the linearisation of the residuals in Eqs. (A.18) yields whereD μμ ,D μu j ,D uμ i , andD uu i j are the tangent moduli of the interface that can be determined from the variation of the total hydrogen concentration and the traction.

A.2.1 The interface concentration and traction
The chemical potential and the tractions at the interface are defined in Eqs. (60), (61) and (65). Using these equations, a set of equations in the incremental forms of the concentration and tractions, at a given values of Δδ n , Δδ t and μ, can be written as whereC p+1 int is the concentration at the end of the increment, i.e.C where the increment of the jacking separation is Δδ s n . In order to solve these equations, we encapsulate the equations in a vector F = F μ F t F n T and define the unknown vector n T . Using Newton-Raphson algorithm, the following linearised form of these equations is iteratively solved for the traction increment δΔT until a convergence condition is achieved where 0 is a zeros vector (3 × 1) and the Jacobian J is defined as where the derivatives are determined from Eq. (A.20) as

A.2.2 The cohesive zone Jacobian
The variation of the hydrogen concentration at the interface is defined bẏ where the diffusion moduli in Eqs. (A.2) are defined bỹ The variation of the traction at the interface is defined bẏ where the deformation moduli in Eqs. (A.2) are defined bỹ The elements of the tangent moduli can be determined from the variations of the incremental forms in Eqs. (A.20). The variation of F i is given by where dF i /dΔδ j denotes the total derivative of the function F i with respect to Δδ j in the unknown vector Δδ = [Δμ int Δδ t Δδ n ] T , the variations of the separation increments δΔμ int , δΔδ t and δΔδ n can be chosen arbitrarily Thus, the unknown derivatives in Eqs. (A.31) can be written in vector form as Hence, the unknown derivatives can be determined by solving the following set of equations

B The Finite Element implementation of coupled mechanical/diffusion framework in Abaqus
In this appendix, we present the implementations of the hydrogen transport, elastic-plastic and interface/cohesive models in Sect. 4 and "Appendix A".

B.1 The hydrogen diffusion Model
The hydrogen transport equations are implemented using Finite Element Method using Abaqus [15]. We use the similarities with the heat equation as discussed in Oh and Kim [21] and Barrera et al. [14]. In these studies, taking advantage of the similar forms of heat and diffusion equations, the temperature T in the heat equation is replaced by the total concentrationC T which yields the diffusion equation. This formulation requires the evaluation of the gradient of the stress tensor, i.e. substitution of Eqs. (50) or (51) into Eq. (5) yields a term that contains ∇σ h . Hence, the Finite Element interpolation functions and its first derivative should exist and be continuous, i.e. C 1 . In other words, at least quadratic shape functions are needed in order for the gradient of the stress tensor to exist. Moreover, the determination of this gradient is not straight forward. In this study, we propose an alternative formulation in which we consider the chemical potential to be the degree of freedom in the diffusion equation. We summarize the analogy of variables between the heat transfer analysis and hydrogen diffusion analysis within Abaqus in Table B.1. The user subroutine UMATHT is used to implement the hydrogen diffusion model. In coupled thermomechanics analysis, the deformation related variables, i.e. σ h and ε pl e , are passed to UMATHT using other user subroutines. The details of passing these variables are explained later in subsequent subsections. In UMATHT, the terms result from the weak form of the heat equation and its variation should   be provided. In particular, the internal energy U q , flux J q and their derivatives with respect to the temperature and temperature gradient, i.e. ∂U q /∂ T , ∂U q /∂∇T , ∂J q /∂ T , ∂J q /∂∇T , have to be defined. Table B.2 shows the analogy between these quantities for the heat and diffusion problems. The total hydrogen concentration can be expressed in terms of the chemical potential μ and the trap densityN D ε pl e . Further, the variation of the hydrogen flux with respect to the chemical potential μ and its gradient ∇μ is required. Thus, the variation of the hydrogen concentration is given in Eqs. (A. 16) and (A.17), and the variation of the hydrogen flux is It should be noted that the hydrostatic stress σ h and effective plastic strain ε pl e are passed into the UMATHT using other user subroutines. In the next section, the details of this procedure is presented

B.2 The elastic-plastic Model
The elastic-plastic model in "Appendix A" can be implemented in Abaqus code using two different user subroutines that are -UHARD and UMAT. In UHARD, the stress integration and the consistent tangent modulus are determined by Abaqus. Derivatives of the yield stress should be provided. In UMAT, the stress integration and the consistent tangent modulus should be provided. It should be mentioned that UHARD is limited to the case of no swelling.

B.2.1 Implementation of the model as an Abaqus UHARD subroutine
In UHARD, the derivatives of the hardening function that are required to evaluate the stress and its variation should be provided. Thus, the variation of the flow stress in Eq. where the initial yield σ H 0 is taken to be a function of the hydrogen concentraion by σ H 0 = Ψ C L σ 0 and Ψ is a monotonically decreasing function of hydrogen concentration in the lattice. Therefore, the derivatives for the power law are The hydrostatic stress σ h and effective plastic strain ε pl e are passed to UMATHT as solution dependent variables at the start of the increment. We use USDFLD subroutine, which allows different field variables to be defined at a material point as a function of time, to obtain the hydrostatic stress and effective plastic strain, and store solution dependent variables (Barrera et al. [14]). More specifically, in USDFLD, the stress field and effective plastic strain are determined at a material point at the start of the increment using utility subroutine GETVRM. The hydrostatic stress is then evaluated using the normal stress components, i.e. σ h = σ kk /3, and together with effective plastic strain are stored into solution dependent variables.

B.2.2 Implementation of the model as an Abaqus UMAT subroutine
The implementation of the elastic-plastic model using UMAT subroutine requires the determination of stress and the consistent tangent modulus at the end of increment. In particular, the rate tangent required by UMAT takes the form which is the rate tangent for the Jaumann derivative of Kirchhoff stress that is scaled by the Jacobian of the deformation. Hence, the stress and the mechanical and chemical tangent moduli C UMAT i jkl and D μu i , respectively, should be provided in the UMAT subroutine.

B.3 The cohesive/interface zone Model
The interface model in Sect. 4.2 and "Appendix A", is implemented into the Finite Element method using a cohesive Fig. B.1 The two-dimensional linear cohesive element i the physical domain and the definition of the local and global coordinates and ii the parent domain element approach and UEL subroutine. The surface-like cohesive formulation is adopted in this work (e.g. see [41][42][43][44]). In this formulation, a cohesive element consists of two surface elements which coincide in the reference configuration. A 4-node two-dimensional linear element is used in this study following the formulation by Elmukashfi and Cocks [45]. The element upper and lower surface elements are denoted S + and S − , respectively, see Fig. B.1. The constitutive behaviour of the cohesive element is defined in the middle surface,S. The middle surface is defined by the two points 1 and 2 which are the mid-points between nodes 1 − −1 + and 2 − −2 + , respectively. Hence, the middle surface is then defined by the interpolation wherex int i I is the current coordinates vector of the middle surface which is related to the current nodal coordinates of the upper and lower surface elements x + i I and x − i I , respectively, bȳ The standard shape functions, N int I , are defined in the local coordinates, ξ ∈ [−1, 1], as It should be mentioned that H I is chosen such that the chemical potential is continuous across the interface, i.e. Eq. (9) , and u int i = Δu int i in Eq. (14). The internal nodal forces and the consistent tangent stiffness matrices across the cohesive surface are determined from the tractions and chemical potential are detailed in the previous sections. In particular, the finite element formulation is detailed in Sect. 3, the interface model is given in Sect. 4.2 and the determination of the nodal forces and the stiffness matrix are illustrated in "Appendix" A.