A cohesive zone model to simulate the hydrogen embrittlement effect on a high-strength steel

The present work aims to model the fracture mechanical behavior of a high-strength low carbon steel, AISI 4130 operating in hydrogen contaminated environment. The study deals with the development of 2D finite element cohesive zone model (CZM) reproducing a toughness test. Along the symmetry plane over the crack path of a C(T) specimen a zero thickness layer of cohesive elements are implemented in order to simulate the crack propagation. The main feature of this kind of model is the definition of a traction-separation law (TSL) that reproduces the constitutive response of the material inside to the cohesive elements. Starting from a TSL calibrated on hydrogen non-contaminated material, the embrittlement effect is simulated by reducing the cohesive energy according to the total hydrogen content including the lattice sites (NILS) and the trapped amount. In this perspective, the proposed model consists of three steps of simulations. First step evaluates the hydrostatic pressure. It drives the initial hydrogen concentration assigned in the second step, a mass diffusion analysis, defining in this way the contribution of hydrogen moving across the interstitial lattice sites. The final stress analysis, allows getting the total hydrogen content, including the trapped amount, and evaluating the new crack initiation and propagation due to the hydrogen presence. The model is implemented in both plane strain and plane stress configurations; results are compared in the discussion. From the analyses, it resulted that hydrogen is located only into lattice sites and not in traps, and that the considered steel experiences a high hydrogen susceptibility. By the proposed procedure, the developed numerical model seems a reliable and quick tool able to estimate the mechanical behavior of steels in presence of hydrogen.


INTRODUCTION
ydrogen embrittlement phenomenon is an issue known since several years in engineering field.Different structural steels and alloys show sensitivity to hydrogen.In particular, when atomic hydrogen gets in contact with these materials they experience a drastically decrease of the mechanical properties that can result in failure H of components.Hydrogen embrittlement phenomenon interests different fields such as mechanical, structural and energetic.For some industrial environment, this problem is widely recognized and studied in literature.For instance, oil&gas industry [1] in which atomic hydrogen is released as product of chemical reactions in environments where the infrastructures operate, pressure vessels for hydrogen storage and transportation [2] and lately even energy devices that use hydrogen as alternative energy carrier.However, other applications in which the presence of hydrogen is less evident have been pointed out recently thanks to the ongoing research on this topic.An example is reported in [3], dealing with wind turbine gearbox bearings, where the hydrogen presence has a deleterious effect in combination with rolling contact fatigue.In this case, it is suggested that hydrogen comes from decomposition of lubricating oil [4] or from water contamination.Scientific literature also reports some failures occurred in threaded fasteners, as in [5] where the possible sources of hydrogen are related to thermal treatments.However, independently on the source that generates atomic hydrogen, the most crucial phase is the diffusion process of hydrogen through the material lattice.According to [6], usually the concentration of hydrogen is split into two parts: the content of hydrogen in the interstitial lattice sites (NILS) driven by hydrostatic pressure, and the amount accumulated in correspondence of the so-called trap sites.In turn, these can be divided in reversible and irreversible based on the hydrogen binding energy (potential energy at microscopic scale).Reversible traps, or low binding energy traps, are mainly related to dislocations and plastic flow.In fact, in [6] the authors showed that plastic strain and hydrogen concentration in reversible traps have similar trends in front of a crack tip.The presence of a crack in a component induces hydrogen atoms to move from the free surface towards the tip.Indeed, crack initiation and propagation are deeply influenced by hydrogen presence and diffusion.In terms of diffusion coefficient the motion of hydrogen through NILS is represented by an ideal lattice diffusivity, D L .The diffusion can be limited or increased by traps and in these circumstances, a trap-affected or apparent diffusivity, D H , is considered.Hydrogen embrittlement is mostly governed by this second diffusion coefficient.Traps effect is not univocal [7,8].Indeed, literature reports that hydrogen in reversible traps is in equilibrium with the one in NILS, and it is an "obstacle" to its transport, thus D H <D L .Anyway, a large population of reversible traps will have a high probability of hydrogen release at room temperature and it can provide a reservoir of mobile (diffusible) hydrogen to supply the crack-tip fracture-process zone and thus increase hydrogen embrittlement.Quenched and tempered high strength steels have microstructural features that increase the strength, but also provide for many hydrogen-trapping sites [9].Contrary, irreversible traps, often saturated even at low hydrogen concentrations, no longer interact with the dissolved hydrogen in the lattice, and D H tends to D L [10].However, a homogeneous distribution of high-energy traps can prevent that hydrogen localize at low energy sites and shield its segregation at the crack tip.Regarding embrittlement mechanisms, several authors proposed possible theories but the most acceptable in literature are three: HEDE (Hydrogen Enhanced DEcohesion, [11]), HELP (Hydrogen Enhanced Local Plasticity, [12]) and hydrides formation and cleavage.These physical models try to explain the embrittlement phenomenon starting from nano and microscopic considerations, and propose that the damage located at crack tip is proportional to hydrogen concentration into the steel [13].Therefore, it is clear that hydrogen embrittlement is a complex mechanism from both mechanical and physical-chemical standpoints.Based, on these considerations, it is hence really important to develop a tool able to assess the sensitivity of steels to hydrogen embrittlement.Literature reports several numerical models that try to simulate properly the embrittlement process considering interacting factors and the involved parameters mainly microstructure dependent.Among these numerical models, cohesive zone models (CZM) seem the most promising to reproduce at macro-scale the effect of hydrogen on mechanical properties of steels/alloys.These finite elements models assemble hydrogen diffusion process together with fracture mechanical behavior through the implementation of cohesive elements.A traction separation law (TSL), whose formulation includes parameters able to describe crack initiation and propagation, governs these particular elements.Therefore, the effect of embrittlement is simulated simply by reduction of the cohesive energy based on the hydrogen content, omitting the micro-mechanisms.The major contributions to this approach come from [14,15,16].Even the model presented in the current manuscript belongs to this category.A complex finite element model including three steps of simulations is developed to reproduce the mechanical behavior of a steel, named AISI 4130, during a toughness tests in hydrogen-contaminated environment.This type of test was carried out during an extensive campaign of experimental tests, being suitable to describe the embrittling effect in steels due to hydrogen.Experimental fracture data are used to calibrate the mechanical behavior of cohesive elements describing crack propagation.Considerations about the implementation of 2D plane stress and plane strain models are proposed.Hydrogen concentration estimated at the crack tip is discussed.

