A phase field formulation for hydrogen assisted cracking
Emilio Martínez-Pañeda Alireza Golahmar Christian F. Niordson
































mail@empaneda.com





[Cam]Department of Engineering, University of Cambridge, CB2 1PZ Cambridge, UK

[DTU]Department of Mechanical Engineering, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark

[cor1]Corresponding author.



We present a phase field modeling framework for hydrogen assisted cracking. The model builds upon a coupled mechanical and hydrogen diffusion response, driven by chemical potential gradients, and a hydrogen-dependent fracture energy degradation law grounded on first principles calculations. The coupled problem is solved in an implicit time integration scheme, where displacements, phase field order parameter and hydrogen concentration are the primary variables. We show that phase field formulations for fracture are particularly suitable to capture material degradation due to hydrogen. Specifically, we model (i) unstable crack growth in the presence of hydrogen, (ii) failure stress sensitivity to hydrogen content in notched specimens, (iii) cracking thresholds under constant load, (iv) internal hydrogen assisted fracture in cracked specimens, and (v) complex crack paths arising from corrosion pits. Computations reveal a good agreement with experiments, highlighting the predictive capabilities of the present scheme. The work could have important implications for the prediction and prevention of catastrophic failures in corrosive environments. The finite element code developed can be downloaded from www.empaneda.com/codes




Phase field Hydrogen embrittlement Stress-assisted diffusion Finite element analysis Fracture

















§ INTRODUCTION



Hydrogen has been known for over a hundred years to cause catastrophic failure in metallic materials <cit.>. Hydrogen atoms enter the material, migrate through the crystal lattice and induce fracture, with cracking being observed in modern steels at one-tenth of the expected fracture toughness. Hydrogen embrittlement is a pervasive problem in the manufacturing, transport, and energy sectors and could compromise the key role of hydrogen as next-generation energy resource. There is a strong need to understand and model how hydrogen degrades the ductility and the fracture resistance of metals.


Modeling efforts in the hydrogen embrittlement community have been largely restricted to discrete numerical methods. Crack initiation thresholds and growth rates have been computed through dislocation-based methods <cit.>, weakest-link approaches <cit.>, and, eminently, cohesive zone models <cit.>, which can easily incorporate the cohesive strength reduction with hydrogen coverage. These methods are however limited when dealing with the complex conditions of practical applications and many important challenges (mixed mode, interacting cracks, complex distribution of initiation sites, etc.) remain unaddressed.


The phase field method for fracture has emerged as a promising variational approach to overcome the limitations of discontinuity-based methods. The method has received increasing attention since the earlier work by Francfort and Marigo <cit.> and important efforts have been devoted to improve solution schemes <cit.>, discretization strategies <cit.> and implementation aspects <cit.>. The model builds upon Griffith's thermodynamics framework and has not only been successfully applied to model brittle fracture but also to ductile damage <cit.>, hydraulic fracturing <cit.>, composites delamination <cit.>, and stress corrosion cracking <cit.>, to name a few. The phase field method is particularly well suited to capture the fracture energy dependence on hydrogen concentration, and will allow addressing the complex cracking conditions intrinsic to hydrogen embrittlement in industrial components.


In this work, we develop a novel variational framework for hydrogen assisted cracking. Hydrogen transport towards the fracture process zone and subsequent cracking are investigated by means of a coupled mechanical-diffusion-phase field finite element framework. The numerical scheme is grounded on: (i) a phase field description of fracture, (ii) an extension of Fick's law for mass diffusion, and (iii) a quantum-based relation between the hydrogen content and the fracture surface energy. We demonstrate the potential of the proposed modeling framework by means of representative case studies.nThe remainder of the paper is organized as follows. Section <ref> describes the hydrogen-dependent phase field formulation. Section <ref> provides details of the stress-assisted impurity diffusion scheme and its coupling with the phase field model. In Section <ref> model predictions are benchmarked against relevant experimental measurements. Finally, concluding remarks are given in Section <ref>.



§ A PHASE FIELD FORMULATION FOR HYDROGEN EMBRITTLEMENT




 §.§ Phase field approximation of fracture



In a one-dimensional setting, the topology of a sharp crack (see Fig. <ref>a) can be described by an auxiliary field variable ϕ(x) ∈[0,1] with

    ϕ(x) =
        
          1       if  x=0
    
          0       if  x ≠ 0

which is referred to as the crack phase field order parameter, with ϕ=0 and ϕ=1 respectively denoting the intact and fully broken states of the material. The non-smooth crack phase field (<ref>) can be approximated by the exponential function,

    ϕ(x)= e^-|x|ℓ

representing a regularized or diffusive crack topology, as illustrated in Fig. <ref>b. The length-scale parameter ℓ determines the width of the smearing function, approaching the discrete crack topology as ℓ→ 0. 


Consider a discrete internal discontinuity Γ in a solid body Ω (Fig. <ref>a), a regularized crack functional Γ_ℓ can be defined (Fig. <ref>b)

    Γ_ℓ(ϕ)= ∫_Ω( 12ℓϕ^2 + ℓ2|∇ϕ|^2)  dV = ∫_Ωγ_ℓ(ϕ,∇ϕ)  dV

