Review of hydrogen-assisted cracking models for application to service lifetime prediction and challenges in the oil and gas industry

2.1 Hydrogen pressure theory (HPT)

This theory considers that steel degradation is related to the increase of hydrogen concentration within internal voids, interstitials, or grain boundaries. Once hydrogen atoms in gaseous form saturate the inner sites, there will be a pressure build-up especially at sharp crack tips, resulting in higher internal stresses and promoting void growth and coalescence, which in turn will lead to cracks propagation. This mechanism is often related to blisters formation and coalescence, leading to increasing pressure or stress to fatal failure (Zapffe & Sims, 1941; Tetelman, 1969). However, many observations of stable crack growth in dry hydrogen gas, as well as in chlorine and other gaseous environments, show that the pressure theory is not always at the origin of fatal cracks. Indeed, it was observed that defects such as dislocations induce higher internal pressure in voids even if the source of hydrogen was at very low fugacity. Dislocations may act as pipelines to the diffusion of hydrogen atoms accelerating HAC with lower hydrogen pressure (Roberston et al., 2009; Birnbaum & Sofronis, 1994; Robertson et al., 2001). However, the theoretical predictions of the internal pressure overestimate the arrival rate of hydrogen flows and underestimates its escape from the voids, hence the pressure buildup at the vicinity of the voids is generally overestimated.

2.2 Hydrogen-enhanced decohesion (HEDE)

The HEDE theory considers that hydrogen diffusing with a large amount into grain boundaries, interstitials, and precipitates is weakening the interatomic bonds in metals under stress leading to breakage of atomic binding to accommodate the slip (Oriani, 1972; Lynch, 2012). This mechanism was first introduced to explain, in one hand, the load relaxation during the tensile tests and, on the other hand, the intergranular failure. According to Lynch (2003), dislocations activity may occur during the decohesion, which will increase the local stresses at specific sites such as atomically sharp crack tips due to adsorbed hydrogen and at the vicinity of cracks where dislocation entanglements result in a tensile-stress maximum.

2.3 Hydrogen-enhanced localized plasticity (HELP)

Beachem (1972) was the first to introduce HELP mechanism. According to his model, the high hydrogen concentration dissolving in the lattice ahead of the crack tip assists in the deformation by dislocations and hence increases plasticity. He also considered that metal microstructure, crack tip stress intensity, and concentration of hydrogen are the key factors favoring such mechanism and influencing the failure mode such as intergranular, quasicleavage, or microvoid coalescence. More precisely, Beachem’s model suggests that hydrogen does not inhibit dislocation motion by locking them but rather allows dislocations to multiply and move inside the steels at reduced stresses. Therefore, hydrogen simply permits or assists the normal fracture processes to become operative at unusually low macroscopic strains. Many experimental evidences conducted in transmission electron microscope (TEM) on various metal structures such as bcc, fcc, and hcp support the fact that the mobility of dislocations is enhanced by the presence of hydrogen (Beachem, 1972; Sirois & Birnbaum, 1992; Sofronis & Birnbaum, 1995; Ferreira et al., 1998; Chateau et al., 2002; Robertson et al., 2009).

Typically, hydrogen will diffuse at dislocations core, increasing dramatically the internal stress, which will boost the dislocations mobility and favor the planar deformation by pile-ups dislocations. The increasing number of dislocations will then apply greater stress at the head of dislocation pile-ups, leading to easy glide at lower macroscopic strains and resulting on softening phenomena observed on macroscopic tensile tests (Ferreira et al., 1998; Jouiad et al., 1998). However, it should be noted that dislocation motion’s enhancement by the diffusion of hydrogen is merely a contributing factor to an overall degradation process and not a mechanism on its own. First, there should be enough plastic flow to create dislocation sites for hydrogen to diffuse and enhance the plasticity, and then the HELP mechanism will be prominent.

It is also worth noticing here that when hydrogen is adsorbed at the crack tip and/or at the voids, adsorption-induced dislocation emission (AIDE) mechanism may occur. AIDE is another complex mechanism of HE that includes both dislocations nucleation and emission, leading to more plastic deformation of the steels. HE due to AIDE is related to the action of adsorbed hydrogen, which is weakening the atomic bounds. It is always attributed to cleavage-like and/or dimpled intergranular failures where diffused hydrogen at those sites is at the origin of dislocations multiplication under stress. There will be then impinging dislocations from those failed sites acting as sources of defects emitting toward healthy material. The reader may refer to more detailed reports on these HE mechanisms in the following references (Oriani, 1972; Robertson, 2001; Lynch, 2003).

2.4 Stress-induced hydride formation