COHESIVE ZONE MODEL SETUP
tarting from experimental results obtained by toughness tests on AISI 4130 steel in hydrogen uncharged and charged conditions [17], a finite element cohesive zone model is built to reproduce the mechanical response of the material.The model, developed using Abaqus software, simulates a 2D C(T) specimen.Exploiting the symmetry of the geometry, only half part of the specimen is considered and a line of cohesive elements is implemented all over the crack propagation path.These zero thickness elements are connected to the continuum elements on the symmetry plane.In this region, the size of the mesh reaches the smallest dimension, around 10 µm, to ensure proper resolution of results.This mesh size is the result of a convergence study.Mechanical properties of continuum elements are assigned according to the data obtained by experimental tensile tests carried out on AISI 4130 steel.Fig. 1 depicts the stress-displacement curve of the base material together with the specific values of Young's modulus, Poisson's ratio and yielding stress.Instead, peculiarity of cohesive zone model is the definition of a phenomenological law, identified as traction separation law (TSL) describing the constitutive behavior of the material inside the cohesive elements.The cohesive behavior is expressed by a stress-displacement function that can present different shapes according to the crack propagation response of the material.Usually, for steels or metals the most noted are Gaussian or trapezoidal shapes in which it is possible to split elastic, plastic and failure parts.For the current case, the TSL presents a trapezoidal shape defined by the four parameters shown in Fig. 1: σ 0, δ 0 , δ N , δ F .Contrary to other authors that customized the user defined elements (UEL) subroutine to define the traction separation law of cohesive elements [18,19], this work presents an easier approach.In fact, the TSL is designed exploiting the traction-separation model available in Abaqus, which assumes an initial linear elastic behavior followed by an initiation and evolution of damage.

S
The damage initiation criterion, by which the process of degradation starts, is the maximum nominal strain criterion according to the definition given in [20].It implies that the damage initiates when the nominal strain ε evaluated at each ith increment of the analysis, reaches the user defined maximum allowable strain value ε 0 .After that, to describe the rate of the degradation process a damage variable D is used.It represents the overall damage in the material and it is related to the stress components as it follows [20]: where σi is the stress tensor computed in the current increment considering absence of damage, and i  is the stress predicted by the elastic traction-separation in the same increment.The values of D run monotonically over a range from zero to one: zero represents the starting condition of damage; one describes a situation of full damage in which the element cannot hold up the applied load anymore.
In order to assign the TSL shape, a tabular function of the effective displacement δ i beyond damage initiation is used.The formulations of damage D as a function of the effective displacement δ i is reported below: Failure part: The values of these parameters are chosen in order to fit the experimental crack propagation data.Therefore, a calibration process of TSL for this specific steel without considering the hydrogen effect is necessary as starting point to develop the more complex model including the presence of hydrogen.