that -converges to the functional of the discrete crack surface for a vanishing length scale parameter ℓ→ 0. Here, γ_ℓ(ϕ,∇ϕ) is the crack surface density function. Consequently, the phase field method for fracture circumvents numerical complications inherent to tracking the evolving discontinuity boundary Γ, enabling to robustly model interactions and branching of cracks of arbitrary topological complexity <cit.>. The fracture energy due to the formation of a crack is then approximated as,

    ∫_Γ G_c ( θ) dS ≈∫_Ω G_c ( θ) ( 12ℓϕ^2 + ℓ2|∇ϕ|^2)  dV

where G_c is the critical Griffith-type energy release rate, which is dependent on the hydrogen coverage θ.



The phase field fracture method resembles traditional continuum damage mechanics models, where the scalar damage field may be interpreted as the phase field ϕ. Accordingly, the loss of stiffness associated with mechanical degradation of the material is characterized as a function of ϕ. This is done by means of the so-called degradation function g ( ϕ), which relates the stored bulk energy per unit volume to the strain energy density of the undamaged solid,

    ψ=g ( ϕ) ψ_0 = [ ( 1 - ϕ)^2 + k ] ψ_0
 
Here, k is a parameter chosen to be as small as possible to keep the system of equations well-conditioned; a value of k=1 × 10^-7 is adopted throughout this work. Integrating over the volume of the body and adding the phase field free energy contribution (<ref>) one reaches the total potential energy functional of the deformation-fracture problem,

    Ψ=∫_Ω{[(1-ϕ)^2+k] ψ_0 + G_c ( θ) (12ℓϕ^2 + ℓ2|∇ϕ|^2) } dV

The variation of (<ref>) with respect to δϕ renders the weak form of the phase field contribution,

    ∫_Ω[ -2(1-ϕ)δϕ ψ_0 +
            G_c ( θ) ( 1ℓϕδϕ
            + ℓ∇ϕ·∇δϕ) ]   dV = 0

and the phase field equilibrium equation in Ω,

    G_c ( θ) (1ℓϕ - ℓΔϕ) - 2(1-ϕ)  ψ_0 = 0

follows immediately upon making use of the product rule and Gauss' divergence theorem.



 §.§ Hydrogen-dependent surface energy degradation



Hydrogen embrittles metallic materials by lowering the bond energy between metal atoms, which translates into a reduction of the fracture resistance. A number of authors have employed Density Functional Theory (DFT) to investigate the decohesion of fracture surfaces with varying hydrogen coverage (see, e.g., <cit.> and references therein). For example, Alvaro et al. <cit.> computed the change in ideal fracture energy in the presence of hydrogen atoms at Σ3 and Σ5 grain boundaries in nickel. Their results, in terms of normalized surface energy γ( θ) / γ( 0 ) versus hydrogen coverage θ, are shown in Fig. <ref>, along with a linear fit to the data.



Equivalently, one can define the critical energy release rate dependence on the hydrogen coverage as,

    G_c ( θ)/G_c ( 0)=1 - χθ


where G_c( 0) is the critical energy release rate in the absence of hydrogen and χ is the damage coefficient that weights the hydrogen-lowering of the fracture energy. χ can be estimated for other materials by fitting DFT data from the literature. The damage coefficient for iron and aluminum is obtained from the work by Jiang and Carter <cit.> - see Table <ref>.



The effect of hydrogen on the fracture resistance can be illustrated by examining the analytical homogeneous solution of a one-dimensional quasi-static problem. Under these circumstances, the Cauchy stress is given by σ=g(ϕ)E ϵ, where E denotes Young's modulus and ϵ the strain. Given a strain energy density ψ_0=E ϵ^2/2, one can readily obtain the homogeneous phase field from the strong form (<ref>),

    ϕ=E ϵ^2 ℓ/G_c ( θ)+E ϵ^2 ℓ


and substituting into the constitutive equation renders the characteristic relation between the homogeneous strain and the homogeneous stress,

    σ=( G_c ( θ)/G_c( θ) + E ε^2 ℓ)^2 E ε


The homogeneous solution for the stress reaches a maximum at a critical stress quantity,

    σ_c= √(27 E G_c ( θ)/256 ℓ)


with the strain counterpart given by,

    ϵ_c = √(G_c ( θ)/3 ℓ E)


Thus, ℓ may be interpreted as a material parameter, which governs the magnitude of the critical stress at which damage initiates. This necessarily implies that the value of G_c is intrinsically related to the choice of ℓ, being equal to the Griffith's critical energy release rate for ℓ→ 0. From (<ref>) an analogy with cohesive zone formulations can be established. With this in mind, it is interesting to note that the characteristic element size h in cohesive zone analyses is typically chosen to be

    h < π/160E G_c/σ_c^2


so as to resolve the fracture process zone with at least 20 elements <cit.>. This implies that h should be at least 5.4 times smaller than ℓ in the phase field problem. The constitutive relation of the homogeneous phase field problem is shown in Fig. <ref> as a function of hydrogen coverage.



As shown in Fig. <ref> for the case of iron, hydrogen significantly degrades the strength and the fracture resistance. The present modeling framework quantitatively accounts for the sensitivity of the surface energy to hydrogen coverage, as observed in first principles calculations; other mechanisms of hydrogen damage can also be incorporated.