This mechanism does not concern iron and nickel alloys. Hydride-forming metals like titanium, zirconium, and vanadium suffer from hydrogen absorption, which causes severe embrittlement. Even at low hydrogen concentration, below the solid solubility limit, stress-assisted hydride formation provokes embrittlement, which is enhanced by slow straining. At higher hydrogen concentrations above the solubility limit, brittle hydrides are precipitated on slip planes and cause severe embrittlement (Pardee & Paton, 1980). This mechanism is also favored by increased strain rates, decreased temperature, and presence of notches in the tested steel. Generally applied stresses to the metal accelerate the formation of hydride. Consequently, hydride formation will recur in the stress field at the crack tip, and the crack will continue to propagate until failure (Kumnick & Johnson, 1980; Gerberich et al., 1986; Huang & Herberich, 1992; Huang & Herberich, 1994; Gerberich, 2012).

Other mechanisms of HE may also occur; however, it is most likely that there is not only one prevailing mechanism responsible for the material failure but rather a real competition between these mechanisms. As the material tends to minimize its internal energy during the failure, all mechanisms may occur simultaneously or alternatively considering the most favorable mechanism with respect to the steel configuration during the damage process, as described by Lynch (2003). According to his assumptions, the steels may undergo a combination of AIDE, HELP, and HEDE mechanisms depending on the fracture mode, which depends on the material structure and other variables such as stress concentration and sample geometry. For example, the dislocations nucleated from the crack tip and explained well via AIDE mechanism may also be explained via HELP mechanism. For both mechanisms, the dislocations movement away from the crack tip will minimize the back-stress on subsequent dislocation emission.

3.1 Model classification

Micromechanical modeling of HAC is a growing area of research and development. Supported by the easy availability of high-performance computing facilities, HAC models are constantly growing in complexity. Far from providing a comprehensive review of the exhaustive literature on micromechanical modeling of HAC (Gangloff, 2003a,b; Gerberich, 2012), the present article focuses on a summary and discussion of key models as well as current emerging trends that might be of benefit for the oil and gas industry.

Models of HAC are varying in complexity in terms of both mechanisms and couplings. Different criteria can be used to compare available models:

1. crack configuration: stationary or growing crack, surface or bulk crack, remote external load (static/dynamic) and/or internal crack pressure build-up (e.g. HIC);

2. source of hydrogen: internal or environmental, H produced at the crack surface or far from the crack surface;

3. hydrogen diffusion: transient or steady state, unidimensional or multi-dimensional, weak or strong coupling with crack tip mechanics, concentration-driven and/or stress-driven, explicit analysis of trapping or analytical approximation, microstructure-dependent traps and/or plastic strain-dependent traps;

4. crack tip mechanics: small deformation or finite deformation theory, linear elastic FM or elastoplastic FM, conventional plasticity (J2) or non-conventional (strain gradient plasticity, dislocation dynamics), weak or strong coupling to local H concentration, explicit stress analysis or simplified analytical solutions;

5. failure mechanism: decohesion-based (HEDE) and/or plasticity-based (HELP) and/or hydrogen-pressure based (HPT).

Available micromechanical and continuum models for HAC are commonly grouped either by failure mechanism (HEDE versus HELP) (Gerberich, 2012) or by type of hydrogen degradation [internal hydrogen assisted cracking (IHAC) versus hydrogen environment assisted cracking (HEAC)] (Gangloff, 2003a,b). In the present review, we adopt another classification aligned with the historical development and evolution of available models: (1) unidimensional models occluding the three-dimensional (3D) distribution of hydrogen concentration and stresses (using semi-analytical approximate solutions) but with the benefit of an explicit relationship between subcritical growth indicators (i.e. K_TH and da/dt), steel properties (yield strength and cohesive stress, H diffusivity), and environment (i.e. H charging pressure/concentration, temperature and applied K) and (2) multidimensional models, which offer a detailed time-dependent analysis of the different chemo-mechanical couplings in the FPZ and leading to advanced two-dimensional (2D)/3D finite element (FE) models.

It is worth noticing that the most available HAC models’ assumptions are based on diffusion-controlled hydrogen permeation to the FPZ. This is quite acceptable when the hydrogen degradation mechanisms occur due to a source of hydrogen far away from the FPZ (i.e. IHAC) because of the relatively large diffusion distance. Fortunately, for the HIC mechanism encountered in the oil and gas industry, which is a bulk cracking phenomenon, HAC models’ assumptions fairly hold. However, for HEAC degradation (such as SCC) where the production of hydrogen happens at the crack surface, the assumption of diffusion-controlled must be carefully verified. To that extent, Turnbull et al. (1996) showed that the use of generalized boundary condition at the crack tip (based on flux and not on constant concentration) allows to properly capture the kinetics of hydrogen at the vicinity of the crack tip and subsequently predict proper crack growth rates. A similar analysis was proposed more recently by Crolet (2016). Analytical models might be quite difficult to revise with the generalized boundary condition proposed by Turnbull et al. (1996); however, this could be introduced in multidimensional numerical models discussed later in the present review.

3.2 Unidimensional approaches to HAC modeling

Table 1 summarizes selected unidimensional models for HAC (with apologies to many authors whose work is not cited). These selected models for HAZ are all based on the decohesion mechanism (HEDE). Note that models of hydrogen diffusion only (i.e. no coupling to any specific hydrogen degradation mechanism) have not been considered.