Calibration process of TSL parameters for AISI 4130 steel
Generally, the cohesive law can be calibrated by a comparison with experimental crack propagation data or it can be theoretically modeled using a numerical method that resolves the fracture process.In this case, we selected the experimental approach.A calibration procedure is developed to estimate the cohesive parameters more suitable to represent the characteristic fracture behavior of AISI 4130, in hydrogen non-contaminated conditions.This calibration compares the trends of the applied force versus the vertical load line displacement (F-VLL) represented in Fig. 2, obtained from experimental toughness tests and from the numerical cohesive model simulating the same loading condition.Force, F, is evaluated at the center of the loading pin, and it corresponds to the load cell measurement.Displacement, V LL , is measured at the node where the clip-on-gauge is experimentally located.In this way, experimental data can be directly compared with the numerical outputs.

Simulation of hydrogen effect
The decreasing of the mechanical performances of a steel due to the embrittlement effect are strictly related to hydrogen content.Therefore, to simulate correctly this phenomenon it is necessary to take into account properly the total amount of hydrogen present into the material.
As already highlighted in the introduction, the concentration of hydrogen inside the steel lattice consists of two contributions: the interstitial lattice sites content and the trapped hydrogen content.To estimate both these contributions, a three-steps model is implemented.Once the total hydrogen concentration is calculated, it is then used to reduce the cohesive law, simulating hydrogen embrittlement effect.The first step is a stress analysis used to calculate the hydrostatic stress field into the specimen.The second step, instead, consists of a mass diffusion analysis in which a redistribution of an initial hydrogen concentration, selected according to the experimental values equal to 1.5 ppm at the free surfaces [17], is computed based on the hydrostatic stress field previously calculated.The equation used by Abaqus to solve the diffusion analysis is the following, with the omission of the temperature term that is not considered: where J * is the flux, s is the solubility, D is the diffusivity of the material, φ is the normalized concentration over the solubility, k P is the term that induces diffusion due to the hydrostatic effect and p is the hydrostatic stress.Table 1: Values of parameters adopted in the three simulations of the cohesive zone model.
Finally, the last step is a stress analysis with cohesive elements implementation to simulate the fracture behavior.In this simulation, also the trap site concentration is computed.
As indicated in the introduction, influence of traps on hydrogen embrittlement is certainly important; however, the characterization of trap type, density, binding energy and occupancy is a demanding task that has to be performed at smaller scales.In order to simplify this problem and to include this contribution in the finite elements codes, we consider here only hydrogen low energy traps and we assume that these are related to the plastic strain.
In fact, starting from results presented in [6], the work by [15] proposed the following relation between the trap content and the plastic strain:   49.0 0.1 This relation was proposed for high and low strength steels, from physical-chemical analytical relations and numerical models [6].It is implemented in the third analysis using the plastic strain calculated during the simulation.The total bulk hydrogen concentration, C, resulted from the sum of interstitial lattice sites CL and the traps sites CT is related to the surface concentration in the cohesive zone through [21]: where θ is the hydrogen coverage factor, 0 b g  is the variation of Gibbs free energy, R is the gas constant and T is the temperature.In turn, this term contributes to the calculation of k factor, by which the cohesive energy is reduced simulating the embrittlement effect [15]: This scaling factor is applied on both stress and displacement of the TSL curve reducing the entire area.
All the relationships and parameters in Eqs. ( 4), ( 5) and ( 6) are calculated in the last analysis, for each increment, in the unique integration point of the continuum element.The parameter, k, is then transferred from the continuum to the corresponding cohesive elements by implementing three subroutines.These are the user subroutines to manage userdefined external databases and calculate model-independent history information (UEXTERNALDB), the user subroutine to redefine field variables at a material point (USDFLD), and the user subroutine to generate element output (UVARM).These subroutines work interactively each other in a Fortran common block.The performed experimental toughness test from which the model have been developed was in accordance with standards, and specimen dimensions ensured the definition not only of J-toughness but also of K-toughness.Therefore, at the specimen core, plane strain condition is likely to occur.On the contrary, at the external surfaces, the stress state would be more similar to plane stress condition.Since the model we are implementing is bi-dimensional, limitations in the analysis should be taken into account.For this reason, the model has been developed both in plane strain and plane stress conditions.Apart from the different type of continuum elements used to describe the material (CPE4 for plane strain and CPS4 for plane stress analysis), the model geometry, the applied displacement and boundary conditions and the followed calibration procedure are coincident.Even if the modelling steps are the same, the plane strain and plane stress models give different results, which will be critically discussed in the following section.