Finally, we use the Langmuir-McLean isotherm to compute the surface coverage from the bulk hydrogen concentration C,

    θ=C/C+exp( -Δ g_b^0/RT)


with C given in units of impurity mole fraction. Here, R is the universal gas constant, T the temperature and Δ g_b^0 is the Gibbs free energy difference between the decohering interface and the surrounding material. As in <cit.>, a value of 30 kJ/mol is assigned to Δ g_b^0 based on the spectrum of experimental data available for the trapping energy at grain boundaries. Thus, the present formulation accounts for the effect of microstructural traps on cracking and can incorporate the influence on mass transport through an effective diffusion coefficient.



§ STRESS-ASSISTED HYDROGEN DIFFUSION COUPLED WITH PHASE FIELD FRACTURE



The formulation presented in this section refers to the response of a solid body Ω with external surface ∂Ω of outward normal n - see Fig. <ref>. With respect to the displacement field u, the outer surface of the body is decomposed into a part ∂Ω_u, where the displacement is prescribed by Dirichlet-type boundary conditions, and a part ∂Ω_h, where the traction h is prescribed by Neumann-type boundary conditions (see Fig. <ref>a). A body force field per unit volume 𝐛 can also be prescribed. The fracture phase field ϕ is driven by the displacement field of the solid. A Dirichlet-type boundary condition can be prescribed at Γ, a given crack surface inside the solid (see Fig. <ref>b). Additionally, a phase field fracture microtraction f can be prescribed on ∂Ω_f. With respect to the hydrogen concentration C, the external surface consists of two parts (see Fig. <ref>c): ∂Ω_q, where the hydrogen flux J is known (Neumann-type boundary conditions), and ∂Ω_C, where the hydrogen concentration is prescribed (Dirichlet-type boundary conditions). Accordingly, a concentration flux entering the body across ∂Ω_q can be defined as q=J·n.





 §.§ Governing balance equations



The principle of virtual work, in the absence of body forces, reads

    ∫_Ω( σ : ∇δu + ω δϕ + ζ·∇δϕ - dC/dt δμ
       + J·∇δμ)  dV = ∫_∂Ω( h·δu + f  δϕ + q  δμ)  dS

where μ is the chemical potential that drives the flux of hydrogen, σ is the symmetric Cauchy tensor, and ω and ζ are microstress quantities work conjugate to the phase field ϕ and the phase field gradient ∇ϕ, respectively. Applying the product rule and Gauss' divergence theorem to the internal virtual work - left hand side of Eq. (<ref>) - renders,

    δ𝒲_i =    ∫_∂Ω( σ n·δu + ζ·n δϕ + J·n δμ)  dS - ∫_Ω( ( ∇·σ) ·δu
       + ( ∇·ζ - ω) δϕ + ( ∇·J + dC/dt) δμ)  dV

Since the second integral on the right-hand side of Eq. (<ref>) should vanish for arbitrary variations, a set of three equilibrium equations can be obtained,

    ∇·σ   = 0onΩ
    ∇·ζ - ω   = 0  onΩ
    dC/dt + ∇·J   = 0  onΩ

The first integral on the right hand side of (<ref>) may be identified as part of the external virtual work. Thus, by accounting for the right hand side of (<ref>), the corresponding set of boundary conditions is obtained,

    h = σ n   on∂Ω_h
    
            f= ζ·n   on∂Ω_f
    
            q= J·n   on∂Ω_q



 §.§ Energy imbalance



The first two laws of thermodynamics within a dynamical process of specific internal energy ℰ and specific entropy Λ read

    d/dt∫_Ωℰ dV = W_e(Ω) - ∫_∂ΩQ·n dS + ∫_Ω Q  dV


    d/dt∫_ΩΛ dV ≥∫_∂ΩQ/Θ·n dS + ∫_ΩQ/Θ dV


where W_e is the external work, Q is the heat flux, Q is the heat absorption and Θ the absolute temperature. Within an isothermal process (Θ=Θ_0) the thermodynamic laws can be readily combined by considering the free energy Ψ=ℰ - Θ_0 Λ, such that

    d/dt∫_ΩΨ dV ≤ W_e(Ω)


which is generally referred to as Clausius-Duhem inequality or the principle of dissipation <cit.>. Given (<ref>), the local free-energy inequality takes the form

    d/dt∫_ΩΨ dV ≤∫_∂Ω_hh·u̇ dS + ∫_∂Ω_f f  ϕ̇ dS + ∫_∂Ω_q q  μ̇ dS


Employing the divergence theorem and considering (<ref>)-(<ref>) one reaches the local dissipation postulate,

    σ : ∇u̇ + ωϕ̇ + ζ·∇ϕ̇ + μĊ + J·∇μ  - Ψ̇≤ 0




 §.§ Constitutive theory



The constitutive prescriptions are constructed, in a thermodynamically consistent manner, by ensuring that the dissipation condition (<ref>) is satisfied. The free energy is defined as,

    Ψ( ϵ,  ϕ,  ∇ϕ,   C ) =   ( 1- ϕ)^2 ψ_0 ( ϵ) - K V̅_H ( C - C^0 ) tr ϵ_Ψ^b
        +  G_c ( θ) ( 1/2 ℓϕ^2 + ℓ/2 |∇ϕ|^2 )_Ψ^s
        + μ^0 C + R T N ( θ_L lnθ_L + ( 1 - θ_L ) ln( 1 - θ_L ) )_Ψ^c