By coupling a steady-state hydrogen decohesion model to an analytical solution of the elastoplastic stress field ahead of a stationary surface crack exposed to hydrogen pressure p_H2, Akhurst and Baker (1981) proposed a simple model to predict the dependence of K_TH on p_H2. The fundamental assumption of this model is that subcritical cracking occurs (i.e. K_I=K_TH) when the tensile stress σ_yy, exceeds the H-reduced cohesive stress σ^* over a cohesive ligament of length x_C. This critical ligament length x_C along which hydrogen damage nucleates is assumed to be a steel property, which is adjusted to fit predictions with experiments. Historically, this failure criterion finds its roots in the model proposed by Ritchie et al. (1973), better known as the RKR model, which was initially proposed for cleavage failure in the absence of hydrogen and used to predict the relationship between critical tensile stress at failure and the fracture toughness K_IC.

In the presence of hydrogen, Akhurst and Baker (1981) suggested the following relationship for the H-degraded cohesive stress ahead of the crack tip:

where σ0∗ is the hydrogen-free cohesive stress, A is an embrittlement factor, and V_H is the partial molar volume of hydrogen in metal. The distribution of tensile stress σ_yy ahead of the crack tip was taken as a function of the yield stress σ_ys and the hardening exponent n, assuming a conventional J2 plasticity (Schwalbe, 1977):

The model was used successfully to predict the decrease in K_TH with pH2 provided that the fitting parameters A and σ0∗ and x_C are properly selected. This is certainly the main drawback of this simple model, as there is no guarantee on the uniqueness of optimum parameter set. Furthermore, the use of conventional J2 plasticity failed in predicting a sufficiently high level of stress concentration, few nanometers ahead of the crack tip, as necessary to break atomic bonds.

To address the shortcomings of conventional J2 plasticity, Huang and Gerberich (1994) and Gerberich et al. (1986) proposed a different mechanistic model, based on dislocation dynamics, which describes the dependency of threshold intensity K_TH on the local hydrogen concentration at the FPZ, noted C_H-FPZ. Their model was inspired by the pioneering work of Thomson and co-workers (in the absence of hydrogen) on dislocation emission and shielding at the crack tip. The proposed model gave a threshold stress intensity in the following form:

where k_IG is the “local” hydrogen-free steel fracture toughness (in the Griffith sense) expressed as function of the surface energy γ_S as kIG=2EγS, whereas k_IH is the local H-reduced fracture toughness varying linearly with the local hydrogen concentration at the crack tip FPZ as follows: k_IH=αC_H-FPZ. For aqueous hydrogen charging, Gangloff (2003a,b, 2008) suggested the following expression for the hydrogen accumulation at the crack tip FPZ:

where C_S is the hydrogen concertation at the crack surface, σ_H is the local hydrostatic stress, and E_B is the trap binding energy. The values of the other model parameters used in Equation 3 (i.e. α, α′ and β′) were determined from computer simulations. Gerberich et al. (1996) reported good model predictions with respect to experiments carried out on high-strength steels when using α=0.5 MPa⋅m^1/2/atom fraction, α′=2×10⁻⁴ MPa⋅m and 1/β′=5 MPa⋅m^1/2. Not only that this model allowed to address theoretically the deficiency of conventional continuum EPFM in predicting high stress levels in the close proximity of the crack tip (in this model, predicted hydrostatic stress can reach up to 50 times σ_ys), but also, it explained the correlation observed experimentally between the threshold stress intensity K_TH and the hydrogen-free fracture toughness K_IC in high-strength steels [note that in Equation 3, K_IC=K_TH(C_H=0)], as well discussed by Gerberich et al. (1996).

As it pertains to predicting stage II crack growth rate, most unidimensional models rely on a simple ingredient – coupling a transient unidimensional diffusion analysis along the cohesive ligament ahead of the crack tip to a local failure criterion (often base on HEDE). Crack growth is assumed to take place when the opening stress over the critical of length x_c exceeds the material cohesive stress (Van Leeuwen, 1979; Gerberich et al., 1986, 1988, 1991, 1996; Chen & Gerberich, 1991; Toribio & Kharin, 1997; Hall and Symons, 2001). The unidimensional diffusion analysis is used to determine the critical time t_c for hydrogen to diffuse from C_S at the crack surface over the critical distance x_c to reach a critical hydrogen concentration C_CRIT in the FPZ. Stage II growth rate is then simply derived as the ratio x_c/t_c. In other words, subcritical crack growth is assumed to take place in the form of successive (instantaneous) crack jumps of length x_c with a period between successive jumps being t_c. Hydrogen equilibrium on the newly created crack surface is assumed to take place immediately upon the crack jump, and therefore, a constant hydrogen concentration C_S is prescribed as a boundary condition at the newly created crack surfaces. Note that we have limited the present review to growth rate models where hydrogen diffusivity is the rate-limiting step for crack growth. This assumption is valid for internal hydrogen embrittlement (IHAC); however, it may not hold for environment hydrogen embrittlement (HEAC), as extensively discussed by Gangloff (2003a,b).

