A Physically Based Model Predicting the Degradation of Hydrogen on Crack Growth Critical Stress Intensity Factor of Metals

: A simple, physically based model is developed to quantitatively predict the degradation of hydrogen on the crack growth critical stress intensity factor (CSIF) of metals. The model is formulated by combining a microscopically shielded Grifﬁth criterion (MSGC) model for plasticity-induced cleavage fracture and thermodynamics decohesion (TDD) theory for hydrogen-enhanced interface decohesion. The hydrogen-inﬂuenced CSIF is described as a function of the intrinsic CSIF (hydrogen-free), initial hydrogen concentration (solubility), hydrogen trap binding energy and crack tip stress. All parameters in the model can be determined with a physical basis and the model is successfully validated by comparison with published experimental data.


Introduction
Crack growth criteria in fracture mechanics have been widely used for the design and safety assessment of metallic components.Particularly, the K criterion, which defines a crack under mode I loading in a component will not grow if K I ≤ K IC , where K I is the stress intensity factor (SIF) of the crack and K IC is the critical stress intensity factor (CSIF, also known as fracture toughness) of material used for the component.However, the components serving in hydrogen-containing environments often suffer from hydrogenassisted fracture (HAF) under lower stresses (K I ≤ K IC ), leading to "brittle" failures.It is necessary to develop a physically based model to accurately predict the degradation of hydrogen on fracture critical parameters, i.e., the CSIF with hydrogen influence (K IH ).
The thermodynamics decohesion (TDD) theory of interface separations induced by adsorbed solute atmospheres was established by Hirth and Rice [25] and Wang [26].It provides criterion for interfacial decohesion from the view of thermodynamics of adsorption.In recent years, the TDD-based criterion has been implemented in coupled mechanicaldiffusion finite element analysis to calculate stresses, strains and hydrogen concentration field around a crack tip or an inclusion or in molecular dynamics simulations of crack tip and has been found capable of reproducing the experimental trends [27][28][29][30].
In this paper, we combine the MSGC model and TDD theory to form a physically based model to study the crack along cleavage planes or grain boundaries in the context of HAF.In this model, all parameters are physically based and can be measured or calculated in advance.More importantly, the model is simple and convenient to provide accurate predictions which are validated by comparing with published experimental data.