where Ψ^b ( ϵ,  ϕ,   C ), Ψ^s (ϕ,  ∇ϕ,    C  ) and Ψ^c ( C ) respectively denote the chemo-elastic energy stored in the bulk <cit.>, the crack surface energy, and the chemical free energy. Here ϵ is the strain tensor, K is the bulk modulus, V̅_H the partial molar volume of hydrogen in solid solution, C^0 the reference lattice hydrogen concentration, μ^0 the reference chemical potential, θ_L the occupancy of lattice sites and N the number of lattice sites. We subsequently develop a constitutive theory from the free energy that couples deformation, fracture and hydrogen diffusion.



  §.§.§ Mechanical deformation


We define ψ_0 as the elastic strain energy density for the undamaged solid,

    ψ_0(ϵ)=12ϵ^T : 𝐂_0 : ϵ

where C_0 is the linear elastic stiffness matrix. A linear elastic description of the material response is considered appropriate to characterize the brittle behavior of metals in the presence of hydrogen. Implicitly multi-scale plasticity formulations (e.g., strain gradient plasticity) have shown that crack tip stresses obtained via linear elasticity are more accurate than those computed by means of conventional plasticity theories <cit.>. We assume small strains and define the strain tensor as,

    ϵ=12[∇𝐮^T+∇𝐮]


Accordingly, the Cauchy stress tensor can be derived from (<ref>) as,

    σ=∂Ψ/∂ϵ = g ( ϕ)  ( C_0 : ϵ) - K V̅_H ( C - C^0 ) I

with the second term having a negligible effect in hydrogen embrittlement phenomena <cit.> and being subsequently disregarded in the analysis. As shown in (<ref>), and discussed in Section <ref>, the phase field model characterizes the loss in stiffness associated with bond breakage by means of the degradation function g ( ϕ).



  §.§.§ Phase field fracture


The scalar microstress ω, work conjugate to the phase field ϕ,

    ω=∂Ψ/∂ϕ = - 2 (1-ϕ)  ψ_0  ( ϵ) + G_c ( θ) 1/ℓϕ

and the vector microstress ζ, work conjugate to the phase field gradient ∇ϕ,

    ζ=∂Ψ/∂∇ϕ = G_c ( θ) ℓ∇ϕ

can be readily derived from the free energy (<ref>). Inserting these constitutive relations in the local balance (<ref>b) renders (<ref>) if the concentration gradient along the small region where ∇ϕ≠ 0 is neglected. On the one hand, the deformation of the solid drives the fracture phase field through the strain energy density ψ_0. On the other hand, the coupling between the phase field fracture and hydrogen transport is introduced by defining a hydrogen concentration-dependent crack surface energy that models decohesion enhancement by hydrogen. We provide this connection in a multi-scale fashion by presenting a quantum-mechanically informed phase field degradation law (see Section <ref>).



  §.§.§ Transport of diluted species


The gradient of the chemical potential is the driving force for hydrogen diffusion. One can readily derive the chemical potential from the free energy (<ref>) by considering the relation between the occupancy and the number of sites, θ_L=C/N,

    μ=∂Ψ/∂ C =μ^0 + RT lnθ_L/1 - θ_L - V̅_H σ_H + G_c' ( θ) θ' (1/2 ℓϕ^2 + ℓ/2|∇ϕ|^2 )

and the constitutive framework is completed by the definition of the mass flux, which follows a linear Onsager relationship,

    J= - D C/R T∇μ

where D is the diffusion coefficient. Thus, hydrogen atoms migrate from regions of high chemical potential to regions of low chemical potential. The coupling with the deformation of the solid is given by the stress-dependent part of μ. As evident from (<ref>), hydrostatic tensile stresses increase hydrogen solubility in the lattice (lattice dilatation) by lowering the chemical potential. Correspondingly, the last term in (<ref>) enhances diffusion from cracked regions to pristine regions. However, the appropriate chemical boundary conditions in a propagating crack require careful consideration. First, one should note that the discussion is only relevant to experiments where the specimen is charged in-situ, as opposed to tests where hydrogen has been dissolved within the material prior to loading. Secondly, one would typically expect that the environment will promptly occupy the space created with crack advance. Consequently, a constant hydrogen concentration C_env or, more appropriately, a constant chemical potential μ_env, should be prescribed on a moving boundary ∂ B_env as dictated by the phase field ϕ. In the finite element model this could be enforced by means of Dirichlet, Neumann or mixed boundary conditions <cit.>. However, we note in passing that crack growth rates da/dt in hydrogen containing environments become load independent due to chemical reaction or mass transport limitation. A constant μ_env prescribed on a moving crack surface will fail to predict da/dt measurements if a critical distance is not introduced in the failure criterion; such scheme could require a highly complex implementation. A physically sound framework could be key to quantitatively capture crack growth rates during the stage II plateau; since this regime has not been addressed in the present work, we assume for simplicity a negligible influence of the phase field on hydrogen diffusion.



 §.§ Numerical implementation



The finite element (FE) method is used to solve the coupled mechanical-diffusion-phase field problem. Using Voigt notation, the nodal values of the displacements, phase field and hydrogen concentration are interpolated as follows,

    u=∑_i=1^m N_i u_i   ,          ϕ=∑_i=1^m N_i ϕ_i   ,            C=∑_i=1^m N_i C_i