Gangloff (2003a,b, 2008) proposed a generalized formalism encompassing the most available diffusion-controlled models for stage II growth rate, as follows:

where D_Heff is the effective (trapping) hydrogen diffusion coefficient, and ξ is a function of the above-indicated variables giving an output from 0.01 to 4, depending on the model. Note that in order to provide a valid formalism for both IHAC and HEAC, Gerberich et al. (1996) replaced the factor 1xc by xaxd2, where x_a is the elementary crack advance and x_d is the diffusion distance (i.e. the average distance travelled by hydrogen to accumulate in the FPZ), which equals x_a in the case of HEAC (and reduces to Equation 5 with x_a=x_c).

Extensive HEAC experiments carried out by Gangloff (2003a,b) at the University of Virginia showed that the fastest measured stage II crack growth rats (da/dt)_II correlate (almost linearly) with the best available measured values of trapped hydrogen diffusivity D_Heff for a wide variety of high-strength steels tested in both aqueous (chloride solutions) as well as gaseous (H₂ and H₂S) environments (see Figure 1). This important finding not only supported the assumption of diffusion control but also provided a mean for estimating the critical “steel” parameter x_c. Using the crack growth rate measurements shown in Figure 2, Gangloff suggested that ξ shall not exceed 3, which yields to x_c=0.7 μm.

Figure 1:

The dependence of (da/dt)_II on D_Heff for high-strength alloys that exhibit HEAC in gases and electrolytes at 250°C. High-strength austenitic stainless steel and nickel superalloys were cracked in high pressure (100–200 MPa) H₂, while maraging and tempered-martensitic steels were cracked in low pressure (~100 kPa) H₂. The dotted line represents TG cracking of Fe-3%Si single crystal in 100 kPa H₂ at 0–125°C. Filled symbols (■) represent the transition from molecular to atomic hydrogen gas. Reprinted from Gangloff (2003a,b), with permission from Elsevier.

Figure 2:

Crack propagation sequence as per the model proposed by Unger (1989). Reprinted from Unger (1989), with permission from Elsevier. The labels A–E represent the steps of the crack opening.

Using a different approach, Neimitz and Aifantis (1985, 1987a,b) proposed an incremental subcritical growth model based on a Barenblatt cohesive zone (ahead of a growing crack) that degrades with time due to the accumulation of hydrogen in the FPZ. Almost in parallel, Unger (1989) proposed an incremental HAC model using a Barenblatt/Dugdale cohesive zone approach similar to that proposed by Neimitz and Aifantis (1985, 1987a,b). Unger’s model, however, introduced a pointwise degradation of the cohesive force σ^*(x, t) ahead of the crack tip by assuming a linear decrease with the local hydrogen concentration (σ∗(x, t)=σ0∗−γC(x, t)), which necessitated a numerical resolution of the transient coupled chemo-mechanical problem. In a subsequent work, Lee and Unger (1988) and Unger (1989) proposed a simplified (analytical) solution to the problem by assuming that the extent of the cohesive zone is small compared to the crack length. Note, however, that Unger’s model used a hydrogen-free critical crack-tip opening displacement (CTOD), in contrast to the model proposed earlier by Neimitz and Aifantis (1985, 1987a,b). Figure 3 illustrates schematically the incremental model proposed by Lee and Unger (1988). Initially (a), a constant hydrogen concentration C₀ is prescribed at the crack tip located at x=a₀ and a remote mode I opening stress σ_∞ is applied, such that the resulting CTOD δ is below the critical CTOD δ_c (or, equivalently, K_I<K_IC). With increasing hydrogen accumulation in the FPZ over time (b), the cohesive stress is locally reduced, which is compensated by an increase in the cohesive zone length, thus increasing the crack opening δ. When the degradation is sufficient (c) to increase δ to its critical value δ_c (assumed to be independent of H content), crack extension occurs and is accompanied by a drop in δ due to the sudden increase of the cohesive stress ahead of the newly created crack front as it enters the healthy portion of the steel. Crack will arrest (d) when most of the hydrogen damage zone has passed and the material ahead the new crack tip is essentially healthy. The crack opening δ drops then to an equilibrium value δ_A<δ_c. The hydrogen concentration at the newly created crack tip (x=a₁) is then fixed to C₀ (e), assuming fast equilibrium takes place immediately after the crack jump. The cycle (e)-(c)-(d)-(e) continues over time causing incremental crack jumps until the length of the crack reaches a critical value where δ=δ_c even in the absence of hydrogen (C₀=0) and unstable failure takes place. This incremental model successfully captured the characteristic three-stage curve of da/dt vs. K using γC₀ as the only fitting parameter. Unger (1995) proposed some insights on how to identify this parameter from standard threshold stress intensity laboratory measurements. The main limitation of this model is the absence of hydrogen trapping as well as the absence of plasticity effects. This necessitates the use of a FE approach as described in the next section.

Figure 3:

Typical calculated crack advancement with time using the cohesive model by Serebrinsky et al. (2004). L_coh=x_v−x_f denotes the length of the cohesive zone. x_f Locates the traction-free point, behind which the crack flanks experience zero traction. x_v Locates the vertex of the cohesive zone, ahead of which there is zero opening displacement. Reprinted from Serebrinsky et al. (2004), with permission from Elsevier.

3.3 Multidimensional approaches to HAC modeling

Selected multidimensional models suitable for the prediction of threshold stress intensity and/or subcritical HAC growth rates are presented in Table 2. The literature review clearly shows that the cohesive zone modeling (CZM) approach has been extensively used over the last decade in HAC modeling. This method, initially introduced by Needleman (1987) to study particle matrix decohesion (in the absence of hydrogen), is particularly well adapted to decohesion/delamination type of failure and has gained popularity for being implemented in most commercial FEM packages (e.g. ABAQUS, ANSYS, COMSOL). In few words, CZM consists of introducing “interface elements” in a thin strip zone of zero (negligible) thickness aligned along the (a-priori known) crack path. Upon applying an external load, interface elements start to separate following a pre-defined relationship called traction-separation law (TSL), which relates the amount of local separation/opening δ to the local cohesive stress in the element σ^*. Complete separation at the interface occurs when a critical opening δ_c has been reached. A total energy E^* is dissipated locally during the separation process, which corresponds to the work of separation: E∗=∫0δcσ∗(s) ds. The main advantage of CZM is its ability to separate the bulk (“global”) steel behavior (as described by the macroscopic stress-strain curve) from the “local” behavior of the cohesive interface (as described by the TSL) while maintaining a clean coupling between both behaviors in a continuum fashion. This makes CZM particularly attractive for HAC modeling, as it provides the localized environment ahead of the crack tip where HE mostly takes place. In addition, the decorrelation between the local and global steel behaviors offered by CZM is quite suitable to capture the local hydrogen effect on atomic bonding energy in the FPZ independently of its impact far from the FPZ (i.e. on the plastic behavior of the bulk steel).

Serebrinsky et al. (2004) introduced a quite comprehensive continuum model for HAC to simulate brittle decohesion in bcc iron systems using quantum-based cohesive elements. The proposed model considers a strong and time-dependent coupling between hydrogen transport and the local degradation of the cohesive zone ahead of a growing crack. The hydrogen transport in the pre-cracked continuum was described by a stress-driven Fickian diffusion (neglecting explicitly trapped hydrogen concentration) as follows:

where V_H is the partial molar volume of hydrogen in the metal solid solution and σ_H is the hydrostatic stress (–σ_ii/3). The effect of trap sites on hydrogen diffusion was implicitly considered through an effective diffusivity coefficient (D_Heff). The value of the hydrostatic stress was derived from a static stress analysis assuming a conventional J2 elastoplastic steel behavior with power law strain hardening. Crack extension was driven by brittle decohesion along {110} planes and was modeled by cohesive elements with a linear hydrogen-affected TSL, fully described by two parameters: the critical cohesive stress  σc∗ and critical opening δ_c. The dependence of the surface energy (ideal work of fracture) on local hydrogen coverage θ was derived from first principles ab initio simulations carried out by Jiang and Carter (2004). Assuming that the critical opening δ_c is insensitive to hydrogen coverage, the critical cohesive stress decreases with hydrogen coverage following the same dependency as the surface energy presented in Equation 7. The hydrogen coverage θ was related to the bulk hydrogen centration C using Langmuir-Mc Lean isotherm.

The model could capture qualitatively and quantitatively stage I and II subcritical growth rates on AISI 4340 steel tested in different aqueous environments. The model results also show some interesting findings, namely, the finite crack jump at initiation (i.e. when K_I=K_TH) as well as the intermittent stepwise-like crack growth (in the order of 40-μm jumps), as shown in Figure 3. This important finding not only supports diffusion control crack growth assumption but is also in agreement with earlier results (discussed in the previous section) regarding the existence of a critical length scale x_c, which identifies to the discrete crack jump in HEAC. Furthermore, the model could reproduce the dependence of threshold stress intensity and stage II crack growth rate on environment, temperature, and yield strength.

Olden and co-workers used the model proposed by Serebrinsky et al. (2004) to simulate HEAC of 25% Cr duplex stainless steels in aqueous environment (Olden et al., 2007, 2009a,b) and more recently (in 3D) for API-5L-X70 steel welds (Olden et al., 2012; Alvaro et al., 2014). To account for the ductile failure in this class of lower strength steels, the authors introduced a polynomial traction separation law as follows:

The interesting novelty in the work of Olden’s group is that the hydrogen coverage θ was calculated from the “total” bulk hydrogen concentration, i.e. the sum of the lattice hydrogen concentration C_L (as computed from the stress-driven diffusion analysis) and the trapped hydrogen concentration C_T. By fitting numerical results reported by Taha and Sofronis (2001) on hydrogen transport in the vicinity of crack tip, Olden et al. (2007) overcome the cumbersome resolution of an explicit hydrogen trapping problem and proposed the following linear relationship that correlates the trapped hydrogen concentration C_T to the local cumulative plastic strain ε_p and lattice concentration C_L:

If the parameters of the TSL (i.e. the initial critical cohesive stress σc∗(0) and the critical opening at separation δ_c) as well as the effective hydrogen diffusivity are well selected, Olden et al.’s (2007) model can properly capture the incubation time to failure, threshold stress intensity, and subcritical growth rates. For 25% Cr duplex stainless steel, Olden et al. (2007) suggested an initial critical cohesive stress in the order of 3.5σ_ys and a critical opening in the order of 20 μm, while for X70 steel (base metal), the proposed values are, respectively, 3.1σ_ys and 0.3 mm. The proposed values resulted in cohesive energies very close to experimental values derived from conventional J-R testing.

To reconcile HEDE and HELP in a unified model and linking the physical mechanisms of HE at the microscale to observable/measurable macroscopic quantities, Sofronis and co-workers introduced the concept of a unit cell model. Liang and Sofronis (2003) proposed a first phenomenological model to study the effect of hydrogen on void nucleation and debonding at the interface of an elastic inclusion embedded in an elastoplastic matrix. The aim of this model was to characterize the effect of hydrogen on the cohesive interface properties and its impact on the macroscopic stress at the onset of decohesion. The phenomenological approach consisted of coupling a transient model for hydrogen transport [accounting explicitly for trapped hydrogen using the trapping model proposed by Sofronis and McMeeking (1989), later modified by Krom et al. (1999)] to material deformation (local softening in the presence of hydrogen) and interfacial decohesion (reduction in the cohesive stress).

To account for the hydrogen-induced shear localization (HELP) and local softening of the elastoplastic matrix in the presence of hydrogen, which is modeled at the continuum scale through a reduction in the local flow stress σ_Y, Liang and Sofronis (2003) used the following relationship developed by the same authors in an earlier work of Sofronis et al. (2001):

where C is the local total hydrogen concentration (lattice and trapped), F_p is the strain hardening function (taken as a power law of the cumulative plastic deformation ε_p), and σ_Y0(C) is the yield stress in the presence of hydrogen. The latter is taken as a decreasing function of the local hydrogen concentration as follows:

where σ_Y0 is the hydrogen-free yield stress and ξ≤1 is a softening factor.

To account for interfacial decohesion at the inclusion/matrix boundary, cohesive FEs were introduced at this interface. The authors used a 2D polynomial TSL (in the sense that both normal and tangential openings were considered) with three parameters δ_n, δ_t and σc∗(Γ) being, respectively, the critical normal opening (in pure normal opening mode), the critical tangential opening (in pure tangential opening mode), and the critical cohesive stress in the presence of surface segregated hydrogen at surface concentration Γ. The cohesive model parameters were then calibrated such that the reversible work of separation is exactly the same as that given by the thermodynamic interfacial decohesion model (in the presence of impurities) proposed by Hirth and Rice (1980) and Rice and Wang (1989). Full details on the model calibration procedure are given by Taha and Sofronis (2001).

In a subsequent work by Liang and Sofronis (2004), the authors introduced another complexity to their original model by adding the effect of grain boundary decohesion. The model was essentially the same – a unit cell system comprising an elastic inclusion embedded in a ductile matrix with a grain boundary passing at the axis of symmetry and loaded in the direction normal to the grain boundary plane. The model was used to simulate intergranular cracking of alloy 690. The simulation results indicated that embrittlement is essentially controlled by the hydrogen-enhanced debonding at the carbide/matrix interface, whereas the effect of hydrogen in the grain boundary decohesion was found to be almost insignificant. Furthermore, the model gave means to estimate the fracture energy at crack initiation in the presence of hydrogen.

The unit cell model was also used by Liang et al. (2004, 2008) to study the influence of hydrogen on void growth during ductile failure. In this study, the model was essentially the same as that used for void nucleation at second phase interface, except that in the study by Liang et al. (2004, 2008), no cohesive elements were used as decohesion along grain boundaries was not considered. The authors simulated the stress distribution in a unit cell comprising a spherical void in a hydrogen-softened elastoplastic matrix, loaded under a remote tensile stress. The main target was to assess the effect of hydrogen on the macroscopic (homogenized) response of the cell as described by the stress-strain curve during ductile failure. The main result was that hydrogen assists void growth regardless of the intensity of trapping and material softening but with some dependency to the stress triaxiality. The main gap filled by the unit void cell model was to provide the theoretical foundations to bridge two scales, i.e. the microstructural degradation and the macroscopic observed embrittlement. As shown in Figure 4, the main idea is that phenomenological traction separation laws issued from a unit cell model can be directly integrated into a larger macroscale model of subcritical crack growth. This has been achieved by Ahn et al. (2007a,b) and successfully used to predict subcritical ductile failure in A533B pressure vessel steel.

Figure 4:

Two-scale modeling of hydrogen-assisted ductile crack propagation proposed by Ahn et al. (2007a,b). (A) Crack tip fracture process zone. (B) Axisymmetric void cell model. (C) Traction-separation law. (D) Cohesive element. (E) Cohesive elements characterized by a traction-separation law based on the void cell model. Reprinted from Ahn et al. (2007a,b), with permission from Springer.

Note that recent investigations by Dadfarnia et al. (ASTM E1681-03, 2003; Dadfarnia et al., 2009, 2014) raised questions around the fidelity and accuracy of HAC results issued from cohesive zone modeling simulations. The authors observed that crack growth required a continuously increasing FPZ and associated computed crack growth velocities were, in general, slower than those observed experimentally. In addition, the method was found not to preserve the accuracy and fidelity of hydrogen concentration and stress fields very close to the crack tip. These challenges led Dadfarnia et al. (2014) to propose a model based on the standard node release technique and the RKR (Gerberich et al., 1986) failure criterion discussed in the previous section.

4.1 First principles approach

In a multiscale computational framework, ab initio is the most basic and proper way in which the treatment of phenomena such as bond-breaking, charge transfer, electronic and optical excitations, etc. can be tackled. This can be accomplished using quantum mechanics (QM), which, owing to its high computational demand, application should be restricted to a relatively small system of a few 100 atoms. In QM computational schemes, the properties of many electron system can be properly addressed by using density functional theory (DFT), most notably in most cases, the Kohn and Sham (1965) DFT. This fact obviously restricts the augmentation of the physical problem at a small number of atoms, which calls for pairing the QM method with molecular mechanics (MM) methods, which, although empirical in nature, offer the capability of describing elastic deformations, small-amplitude vibrations, etc. The MM methods can tackle millions of atoms or more, therefore offering significant advantage over employing single QM computational schemes. A noteworthy example of empirical potential for the MM method is the embedded-atom method (EAM). Therefore, computational schemes that combine QM and MM methods offer promising solutions to challenging problems requiring atomistic simulations of materials (Lin & Truhlar, 2006; Bermstein et al., 2009; Zhang et al., 2012). In this context and employing coupled DFT and EAM potentials, Lu and Kaxiras (2004) and Choly et al. (2005) studied a different dislocation behavior (Figure 6).

Figure 6:

In hybrid QM/MM modeling, the MM domain encompasses the QM domain in a seamless fashion.

Many works have been performed using first-principles approach to model hydrogen effect on Fe systems. For the first time, the reduction of cohesive energy in iron grain boundaries was reported by Zhong et al. (2000), using first-principles calculations. Other works (Tateyama & Ohno, 2003; Lu and Kaxiras, 2005) that are based on first-principles approach study the vacancy-hydrogen interactions and correlate the findings with HE. In a more recent first-principles study, Yamaguchi et al. (2011) argued that cohesive energy of bcc Fe Σ3(111) and fcc Al Σ5(012) grain boundaries can be significantly reduced if a large amount of hydrogen atoms segregate in those grain boundaries. Whereas in the study by Yamaguchi et al. (2011) and Zhong et al. (2000), the cohesive energy of grain boundaries was studied, Itakura et al. (2013) focused on assessing the effect of hydrogen on the mobility of a screw dislocation in bcc Fe using first-principles approach. Counts et al. (2010, 2011) provide a series of DFT computations to investigate the interaction of single and multiple H atoms in a body-centered cubic (bcc) Fe. They show that H-H interactions are weak as well as their maximum binding energy, while they investigate the binding energy of H to defects.

Similarly, Desai et al. (2010) studied the mechanistic aspects of hydrogen trapping in a Fe-Y alloy, using first-principles approach. Extending the work performed by the aforementioned authors, Matsumoto et al. (2014) performed MM simulations involving crack growth, nanoindentation, and tensile loading of a poly crystalline nanorod. They showed that in the case of crack growth and incorporating pseudo-hydrogen effects, the steels exhibited brittle behavior in nanoindentation simulations. Again, incorporating hydrogen pseudo-effects, the mechanism that was identified to be in play was HELP, and lastly, for tensile loading simulations, incorporating hydrogen pseudo-effects, the specimen fractured along the grain boundaries.

One of the most cited work regarding first-principles computation of fracture energy of Fe and Al due to hydrogen uptake is that of Jiang and Carter (2010) and Horstemeyer (2009). The authors assessed the ideal fracture energies of iron and aluminum in the presence of varying amounts of hydrogen, along with providing some insight into the cohesion-reduction mechanism of HE in metals. This work has been the foundation of many later works (Tadmor et al., 1996; Olden et al., 2008a,b, 2009a,b; Alvaro et al., 2014), which develop a macro-model, incorporating the results produced by Horstemeyer (2009).