Methodology
In order to evaluate the role of crack tip emission and shielding of near-tip dislocations in cleavage fracture of semibrittle materials, Gerberich et al. [31][32][33] developed a discretized dislocation model consisting of five single dislocations near the tip and a giant superdislocation for far field.Based on two-dimensional computational simulations using the microscopic Griffith fracture criterion, the MSGC-based model is developed as [31-33]: where K IC is the CSIF under mode I loading, k IG is the local microscopic Griffith stress intensity for fracture in plane strain, σ ys is yield strength and α and β are crack tip stress field constants, which can be determined from the computational simulations.Thus, the interplay between crack tip plasticity and stress-based decohesion is captured by the MSGC model to give accurate an prediction of K IC [1][2][3][31][32][33].In order to extend the model to include hydrogen influence, Gerberich et al. [1][2][3] assumed that the local Griffith stress intensity decreases linearly with the increment of local hydrogen concentration, as: where K IH is the CSIF with hydrogen influence, C H is the hydrogen concentration at the potential cracking site ahead of the crack tip, and α is a constant representing the decreasing rate of k IG with the increment of hydrogen concentration.That is, hydrogen-influenced local microscopic Griffith stress intensity k IGH is assumed as k IGH = k IG − αC H .It is worth nothing that α is not physically determined and is often obtained by fitting against experimental data for a given material [1,2,[6][7][8][9][10]12,[16][17][18][19][20][21][22].
According to the microscopic Griffith fracture criterion, for the separation of an interface, the k lG can be expressed as a function of the effective work of decohesion γ eff of the interface, which is the sum of the reversible work of decohesion 2γ 0 int and plastic work γ P .If plastic work is not considered for semibrittle materials, γ eff equals to 2γ 0 int , therefore [32] k where E is Young's modulus, ν is Poisson rate and 2γ 0 int is the reversible work of decohesion in absence of hydrogen.According to the Langmuir-McLean adsorption isotherm, Novak et al. [28] and Jung et al. [29] proposed a physically based expression for the TDD theory built by Hirth and Rice [25] and Wang [26], describing the degradation of cohesive energy when hydrogen is introduced: where 2γ H int is hydrogen-influenced cohesive energy, ∆g 0 int and ∆g 0 FS (∆g 0 = ∆h 0 − T∆S 0 , ∆h 0 is enthalpy change, ∆S 0 is entropy change, and T is temperature) are Gibbs free energy excesses when hydrogen is segregated to the interface and absorbed onto the free surfaces (FS) created by separation of the interface, respectively, and Γ is hydrogen coverage of the interface defined as the number of atoms per unit area of interface.With the consideration of a monolayer hydrogen coverage, the hydrogen coverage Γ can be related to hydrogen occupancy θ H of the interface (taken as hydrogen traps) by Γ/Γ max = θ H [25][26][27][28][29], where Γ max is the saturated hydrogen coverage.Therefore, Equation (4) can be rewritten as where Equation ( 5) can be verified by first-principle calculations (density functional theory methods).As shown in Figure 1, Jiang and Carter [34] obtained the cohesive energy of Fe(110) and Al(111) crystal interface using first-principle calculations, and Jensen et al. [35] calculated the cohesive energy of Ni-∑5(012) grain boundary.It is shown that when θ H < 1 (monolayer coverage), cohesive energy is reduced approximately linearly by increase of θ H .The linear coefficient η obtained by fitting the data using Equation ( 4) is shown in Figure 1.
.If plastic work is not considered for semibrittle materials,  equals to 2 , therefore [32] where  is Young's modulus,  is Poisson rate and 2 is the reversible work of decohesion in absence of hydrogen.According to the Langmuir-McLean adsorption isotherm, Novak et al. [28] and Jung et al. [29] proposed a physically based expression for the TDD theory built by Hirth and Rice [25] and Wang [26], describing the degradation of cohesive energy when hydrogen is introduced: where 2 is hydrogen-influenced cohesive energy, Δ and Δ ( Δ = ∆ℎ − ∆ , ∆ℎ is enthalpy change, ∆ is entropy change, and  is temperature) are Gibbs free energy excesses when hydrogen is segregated to the interface and absorbed onto the free surfaces (FS) created by separation of the interface, respectively, and  is hydrogen coverage of the interface defined as the number of atoms per unit area of interface.With the consideration of a monolayer hydrogen coverage, the hydrogen coverage  can be related to hydrogen occupancy  of the interface (taken as hydrogen traps) by / =  [25][26][27][28][29], where  is the saturated hydrogen coverage.Therefore, Equation ( 4) can be rewritten as where Equation ( 5) can be verified by first-principle calculations (density functional theory methods).As shown in Figure 1, Jiang and Carter [34] obtained the cohesive energy of Fe(110) and Al(111) crystal interface using first-principle calculations, and Jensen et al. [35] calculated the cohesive energy of Ni-∑5(012) grain boundary.It is shown that when  < 1 (monolayer coverage), cohesive energy is reduced approximately linearly by increase of  .The linear coefficient  obtained by fitting the data using Equation ( 4) is shown in Figure 1.[34] and Ni-∑5(012) grain boundary [35] (Ni-∑5(012) grain boundary has three sets of data due to different calculation schemes used by the authors).According to Equation (3), the hydrogen-influenced local microscopic Griffith stress intensity k IGH can be expressed as k IGH = k IG (1 − ηθ H ). Therefore, the hydrogen- influenced CSIF K IH can be formulated as: Incorporating Equation ( 1), K IH can be scaled by K IC as: In addition, taking the logarithms of Equations ( 1) and ( 7) leads to ln(β Equation ( 9) is the proposed model in this paper, which is very simple and describes the K IH as a function of K IC , η and θ H . Based on the Langmuir-McLean isotherm for surface adsorption [36], hydrogen occupancy of interface (θ H ) can be related to the concentration in lattice sites [26,28,37,38] (note that, ), as: where E b is hydrogen binding energy with the interface, R = 8.314 J/(mol•K) is gas constant, V H is the polar molar volume of hydrogen, C 0 is the initial lattice hydrogen concentration (or the hydrogen solubility in equilibrium with the environment in atom fraction in absence of stress), and σ H is the hydrostatic stress at crack tip.Equation ( 10) considers the enhancement of stress and hydrogen trapping on local hydrogen concentration.For a given material, β and K IC are constants, η and θ H have their physical origins and can be calculated by Equations ( 6) and (10), respectively.Therefore, the model can be used to predict the hydrogen-influenced CSIF directly.

Model Validation
In order to validate the model in predicting K IH , parameters including ∆g 0 int , ∆g 0 FS , Γ max and E b , etc., needed to be determined beforehand.Although all these parameters are physically based, not all of them can be collected precisely for a given material at present due to inadequate experimental work.Table 1 shows the parameters recommended for Inconel 690 alloy by Liang et al. [27] and high strength martensitic steels [28,30] when the TDD theory was implemented in their HAF numerical simulations.These values were obtained by either experimental measurements or theorical calculations [27,28,30].Using Equation ( 6), the corresponding parameters can be calculated, i.e., η = 0.47 and 1.36, respectively, for Inconel 690 alloy and martensitic steels.The calculated η of Inconel 690 alloy was very close to the one obtained by first-principle calculations for Ni-∑5(012) grain boundary as shown in Figure 1.However, the martensitic steels had a higher η (>1), due to a relatively higher Γ max (8.08 × 10 −5 mol/m 2 ) [30].

Parameter
Alloy 690 Martensitic Steel

36
Note: γ FS -surface energy, γ int -interface energy, 2γ 0 int -cohesive energy, ∆h 0 int -enthalpies of hydrogen segregation to an interface, ∆s 0 int -entropies of hydrogen segregation to an interface, ∆g 0 int -Gibbs free energy excess of hydrogen segregation to an interface, ∆h 0 FS -enthalpies of hydrogen absorption at a free surface, ∆s 0 FS entropies of hydrogen absorption at a free surface, T-temperature, Γ max -saturated hydrogen coverage, E b -hydrogen trap binding energy, η-calculated by Equation (6).
Hydrostatic stress ahead of a mode I crack tip can be given by σ H = nσ ys .For a blunt crack tip, the classical continuum fracture mechanics theory (J 2 ) predicts n = 3 ∼ 5, depending on the corresponding work-hardening behavior.For an atomistically-sharp intergranular crack tip, n is calculated between 4 to 7 in a dislocation-free zone ahead of crack tip [39].Gerberich et al. [40] suggested a higher n (8~13) through the computational simulations based on discretized dislocation model.In this work, we took n = 4, since a higher Γ max had already been used.On the other hand, it was proved that n ≈ 4 predicted well in finite element simulations [41] and other models [12,16,42] compared with experimental data.Moreover, according to HELP theory, hydrogen can reduce the initial yield stress [14,15], which is equivalent to the drop of n.
Experimentally measured CSIF of Inconel 690, 625 and 718 alloy as a function of different initial hydrogen concentrations have been reported [43,44].Using the parameters recommended by Liang et al. [27] for Inconel 690 alloy and crack tip hydrostatic stress σ H = 4σ ys , the CSIF values of these three nickel-based alloys were also predicted directly by Equation (9). Figure 2 shows the comparison between the predicted CSIF and the measured results reproduced from Refs.[43,44].It is obvious that the predicted values using Equation ( 9) are consistent with the experimental results of the three nickel-based alloys, validating the applicability of the model.Note:  -surface energy,  -interface energy, 2 -cohesive energy, ∆ℎ -enthalpies of hydrogen segregation to an interface, ∆ -entropies of hydrogen segregation to an interface, ∆ -Gibbs free energy excess of hydrogen segregation to an interface, ∆ℎ -enthalpies of hydrogen absorption at a free surface, ∆ entropies of hydrogen absorption at a free surface, temperature,  -saturated hydrogen coverage,  -hydrogen trap binding energy, -calculated by Equation ( 6).  9) and experimental values of critical stress intensity factor (CSIF) for three nickel-based alloys reproduced from [43,44] (Inconel 690, 718 and 625) as a function of initial hydrogen concentration (the values of Inconel 690 alloy in the figure were obtained by  = /(1 −  ) / ,  is the J-integral given in the reference [43], "solution annealed" means the materials were solution annealed at 1100 °C for 4 h followed by air cooling, and "aged" means aged at 718 °C for 10 h).In order to keep consistent with the references,  is in weight ratio in the figure, which should be tranformed to atomic ratio when using Equation (9).9) and experimental values of critical stress intensity factor (CSIF) for three nickel-based alloys reproduced from [43,44] (Inconel 690, 718 and 625) as a function of initial hydrogen concentration (the values of Inconel 690 alloy in the figure were obtained by K = JE/ 1 − ν 2 1/2 , J is the J-integral given in the reference [43], "solution annealed" means the materials were solution annealed at 1100 • C for 4 h followed by air cooling, and "aged" means aged at 718 • C for 10 h).In order to keep consistent with the references, C 0 is in weight ratio in the figure, which should be transformed to atomic ratio when using Equation ( 9).
As aforementioned, many previous models introduced the empirical coefficient which cannot be determined in advance, therefore a fitting method was often used.Clearly, fitting is unnecessary for the model of Equation ( 9) as each parameter in the model corresponds to a specific physical variable which can be determined in advance by experimental measurements or theorical calculations.However, in order to avoid the complex processing of all parameters for a given metal, Equation ( 9) was also validated by fitting a series of experimental data [19,21,[45][46][47][48][49] reported for different materials by selecting E b and η as the fitting parameters.The other parameters are set as constants using the values reported in the literature [1,2,[18][19][20][21]27,50], as listed in Table 2. Figure 3 shows the fitting curves and obtained E b and η values.It is shown that all experimental data can be fitted very well by Equation ( 9), which further verifies that Equation ( 9) is applicable and accurate even though some parameters are assigned as constant to reduce the number of variables to be fitted.It was also found that the η and E b values obtained by fittings were in the common range of these two parameters (η ranges from 0.2 to 1 and E b ranges from 10 to 25 kJ/mol for these materials).Moreover, Figure 3e,f indicate that η increased with decreasing temperature, in line with the experimental observation that decreasing temperature enhances HAC [48,49].
very well by Equation ( 9), which further verifies that Equation ( 9) is applicable and accurate even though some parameters are assigned as constant to reduce the number of variables to be fitted.It was also found that the  and  values obtained by fittings were in the common range of these two parameters ( ranges from 0.2 to 1 and  ranges from 10 to 25 kJ/mol for these materials).Moreover, Figure 3e,f indicate that  increased with decreasing temperature, in line with the experimental observation that decreasing temperature enhances HAC [48,49].(d) A516 carbon steel [46] and X42 pipeline steel [47]; (e) AISI 4130 high strength steel [48] with temperature effect; (f) AISI 4340 high strength steel [49] with temperature effect.In (a,b,c),  is in weight ratio, which should be transformed to atomic ratio when using Equation (9).In (d,e,f),  =  √, where  = 180 •√ exp . / [12] is the Sieverts solubility constant,  is hydrogen gas pressure.
It has been suggested that, for precipitation-hardened alloys containing cleavage planes or grain boundaries decorated by precipitates, the precipitates are responsible for the decohesion of the interfaces, as they constitute the sites of decohesion due to increase of interfacial stress caused by impingements of dislocation pile-ups onto the precipitates  [21]; (b) AERMET 100 ultrahigh-strength steel [19]; (c) AISI 4340 high strength steel [45]; (d) A516 carbon steel [46] and X42 pipeline steel [47]; (e) AISI 4130 high strength steel [48] with temperature effect; (f) AISI 4340 high strength steel [49] with temperature effect.In (a-c), C 0 is in weight ratio, which should be transformed to atomic ratio when using Equation (9).In (d-f), C 0 = S 0 √ P, where

RT
[12] is the Sieverts solubility constant, P is hydrogen gas pressure.
It has been suggested that, for precipitation-hardened alloys containing cleavage planes or grain boundaries decorated by precipitates, the precipitates are responsible for the decohesion of the interfaces, as they constitute the sites of decohesion due to increase of interfacial stress caused by impingements of dislocation pile-ups onto the precipitates and the reduction of interface cohesive strength caused by the presence of hydrogen [27][28][29][30].Liang and Sofronis [27] have reported that carbide-matrix interface decohesion can produce microcracks propagating along a grain boundary.Ai et al. [51] identified an E b value of 10.2 ± 4.6 kJ/mol for the Ni 3 Al precipitates, which are major reversible trap sites decorated along various interfaces in aged Monel K-500 alloy, while Li et al. [52] obtained the value of 11.5 kJ/mol for the M 2 C precipitates in AERMET 100 steel.Clearly, these two experimental E b values are very close to the ones obtained by the fittings, i.e., 14.2 kJ/mol for Monel K-500 alloy (Figure 3a) and 10.0 kJ/mol for AERMET 100 steel (Figure 3b).This indicates that these E b values obtained by fittings in Figure 3 are credible.

Discussion
Assuming that microscopic fracture parameters such as local fracture stress/energy decreased linearly as local hydrogen concentration increased, Oriani and Josephic [7], Akhurst and Baker [8], Lee and Unger [9], Unger [10] and Wang et al. [12] proposed their models to evaluate the macroscopic fracture properties such as the crack growth CSIF.However, the linearly decreasing coefficients they introduced cannot be known in advance due to the lack of physical origins.In this paper, we proposed a physically based simple model by combining the MSCG model and TDD theory.The model can provide a direct prediction for the macroscale CSIF of materials with hydrogen influence, as all parameters in the model have clear physical origins and are measurable or calculable.
In Equation ( 9), η is a unitless parameter and C 0 is in atom ratio.Analyzing Equation (10), it can be further found that when the initial hydrogen concentration C 0 is small, Equation (10) can be simplified as θ H = C 0 exp(σ h V H /RT) exp(E b /RT), at this point, θ H is proportional to C 0 .However, when C 0 increases, θ H approaches gradually to 1.
In Gerberich's model, Equation (2), the unit of coefficient α depends on the unit of C H .There is no limit on the selection of the unit of C H . Taking account of the hydrostatic stress effect and hydrogen trapping effect on hydrogen concentration, C H = C 0 exp(σ h V H /RT) exp(E b /RT).However, if we divide Equation (2) by Equation (1), we can obtain the following relationship, as: If the hydrogen concentration is low, i.e., C H < 1, we have (αC H ) 2 2αk IG C H .In fact, it has been reported that a few ppm of hydrogen concentration is enough to induce HAF for high-strength steels, particularly for ultrahigh-strength steels.Thus, ignoring (αθ H ) 2 , Equation (11) can be approximated as: Clearly, if the hydrogen concentration C H in Equation (12) also uses the unit of atom ratio, C H is approximately equal to θ H in amount at the low level of hydrogen concentration.At this point, replacing the local hydrogen concentration term C H with θ H , Equation (12) shows the same form with Equation (8) and It is indicated that model Equations ( 2) and ( 8) are equivalent to the condition of low hydrogen concentration level and the empirical parameter α postulated by Gerberich et al. [1,2]  physically based, which can be measured or calculated.In fact, Equation (8) can be approximated by its first order of Taylor series expansion, given as: where λ = −ηk 2 IG α σ ys is a constant for a given material.As known, σ H /σ C = K IH /K IC , where σ H and σ C are the local fracture stress with and without hydrogen influence, respectively.The analysis explains why the models based on the assumption that fracture stress is linearly degraded by concentration of hydrogen can lead to good fitting of experimental data.In other words, the TDD theory can provide a physical basis for this assumption and the model is derived from it.
However, if the hydrogen concentration level is gradually increased, θ H will approach to 1. Therefore, there is a limit value for hydrogen-influenced K IH : lim This is consistent with experimental trends shown in Figure 3a-c,e,f.It is shown that when hydrogen concentration or hydrogen pressure is increased to a high level, K IH is no longer decreased and remains at a minimal value, K IHmin .In fact, this seems to be a common phenomenon.For example, Barthelemy [53] presented the effect of hydrogen partial pressure on rupture pressure in disk rupture tests where hydrogen pressure was increased to a specific level and the helium gas was added until the rupture of the disk.It is shown that the embrittlement index defined as the ratio of helium pressure to hydrogen pressure (P He /P H2 ) increased quickly at low hydrogen partial pressure (before reaching a constant level near 25 bar).However, further increase of hydrogen pressure until 200 bar had no effect on the embrittlement index.Elices and Gutierrez-Solana [47] performed tensile tests on double-notched X42 specimens and showed a ~40% reduction in area loss at hydrogen pressure near 7 MPa.However, increasing the hydrogen pressure to 35 MPa only resulted in a ~50% reduction in area loss.Although the reduction in area loss appears to show a continual increase, it is possible to obtain a critical pressure.Gerberich's model, Equation (2), also shows that K IH decreases with increasing hydrogen concentration; however it will reach to its limit when k IG αc H and the limit is 0. This clearly make no sense in physics.Therefore, the model, Equation (9), proposed in this paper is more consistent with the experimental data.
It should be noted that, according to Equation (15), the K IHmin values are 15, 10, 11, 60 and 8 MPa• √ m, respectively, for Monel K-500 alloy, AERMET 100 ultrahigh-strength steel, AISI 4340 high-strength steel, A516 carbon steel and X42 pipeline steel.Different materials have different K IHmin values.The K IHmin value can be used as a criterion to quantitate and classify the HAF susceptibility of a material.A material with a higher K IHmin is preferred for hydrogen service.Moreover, for a conservative design or safe assessment of a component exposed to a hydrogen-containing environment.It can be believed no matter how much hydrogen is introduced into the component, HAF will not occur when K I < K IHmin .
In addition, as aforementioned, the model does not consider the role of plastic work for simplicity.In fact, plastic work can be involved by using the expression proposed by Jokl et al. [54], which describes the plastic work γ P as a function of 2γ 0 int by γ P = A 2γ 0 int q , where A and q are constants to be determined.At this point, as plastic work is much larger than the reversible work of decohesion, that is, the latter can be ignored (γ eff = 2γ 0 int + γ P ≈ γ P ), the k IGH can be expressed as k IGH = k IG (1 − ηθ H ) q and thus the model Equation ( 9) becomes Using Equation ( 16) to fit the experimental data of Monel K-500 alloy with different q values, one can see that the data can also be fitted very well with almost no gap between fitting curves (q = 1, 2, 4 and 6).With the increase of q, the obtained η and E b decreases but E b tends to plateau after q > 6, as shown in Figure 4. Therefore, whether plastic work is considered will not affect the applicability of the model.However, to understand the role of plastic work and quantify its contribution (for example, determining the value of q for different materials), more experimental and theoretical work is needed.

Conclusions
This paper proposed a simple physically based model to predict the crack growth critical stress intensity factor (CSIF) with hydrogen influence by combining the microscopically shielded Griffith criterion (MSGC) model and thermodynamics decohesion (TDD) theory.The model predicts the hydrogen-influenced CSIF as a function of the CSIF without hydrogen influence, initial hydrogen concentration or hydrogen solubility, hydrogen trap binding energy and crack tip stress level.The prediction from the model shows good agreement with experimental data reported in the literature.The advantage of the model is that all parameters in the model have physical origins and can be measured or calculated in advance.It can promote the application of stress intensity factor-based crack growth criteria in the design and safety assessment of a metallic component exposing to hydrogen.

Figure 2 .
Figure 2. Comparisons between predicted values by Equation (9) and experimental values of critical stress intensity factor (CSIF) for three nickel-based alloys reproduced from[43,44] (Inconel 690, 718 and 625) as a function of initial hydrogen concentration (the values of Inconel 690 alloy in the figure were obtained by  = /(1 −  ) / ,  is the J-integral given in the reference[43], "solution annealed" means the materials were solution annealed at 1100 °C for 4 h followed by air cooling, and "aged" means aged at 718 °C for 10 h).In order to keep consistent with the references,  is in weight ratio in the figure, which should be tranformed to atomic ratio when using Equation (9).

Figure 2 .
Figure 2. Comparisons between predicted values by Equation (9) and experimental values of critical stress intensity factor (CSIF) for three nickel-based alloys reproduced from [43,44] (Inconel 690, 718 and 625) as a function of initial hydrogen concentration (the values of Inconel 690 alloy in the figure

Figure 4 .
Figure 4.The fitting curves based on Equation (10) for critical stress intensity factor (CSIF) of Monel K-500 alloy using different  values.
can be replaced by ηk IG /2.α has no clear physical origin, but η and k IG are