where m is the number of nodes and N_i are the interpolation matrices - diagonal matrices with the nodal shape functions N_i as components. Accordingly, the corresponding gradient quantities can be discretized by,

    ε=∑_i=1^m B^u_i u_i   ,          ∇ϕ=∑_i=1^m B_i ϕ_i   ,          ∇ C=∑_i=1^m B_i C_i


Here, B_i are vectors with the spatial derivatives of the shape functions and B^u_i denotes the standard strain-displacement matrices.



  §.§.§ FE discretization of the deformation-phase field problem

The weak form for the deformation problem can be written, in its most general form, as

    ∫_Ω( g ( ϕ) σ_0 : δϵ - 𝐛·δ𝐮)  dV - ∫_∂Ω_h( 𝐡·δ𝐮)  dS = 0

where σ_0 is the Cauchy stress tensor of the undamaged solid. Making use of the finite element discretization outlined above and considering that (<ref>) must hold for arbitrary values of δu, the discrete equation corresponding to the equilibrium condition can be expressed as the following residual with respect to the displacement field,

    r_i^u =∫_Ω[(1-ϕ)^2+k] (B_i^u)^Tσ_0 dV - ∫_ΩN_i^T b dV - ∫_∂Ω_hN_i^T h dS


Similarly, the out-of-balance force residual with respect to the evolution of the crack phase field is obtained by discretizing (<ref>) and considering the contribution from the scalar phase field microtraction (<ref>),

    r_i^ϕ= ∫_Ω[ -2(1-ϕ) N_i  H + G_c ( θ) (1ℓ N_i ϕ
        + ℓB_i^T ∇ϕ) ]  dV - ∫_∂Ω_f N_i f  dS

where we have introduced the so-called history variable field H to ensure irreversibility,

    H =
        ψ_0(ϵ)       ifψ_0(ϵ) > H_t
    
          H_t     otherwise

Here, H_t is the previously calculated energy at time increment t. Thus, the history field satisfies the Kuhn-Tucker conditions. The components of the consistent stiffness matrices can be obtained by differentiating the residuals with respect to the incremental nodal variables:

    K_ij^u = ∂r_i^u/∂u_j = ∫_Ω[(1-ϕ)^2+ k] (B_i^u)^T C_0 B_j^u dV


    K_ij^ϕ = ∂ r_i^ϕ∂ϕ_j = ∫_Ω[ ( 2H + G_c ( θ)ℓ) N_i N_j + G_c ( θ) ℓB_i^T B_j ]  dV




  §.§.§ FE discretization of mass transport


Substituting the constitutive equation for the chemical potential (<ref>) into the Fick-type relation for the mass flux (<ref>), one reaches

    J=- DC/(1 - θ_L )(∇ C/C - ∇ N/N) + D/RT C V̅_H ∇σ_H

which, after making the common assumptions of low occupancy (θ_L ≪ 1) and constant interstitial sites concentration (∇ N=0), renders

    J=-D ∇ C + D/RT C V̅_H ∇σ_H

Substituting in (<ref>), the hydrogen transport equation becomes

    ∫_Ω[ δ C ( 1/DdC/dt) + ∇δ C  ∇ C  - ∇δ C ( V̅_H C/RT∇σ_H ) ]  dV = - 1/D∫_∂Ω_qδ C q  dS


A residual vector can be readily obtained by discretizing (<ref>), given that δ C indicates an arbitrary virtual variation of the hydrogen concentration: 

    r_i^C= ∫_Ω[ N_i^T ( 1/DdC/dt) + B_i^T  ∇ C  - B_i^T ( V̅_H C/RT∇σ_H )  ]  dV + 1/D∫_∂Ω_q N_i^T q  dS


From which a diffusivity matrix can be defined,

    K_ij^C = ∫_Ω( B_i^T  B_j  - B_i^T V̅_H/RT∇σ_H N_j   )  dV


where the discretization given in Eq. (<ref>) has also been employed to interpolate the time derivatives of the nodal concentrations. The diffusivity matrix is affected by the gradient of the hydrostatic stress, σ_H, which is computed at the integration points from the nodal displacements, extrapolated to the nodes by means of the shape functions, and subsequently multiplied by B to compute ∇σ_H. We employ 8-node quadrilateral finite elements with C_0 continuity.


At the same time, one can readily identify a concentration capacity matrix,

    M_ij = ∫_Ω N_i^T 1/D N_j  dV


and a diffusion flux vector,

    F_i = - 1/D∫_∂Ω_q N_i^T q   dS


Accordingly, the discretized hydrogen transport equation reads,

    K^C C + MĊ = F




  §.§.§ Coupled scheme


The deformation, diffusion and phase field fracture problems are weakly coupled. First, mechanical deformation impacts diffusion through the stress field, governing the pressure dependence of the bulk chemical potential. Secondly, mass transport affects the fracture resistance via hydrogen buildup in the fracture process zone, reducing the critical energy release rate. And thirdly, the hydrogen-sensitive phase field degrades the strain energy density of the solid.


We solve the linearized finite element system,

    [ K^u   0   0;   0 K^ϕ   0;   0   0 K^C ][ u; ϕ; C ] + 
     [ 0 0 0; 0 0 0; 0 0 M ][ u̇; ϕ̇;  Ċ ]= [ r^u; r^ϕ; r^C ]