Many times, ab initio calculations or QM/MMs are correlated with in situ testing and observation to validate the results produced by the multiscale computations in the lower scales. Tests like in situ transmission electron microscopy (TEM), where tests can be performed and observed at the nanoscale (dislocations, atoms, etc.), electrochemical nanoindentation, micropillar compression testing, which offers the possibility of quantifying the slip behavior of both single grains and grain boundaries, etc. They all can be correlated to diffusion and trapping parameters using an electrochemical permeation cell and a thermal desorption spectroscopy (TDS) system, respectively.

The ab initio and MM methods are well suited for studying phenomena at the nanoscale. It is very hard to link all the scales from QM scale to macroscale, so there is a real need to extend these methods to account for equipment scale engineering problems and structures.

4.2 Quasi-continuum density functional theory

A successful concurrent multiscale method that can address the phenomena that occur at the atomic scale with the accuracy of DFT and at the same time tackle length scales relevant to real-size structures or specimen testing is the quasi-continuum density functional theory, or simply QC, which was originally proposed by Tadmor et al. (Lu et al., 2006). The QC method combines atomistic methods with continuum models, which, compared to standalone atomistic or continuum simulations, offer great advantages in terms of computational efficiency and accuracy. The concept behind the QC method is that the processes that occur in the atomic level often occur in a very small domain, while the behavior of the rest of the atoms in the bulk of the steel can be modeled using well-established continuum theories. This underlying concept is exploited in the QC method by using atomistic modeling in the areas of interest while modeling the rest of the steel using continuum FE description. The present QC-DFT approach (Fivel, 2008; Peng et al., 2008; Peng & Lu, 2011) employs the combined QM and MM treatment, which was previously discussed and exhibited remarkable adaptability to changing circumstances, such as nucleation of new defects or migration of existing ones.

4.3 Discrete dislocation dynamics

While MM or MD methods can be employed to study one or a few dislocation lines and their interaction with other phenomena, e.g. hydrogen, the scale of the domain is small and MD is mostly used to study the properties of the dislocations in a model with volume not greater than 200 nm³ (Devincre et al., 1992). Many times, there is a need to study a representative number of dislocations present in crystal and their interaction to derive useful conclusions regarding, e.g. dislocation mobility during plastic deformation, dislocation mobility with hydrogen residing in the dislocations, etc. Discrete dislocation dynamics (DDD) can tackle these problems by treating each dislocation explicitly and having elastic properties embedded in a crystal with elastic properties. External loadings are passed on to the dislocations through the crystal they are embedded in. In this sense, length scales up to 50 μm³ can be simulated using DDD and a large number of dislocations and their interaction can be modeled (Devincre et al., 1992). At higher length scales, CM models are used to solve the engineering problem, while DDD simulations are used to feed constitutive equations to the CM models. Three-dimensional discrete dislocation dynamics is the brain child of Kubin and Canova (1992) and Li (1964). The first example dealt with an fcc single crystal containing sets of edge and screw dislocations segments in a continuum media. Each finite dislocation segment generates a long-range elastic stress field, which, for the case of isotropic elasticity, was evaluated following the works of Li (De Wit, 1967) and De Wit (Willis, 1970). When anisotropy is at play, the line integrals proposed by Willis (Jothi et al., 2015a,b,c) are used. In this framework, the stress applied at each dislocation is effectively considered at the middle point of each segment. This stress is evaluated as the superposition of the internal stress by all other segments in the modeled 3D volume and the applied stress induced by the external loading. In fact, this induces a force that is expressed by the Peach-Koehler equation. DDD simulations can offer statistical data regarding dislocation densities, cumulated shear strain, stored energy, etc. (Devincre et al., 1992).

4.4 Microscale to macroscale modeling

A more computationally efficient approach, which many times is used in polycrystalline steels, e.g. steels, Al alloys, etc., is the coupling of a microscale model, with length scales comparable to the grain size of the studied steels (from a few to hundreds of micrometers) and the component size macroscale model. At the microscale, a polycrystalline unit cell, or representative volume element (RVE), is modeled, which contains grains and grain boundaries as well as other potential features such as crack, defects, voids, etc. (Benedetti & Aliabadi, 2015) (Figure 5). RVE can play a two-fold role. Either it can be used to determine the steel properties at the macro scale or it can be used to model the region of interest, e.g. a propagating crack, where the phenomena can be modeled at the microscale. The coupling of the two scales can be achieved, e.g. by transferring macro strains to the RVE as periodic boundary conditions while transferring the averaged micro stresses as macro stresses (Groeber et al., 2006). Another important feature of this method is the fact that both models are relatively easy to construct; the macro model mostly from CAD models or drawings of the specimens, and the microscale model can be constructed using a Dual-beam focused ion beam scanning electron microscope (FIB/SEM) combined with an electron backscatter diffraction (EBSD) system (Jothi et al., 2014a,b). Several works have focused on hydrogen-related phenomena, using the aforementioned method, such as the works of API 579-1 (2007), Rimoli and Ortiz (2010), Ebihara et al. (2011), Jothi et al. (2014a,b, 2015a,b,c), and Benedetti and Aliabadi (2015).

6.1 Model parameter identification

6.2 Limitations of fracture mechanics similitude (specimen versus structure)