RESULTS AND DISCUSSION
irst, a TSL calibration for cohesive elements is developed to find out the cohesive parameters that replicate the behavior of the material in absence of hydrogen for both plane stress and plane strain configurations model.Tab. 2 shows the resulting TSL parameters.The TSLs have a very small flat part, i.e. (δN -δ0); on the contrary, they present a long trend describing the damage to failure, i.e. (δF -δN).This means that the TSL is more similar to a triangle than a trapezoid.The TSL plateau stress, σ 0 is between 1.75 and 1.85 the yielding stress.The area below the TSL plot corresponds to the separation work needed for crack initiation and propagation.Values of this area are deeply different between plane strain and plain stress cases, accounting for the stress combination at the separation surface.It should be mentioned that data in Tab. 2 are results of a trial and error procedure, and they correspond to the best combination of TSL parameters able to fit experimental data.Table 2: Characteristic values for TSL parameters for the hydrogen-free model, for plane strain and plane stress models.Displacement and stress notation is referred to Fig. 1; the area is below the TSL curve.

Hydrogen
Now, considering the model that predicts the response of the material in presence of hydrogen, the first step of analysis evaluates the hydrostatic stress field.This is imported as predefined field in the second mass diffusion analysis.Here, the hydrostatic stress drives the flux of hydrogen computing the interstitial lattice sites content, CL.Fig. 3 shows CL F concentration field around the crack tip, evaluated during the second step of the analysis for plane strain model.In accordance with literature [14], the peak of hydrostatic stress and then of hydrogen NILS concentration is not located at the crack tip, but slightly forward.In our simulations, this peak is located at the end of the first element ahead the tip.In this node, the model estimates a quite high CL concentration, equal to 4.43 ppm.It is around three times the initial imposed concentration at free edges (1.5 ppm).For plane stress model, similar contour is found, with a peak value for hydrogen concentration of 2.03 ppm.The last step calculates the current plastic strain, and consequently the contribution of the trapped sites, CT.Then, based on the total hydrogen concentration, NILS and traps, the decreasing factor k is applied to the previous TSL to reduce the area below the curve, thus the energy adsorbed by the cohesive elements.Thanks to this factor, it is possible to simulate a lower resistance to crack propagation and to evidence the embrittlement of the steel.
In order to get F-VLL curves similar to experimental data, it results that k has to be applied to both cohesive stress and displacement of TSL.In particular, for plane strain model in presence of hydrogen, we found that the reduced TSL stress and displacement are respectively σ0H=k•σ0 and δH=k 3/2 •δ, where δ is the generic displacement of the TSL curve.While the reduction in stress is in agreement with literature results [15,16], the further reduction in displacement was quite unexpected.However, it can be justified by the high hydrogen sensitivity observed during the experimental fracture toughness tests for the current steel.In fact, it showed a more brittle behavior compared to other steels tested in the past [22].For plane stress model, the same ratio was derived for the reduced TSL stress σ0H=k•σ0, while the reduced displacement is slightly lower: δ H =k•δ. Fig. 4 shows contours of the decreasing factor k around the crack tip for plane strain model.The maximum is located in front of the crack tip and it is equal to 0.26.On the tip, it is slightly lower and equal to 0.29.Instead, the maximum value of k factor for plane stress model is 0.36.This means that the embrittlement effect is more evident in case of plane strain than plane stress.Moreover, both models agree in the prediction of absence of plastic strain measured at the crack tip during the last step.Therefore, from the numerical simulations, it seems that no hydrogen is located in traps, i.e.C T =0, and the only contribution for k factor calculation is provided by the CL hydrogen amount.This means that the embrittlement is only due to hydrogen in the lattice and suggests a very brittle behavior of AISI 4130 in presence of hydrogen.As reported in literature [7,8], AISI 4130 is supposed to have many trap sites, mainly due to its martensitic structure, resulting in high strength, rather than to its plastic strain.
Although, the study should be further deepened, these numerical finding can be supported by the first experimental observations of the fracture surface of the specimens tested in hydrogen environment.In fact, for most of the steels, on the fracture surfaces vast brittle regions and smaller ductile areas appeared, likely related to the local plastic strain and therefore to hydrogen in traps [22].Instead, in the case of AISI 4130, from the images collected after toughness tests, most of the failure surfaces showed brittle features, as in Fig. 5.Some cracks also propagated from lateral groves.Surely, to confirm the fully brittle behavior of this steel in presence of hydrogen more data would be necessary.In literature, there is no analytical law to calculate the hydrogen concentration in traps CT, to be compared with this numerical result.Moreover, the diffusion coefficient of this steel should be experimentally measured, including its effective and apparent contributions (DL and DH).If these two diffusivity values were available, their ratio would allow an estimation of CT, according to [7].
As a final comment to our result, it is even possible that the law expressed in Eq. 4 is not completely able to describe the mechanical behavior of AISI 4130, even if it included a wide class of materials.Fig. 6 compares force displacement (F-VLL) curves from experimental tests and numerical simulations (plane strain and plane stress).The plot shows the mechanical responses of both models in hydrogen uncharged and charged conditions.It is evident the difference between the two behaviors without and with hydrogen.The macro-mechanical response of the specimen during toughness test is deeply influenced by hydrogen.The curve is decreased not only in its load peak, but also in the maximum displacement and in the damage history, as confirmed from TSL data in Tab. 2. By observing Fig. 6, it is possible to note that in plane stress condition the calibration process does not allow a good overlapping with the experimental results.Despite many simulations, it appeared that the combination of TSL parameters for this model never resulted in a good fitting with the experimental data.
In terms of fracture toughness J value, that represents the area below the F-VLL curve, the value obtained by the calibration process simulating the test without hydrogen is a bit smaller than the experimental one (experimental fracture toughness without hydrogen = 215 N/mm).Then, this tiny discrepancy is reversed comparing experimental and numerical curves related to the hydrogen environment conditions (experimental fracture toughness with hydrogen = 22 N/mm).On the contrary, the plane strain model better fits experimental data in both hydrogen and no-hydrogen conditions.Therefore, between the two developed models, the 2D plane strain one seems to represent better the experimental evidence.However, as a future development, a 3D numerical model would give more information about the embrittlement depending on the local stress condition.