by means of a time parametrization and an incremental-iterative scheme in conjunction with the Newton-Raphson method. Due to its robustness, the staggered solution scheme proposed by Miehe et al. <cit.> is adopted to solve the deformation-phase field coupling; a time increment sensitivity analysis is conducted in all computations.


The modeling framework is implemented in the commercial finite element package ABAQUS via a user element subroutine. The code can be downloaded from www.empaneda.com/codes, including documentation describing the details of the implementation in ABAQUS, verification case studies and integration strategies. Post-processing of the results is carried out by means of Abaqus2Matlab <cit.>.



§ RESULTS



A number of case studies of particular interest are addressed to show the capabilities of the modeling framework presented. Firstly, we model the benchmark problem of a cracked square plate subjected to tension, verifying our implementation with the results by Miehe et al. <cit.>, and we subsequently assess the role of hydrogen (Section <ref>). Then, in Section <ref>, hydrogen assisted failure in notched components is investigated, and model predictions are compared with the experimental results of Wang et al. <cit.>. The experimental work by Olden et al. <cit.> is reproduced in Section <ref>, where the model is used to estimate cracking thresholds under constant loading conditions. We investigate internal hydrogen assisted cracking and compare with the compact tension experiments by Thomas et al. <cit.> in Section <ref>. Finally, Section <ref> deals with complex crack propagation from defects that typically arise due to pitting corrosion.



 §.§ Cracked square plate subjected to tension in a hydrogenous environment



We model a square plate with a horizontal crack going from the left side to the center of the specimen. The geometric setup, dimensions and boundary conditions are given in Fig. <ref>. We load the plate by prescribing the vertical displacement in the upper edge, and restrict both vertical and horizontal displacements in the lower boundary. A uniform distribution of hydrogen concentration is assigned throughout the specimen as initial boundary condition: C(t=0)=C_0. In addition, a constant hydrogen concentration C=C_b is prescribed at the boundaries during the entire numerical experiment. All the outer boundaries of the specimen, including the crack faces, are in contact with the environment. These boundary conditions mimic the usual scenario in field applications and laboratory experiments, where the specimen is pre-charged before loading in the same aggressive environment; accordingly, we assume C_b=C_0. For materials with high hydrogen diffusivity, the use of generalized Neumann-type boundary conditions <cit.> or σ_H-dependent Dirichlet boundary conditions <cit.> could provide a more accurate description of the electrochemistry-diffusion interface. As in <cit.>, we adopt the following material properties: Young's modulus E=210 GPa, Poisson's ratio ν=0.3, and critical energy release rate G_c( 0 )=2.7 MPa mm. We assume an iron-based material and consequently adopt a hydrogen damage coefficient of χ=0.89 and a partial molar volume of V̅_H=2000 mm^3/mol <cit.>. For this benchmark study, we assume that the load is prescribed at a sufficiently low rate such that hydrogen equilibrium can be achieved. Thus, from the diffusion coefficient and the test time t_f one can define a diffusion length ℓ_C = √(D t_f), which is taken to be significantly larger than the characteristic length of the fracture problem (ℓ_C >> ℓ). A very fine mesh is employed in the expected crack propagation area after a sensitivity study, ensuring that the characteristic element length h is sufficiently small to resolve the fracture process zone (h< ℓ/7.5). Namely, the square plate is discretized by means of 43110 plane strain quadratic elements.



To assess the capabilities of the model in capturing hydrogen-induced damage we compute the mechanical response under three different environments that range from 1 wt ppm, which corresponds to a 3% NaCl aqueous solution, to 0.1 wt ppm. Fig. <ref> shows the results computed in terms of the load-displacement curve. As shown in the figure, the model shows a strong sensitivity to the environment, with the peak load decreasing with the hydrogen content, in agreement with experimental observations. The numerical results in the absence of hydrogen quantitatively agree with those by Miehe et al. <cit.>, verifying the numerical implementation.



As shown in Fig. <ref>, hydrogen accumulates in the vicinity of the crack as the external load and the hydrostatic stress increase.



Atomic hydrogen concentration lowers the critical energy release rate of the material, facilitating crack initiation and subsequent propagation. Fig. <ref> shows the resulting fracture patterns at different load steps with the red color denoting the case where ϕ=1. The trend depicted by the crack path is the same in all the environments considered. 





 §.§ Hydrogen assisted failure in notched cylindrical bars



We now turn our attention to hydrogen assisted failure from notched specimens under rising load and different environmental conditions. We reproduce the experimental work by Wang et al. <cit.> on cylindrical bars. The plane strain implementation can be readily modified to deal with axisymmetric samples by changing the strain-displacement matrix and integrating the discretized equations in polar coordinates. Fig. <ref> shows the specimens geometry, the mesh employed and the path depicted by the crack that initiates at the notch tip. A total of 52955 axisymmetric 8-node elements are used to discretize the bar. The length of the characteristic element is chosen to be 10 times smaller than ℓ to assure mesh-independent results and resolve the fracture process zone.



We mimic the conditions of the experiments on AISI 4135 steel bars by Wang et al. <cit.>. Uniaxial tension tests are conducted at a nominal strain rate of ε̇=8.3 × 10^-7 s^-1. The experiments are conducted at room temperature. The material properties for AISI 4135 steel are: Young's modulus E=210 GPa, Poisson's ratio ν=0.3, partial molar volume of V̅_H=2000 mm^3/mol, and diffusion coefficient D=3.8 × 10^-5 mm^2/s. The specimens are precharged and subsequently loaded in the same environment; hydrogen concentrations range from 0.08 to 2.2 wt ppm. For each case, the maximum net section stress is computed at the remote boundary. Fig. <ref> shows the failure stress obtained in each numerical experiment as a function of the hydrogen content. The results from Wang et al. <cit.> are also shown; a critical energy release rate of G_c=64 KJ/m^2 provided the best fit to the data.



Phase field predictions reveal a good agreement with the experimental data. Increasing the hydrogen content translates into a reduction of the critical energy release rate, which in turn lowers the strength. Hydrogen accumulates in the vicinity of the notch tip, where tensile hydrostatic stresses are larger. The strength changes drastically within low hydrogen levels but eventually saturates: increasing the hydrogen content beyond 0.5 wt ppm leads to minor changes in the failure stress. The saturation stage is captured due to the presence of the Langmuir-McLean isotherm (<ref>). The embrittlement law seems to capture the experimentally observed response well. Discrepancies are larger in the experiments with very low hydrogen concentrations, where grain boundary decohesion may not be the only mechanism assisting failure.



 §.§ Cracking threshold for duplex steels in sea water under constant load



We investigate the capabilities of the present model to capture hydrogen damage under constant mechanical load by reproducing the experiments by Olden and co-workers <cit.>. In their work, single-edge notched tensile (SENT) specimens were immersed in slowly circulating artificial sea water at 4^∘ C under cathodic protection conditions. The material under consideration is a 25%Cr duplex stainless steel with E=200 GPa and ν=0.3; despite its wide use in sub-sea applications, duplex steels are sensitive to hydrogen embrittlement at low corrosion protection potentials <cit.>. The geometry of the fatigue pre-cracked SENT specimens is shown in Fig. <ref>, along with the finite element mesh employed. A total of 98748 plane strain quadratic elements are used in the discretization. A constant net stress is prescribed at the remote boundaries and a subsurface concentration of 1 wt ppm is assumed in the outer surfaces, based on the experimental conditions (3.5% NaCl solution and a corrosion protection potential of -1050 mV_SCE). The tests have been conducted at room temperature and the diffusion coefficient has been measured as D=3.7 × 10^-6 mm^2/s.



The initial mechanical load introduces a certain amount of damage that is insufficient to propagate the crack. The magnitude of the phase field variable increases and eventually triggers crack growth as time goes by and hydrogen diffuses towards the fracture process zone. Since the mechanical load is kept constant, damage is caused by the reduction in the critical energy release rate with time due to hydrogen diffusion.


We compare our cracking threshold predictions with the experimental data by Olden et al. <cit.>. Results are shown in Fig. <ref> in terms of net stress versus time. As shown in the figure, we are able to reach a good agreement with the tests and extend our predictions over time scales that are impractical in laboratory experiments. This allows us to identify the net stress level below which cracking will not occur: σ_TH=510 MPa. Thus, safe regimes of operation can be efficiently identified by means of the present phase field formulation. However, one should note that the lack of experimental data at larger time scales introduces a degree of uncertainty in the numerical predictions of this case study. Specifically, more tests within the range of 1000-2000 h would have been desirable given the scatter observed for net stresses of 540-560 MPa.










 §.§ Internal hydrogen assisted cracking in compact tension specimens



Hydrogen assisted cracking frequently takes place from hydrogen dissolved within the material prior to loading. This phenomenon is termed internal hydrogen assisted cracking, and atomic hydrogen is typically introduced by manufacturing operations, such as casting, welding, electrochemical machining, or heat treatments. Mechanical loading redistributes the internal hydrogen content to the fracture process zone and promotes cracking.
 

We show the capabilities of the present phase field formulation to model internal hydrogen assisted cracking by mimicking the experiments by Thomas et al. <cit.> on AerMet100. AerMet100 is secondary-hardening quenched and tempered martensitic steel developed for high-performance aerospace applications. An electroplated coating is used for corrosion resistance, but atomic hydrogen can be codeposited in the plating and steel substrate <cit.>. In this case study, we model pre-cracked Compact Tension specimens that have been pre-charged prior to loading. Our aim is to predict the threshold stress intensity for internal hydrogen embrittlement K_TH as a function of the hydrogen content. Fig. <ref> shows the geometry and dimensions of the specimens, following <cit.>, and in agreement with the ASTM E 1820 Standard. Fig. <ref> also shows the mesh (Fig. <ref>) and the characteristic mode I crack propagation path (Fig. <ref>). A total of 34792 quadratic plane strain elements have been used and the stress intensity factor at crack initiation is computed following the ASTM E 1820 Standard. The material properties related to the mechanical response are Young's modulus E=194400 MPa, and Poisson's ratio ν=0.3. The experiments have been conducted at room temperature and the diffusion coefficient for AerMet100 has been measured as D=2 × 10^-4 mm^2/s <cit.>. Mimicking the experiments, a loading rate of K̇=2.2 × 10^-4 MPa√(m)/s is employed. The length parameter is chosen to be ℓ=0.3 mm; 7 times larger than the characteristic element length along the extended crack plane.