CONCLUSIONS
he current work presents a cohesive zone model reproducing a toughness test in presence of hydrogen for both plane stress and plane strain configurations.The following conclusions can be drawn:  firstly a calibration of traction-separation law (TSL) that defines the constitutive properties of cohesive elements is developed for the steel AISI 4130;  a model running on three steps of simulations calculates the total hydrogen concentration present into the specimen including the interstitial lattice sites and the trapped amount.Based on the total content of hydrogen a scale factor is computed and used to reduce the TSL area previously calibrated in order to simulate the embrittlement effect of hydrogen;  for the present steel, AISI 4130, the considered relation to calculate trapped hydrogen concentration CT depending on plastic strain seems not appropriate.By the implemented equations, the numerical model is not able to estimate the hydrogen concentration in reversible traps.Therefore, the model computes hydrogen embrittlement only as a function of hydrogen in interstitial lattice sites.This is a consequence of a brittle behavior of the steel, according to the experimental test; T  between the plane strain and plane stress modelling configurations the more suitable to represent the experimental results is the plane strain;  the developed model seems a tool able to provide a proper description of the behavior of the material in presence of hydrogen, but still needs improvements to clarify the influence of trap sites to the embrittlement phenomenon.

Figure 1 :
Figure 1: Experimental tensile curve and values of mechanical parameters used for continuum elements.TSL shape as stressdisplacement trend.

Figure 2 :
Figure 2: Scheme of the specimen and the F -V LL quantities used for the calibration process.
Tab. 1 summarizes the values of the parameters used in the model.It should be pointed out that the scientific literature only reports experimental values of the apparent diffusion coefficient DH for the current steel.Thus, this was the value implemented into the model.The concentration obtained from this analysis refers only to the "free" hydrogen present in the interstitial lattice sites, C L .ppm mm N -1/2 1.62•10 -5 mm 2 /s 30 kJ/mol 8.314 J K -1 mol -1

Figure 3 :
Figure 3: Contours of C L concentration in lattice interstitial sites at crack tip region, plane strain model.

Figure 5 :Figure 6 :
Figure 5: Fracture surface of a specimen after the toughness test.The specimen was hydrogen pre-charged before the test.Typical brittle surface is evidenced by facets (black arrow), intergranular cracks (white arrow) along grain borders and secondary cracks (yellow arrow).