Fig. <ref> shows the results obtained in terms of K_TH versus the uniform hydrogen concentration attained prior to loading. A value of G_c=30 KJ/m^2 is assigned to the critical energy release rate to fit the experiments in air. As shown in Fig. <ref>, the numerical predictions qualitatively capture the experimental trends. Namely, the cracking threshold decreases significantly with increasing hydrogen content. Internal hydrogen embrittlement is particularly severe for this very high strength steel, with K_TH being roughly a tenth of the fracture toughness in inert environments. While the agreement with experiments is considered satisfactory, finite element results fail to quantitatively predict the degree of embrittlement at large hydrogen concentrations. Plausibly, the differences are a consequence of fitting G_c with stress intensity measurements that pertain to fracture processes where nucleation, growth and coalescence of microvoids are the mechanisms governing fracture. The value assigned to the trapping energy at grain boundaries may constitute another source of discrepancy; accurate measurements of trap characteristics remain a challenging task.





 §.§ Crack coalescence from corrosion pits



Hydrogen assisted cracking frequently manifests in combination with other corrosion damaging effects. For example, this is often the case in aqueous solutions where cathodic protection has been used to enhance corrosion resistance <cit.>. Cracking in pipelines and other structural components often commences at the sites of corrosion pits, and is accelerated by hydrogen embrittlement. We investigate the capabilities of the present phase field formulation to predict complex crack propagation due to existing defects. This capability could be a game-changer in the energy industry, where structural integrity assessment often involves estimating crack growth and coalescence from well-mapped corrosion pits as a function of time.


We consider a plate with three preexisting defects, as outlined in Fig. <ref>. The defects, which typically originate due to anodic dissolution, are defined in the numerical model by prescribing an initial history field ϕ=1. Introducing pre-existing cracks as heterogeneities in the phase field order parameter (as opposed to discrete geometrical discontinuities) has proven to yield predictions closer to analytic fracture mechanics solutions <cit.>. We assume E=200 GPa, ν=0.3, and energy release rate G_c=90 KJ/m^2. The length parameter is 6 times larger than the characteristic element size, ℓ=0.6 mm. A mesh of 30908 quadratic quadrilateral plane strain elements is employed to model the plate. We consider a diffusion coefficient on the order of that of fcc iron at room temperature, D=1 × 10^-8 mm^2/s, and prescribe an initial hydrogen concentration of C_0=1 wt ppm throughout the specimen. This hydrogen concentration is maintained during the test at the sides of the specimen and in the corrosion pits. As outlined in Fig. <ref>, we prescribe a uniform displacement at the sides, being the loading rate u̇=0.0416 μm/s.



The crack propagation contours are shown as a function of time in Fig. <ref>. Blue and red colors correspond to the completely intact and the fully broken state of the material, respectively. First, damage is restricted to the corrosion pits, which act as stress risers and attract the hydrogen (Fig. <ref>). Then, cracking initiates at the largest defect (Fig. <ref>) and eventually coalesces with the central hole (Fig. <ref>). After arresting at the center defect, the crack reappears at its left corner (Fig. <ref>) and continues to propagate. At the same time, cracking initiates at the smaller pit (Fig. <ref>) and both propagating cracks eventually coalescence (Fig. <ref>). Crack initiation takes place after 5040 s while the final coalescence occurs 2520 s later. This benchmark problem shows the robustness and capacity of the present phase field scheme to model complex cracking conditions that may arise in practical applications; crack deflection, arrest, restart and coalescence are successfully captured.





§ CONCLUSIONS



A phase field formulation for the long-standing problem of hydrogen assisted cracking has been presented. A coupled deformation-diffusion-phase field scheme has been developed in the framework of the finite element method. The model builds upon a hydrogen-dependent surface energy degradation law based on Density Functional Theory that adequately characterizes the sensitivity of the fracture energy to atomic hydrogen concentration. A wide variety of case studies have been addressed to benchmark the potential of the present phase field formulation for hydrogen embrittlement. First, the cracked plate benchmark is modeled to verify the numerical implementation and examine the influence of hydrogen in the load-displacement curve. Then, the failure stress in notched cylindrical bars is computed for a very wide range of environmental conditions, with finite element results showing a good agreement with the experiments. The model is also used to estimate exposure times for fracture in sea water under constant loads; numerical predictions enable to identify the cracking threshold without the need for time consuming laboratory tests. Numerical experiments in Compact Tension specimens are then conducted to compute internal hydrogen stress intensity thresholds as a function of hydrogen content. Finally, the capabilities of the model to capture complex crack propagation arising from defects intrinsic to corrosive environments are investigated. Computations show that phase field methods for fracture are a suitable tool to achieve the elusive goal of lifetime prediction of engineering components undergoing hydrogen assisted cracking.
 


§ ACKNOWLEDGMENTS



E. Martínez-Pañeda acknowledges financial support from the Ministry of Economy and Competitiveness of Spain through grant MAT2014-58738-C3 and the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n^∘ 609405 (COFUNDPostdocDTU). C.F. Niordson acknowledges financial support from the Danish Independent Research Fund (DFF-7017-00121).













elsarticle-num


