Analysis of Hydrogen-Assisted Brittle Fracture Using Phase-Field Damage Modelling Considering Hydrogen Enhanced Decohesion Mechanism

: This study proposes a hydrogen-assisted fracture analysis methodology considering associated deformation and hydrogen transport inside a phase-ﬁeld-based formulation. First, the hydrogen transport around a crack tip is calculated, and then the effect of hydrogen enhanced decohesion (HEDE) is modeled by deﬁning the critical energy release rate as a function of hydrogen concentration. The proposed method is based on a coupled hydrogen mechanical damage under phase-ﬁeld and implemented through a user subroutine in ABAQUS software. The test using compact tension (CT) sample is investigated numerically to study the hydrogen embrittlement on 45CrNiMoVA steel. Experimentally, the microstructural fracture presents a mixed brittle fracture mode, consisting of quasi-cleavage (QC) and intergranular (IG) fracture with hydrogen. This fracture mode is consistent with the suggested HEDE mechanism in the model. The simulation results show that hydrogen accumulates at the crack tip where positive hydrostatic stress is located. Moreover, the model estimates the initial hydrogen concentration through iterations. The simulated load-line displacement curves show good agreement with the experimental plots, demonstrating the predictive capabilities of the presented scheme.


Introduction
During the material forming processes, the atomic hydrogen tries to penetrate the metal and causes degradation of mechanical properties.Due to the degradation of mechanical properties, the metal may lose its strength and ductility, and an embrittlement fracture will occur.This phenomenon is called hydrogen embrittlement (HE).According to reported studies, Johnson [1] was the first who observed the phenomenon of hydrogen embrittlement in steel in 1875 and concluded that hydrogen may lead to significant deterioration of mechanical properties of steel and could produce brittle fracture without warning.Since then, the hydrogen-assisted fracture properties of the metals have attracted much attention.The related experimental research shows that the fracture micro-mechanisms transformed from dimples ductile fracture (microvoid coalescence (MVC)) without hydrogen to brittle fracture (quasi-cleavage (QC) or intergranular (IG) fracture + transgranular (TG) fracture mechanism) with hydrogen [2,3].For various hydrogen-induced cracking phenomena, there is no unified theory that explains the mechanism underlying HE.To date, there are two generally accepted HE mechanisms: hydrogen-enhanced localized plasticity (HELP) and hydrogen-enhanced decohesion (HEDE).Birnbaum et al. [4], Sofronis et al. [5], and Liang et al. [6] introduced the HELP mechanism, which was subsequently applied to the unit cell model by Ahn et al. [7] to explore the influence of HELP on void expansion and coalescence.However, the HELP mechanism has no effect on the material's yield strength.To fully characterize hydrogen embrittlement, it must be paired with other theories such as HEDE.At first, Pfeil [8] introduced the HEDE mechanism, which was further investigated by Troiano [9], Gerberich et al. [10], and Oriani [11].The HEDE assumes that the presence of hydrogen in metal reduces the interatomic bonding force between adjacent crystallographic planes.During the fracture propagation, ductility tends to convert into brittleness when hydrogen concentration near the crack tip reaches a certain amount, lowering cohesive strength and surface energy along the grain boundary.
Experiments can be performed to determine the macroscopic characteristic curve and microscopic fracture mechanism of hydrogen embrittlement.A computational model can be developed by considering the hydrogen influence in the HEDE mechanism to predict hydrogen diffusion around a crack tip and to understand hydrogen embrittlement.Jiang and Carter [12] evaluated the fracture energy of iron that was degraded due to a hydrogen attack.They used the Born Haber thermodynamic cycle and discovered that the fracture energy reduces with an increasing hydrogen coverage.In another study, Serebrinsky [13] proposed a hydrogen degradation law based on decreasing the cohesive strength with increasing hydrogen concentration.Jembieet al. [14] and Sobhaniaragh et al. [15] simulated hydrogen-assisted fracture by using the cohesive zone model with the HEDE mechanism.Using a similar methodology, Li et al. [16] simulated hydrogen-assisted brittle fracture initiation in low alloy 30CrMo steel using the so-called three-step FE procedure coupled with hydrogen-informed CZM analysis for the HE problem.In recent studies, the phasefield model is employed by researchers to model brittle fractures.This model offers certain advantages, such as, it is easy to implement, capturing the original finite element model, and can mimic complex cracking processes (such as branching and merging of cracks).The Griffith theory [17] serves as the foundation for the phase-field model.Francfort et al. [18] used the quasi-static crack evolution of establishing the variational model that does not set a predefined crack path, hence bypassing the limitations of the Griffith theory.Bourdin et al. [19,20] investigated the viability of numerical approaches.They considered a diffuse crack region that was dominated by a phase-field variable to characterize the material's sharp fracture surface topology, which was interpolated between intact and shattered materials.Miehe et al. [21,22] proposed an operator splitting approach that updates the crack phase and displacement fields progressively to predict brittle fracture.Miehe et al. [23] further developed the phase-field modeling of ductile fracture by including the plastic behavior.Ambati et al. [24,25] applied the follow-up model by integrating the degradation function and the plastic strain state to forecast crack initiation in the plastic zone.Their study compared the numerical results to experimental data.Based on the fracture criterion of the phase-field model and its advantages to simulate the fracture mechanism, many researchers have adopted this to simulate the phenomenon of HE.In this method, hydrogen diffusion and the subsequent fracture are studied within the framework of a phase-field model.Philip et al. [26] introduced hydrogen diffusion, plastic strain gradient, and the phase-field fracture formula relating the hydrogen sensitivity.The hydrogen sensitivity is calculated by the first principles into the framework of the phase-field model to predict hydrogen-assisted fracture.Martnez-Paeda et al. [27] presented the hydrogen-assisted cracking phase-field formula using the finite element approach (FEM).Additionally, Wu et al. [28] introduced the phase-field regularized cohesive model (PF-CZM) to simulate hydrogen-induced fractures.
The presented study proposes a hydrogen-assisted fracture phase-field model that incorporates the coupled problems of hydrogen transport and the HEDE mechanism to capture hydrogen's effect on damage evolution.The remainder of the paper is structured in the following manner.Section 2 provides an overview of the phase-field model.Following that, the link between hydrogen transport and hydrogen-assisted deterioration is established.Section 3 explains the implementation of finite element analysis in ABAQUS via user subroutine UEL.Section 4 explains the hydrogen-charged (CT) specimens of 45CrNiMoVA steel to demonstrate the capability of the numerical model.This section also covers the analysis and discussion of simulation results.Finally, Section 5 draws important conclusions.Where φ = 0 represents the material in an entirely intact state and φ = 1 represents the material in a completely broken state.Assuming the bar is completely open at x = 0, the non-smooth crack topology can be modeled using the exponential function:

Hydrogen Assisted
where l is a length scale parameter and φ(x) is the phase-field value of the diffusive crack topology surface.By assuming that the length scale l → 0 , then sharp case is approached and Equation (1) decays monotonically to zero because of moving away from the crack location.
Equation ( 1) is the solution of the ordinary differential equation for x = 0 as follows: Equation ( 2) follows the Dirichlet-type boundary conditions: φ(0) = 1 and φ(±∞) = 0.This can be transformed into the Euler-Lagrange equation of the variational problem tional is then given by: The integration over the volume domain dV = Γdx gives I(φ = e −|x|/l ) = lΓ; therefore, the crack surface Γ is related to the crack length scale parameter l.Consider a discrete internal discontinuity Γ in a solid body Ω, the crack surface density Γ l (φ) is introduced by means of the regularized crack functional as: where r l (φ, φ ) denotes the one-dimensional crack surface density function.Similarly, it can be stated as follows in higher dimensions: where ∇φ represents the spatial gradient of φ.
The fracture energy dissipated by the formation of a crack can be calculated as [29]: where G c is the critical energy release rate.The total potential energy of a solid body can be used to relate the bulk energy Ψ b (u, φ) to the fracture energy Ψ s (φ) as [29]: where u is the vector of displacements and ψ 0 (ε) is the elastic strain energy density, is defined as a function of the standard strain tensor ε and the fourth-order linear elastic tensor C 0 .This can be expressed as [30]: Such that, assuming a small strain theory, the standard strain tensor ε can be calculated as: ε = 1 2 (∇u) T + ∇u .

Governing Balance Equations
By taking the variation of the total potential energy (Equation ( 7)) can be written (weak form) with respect to ε and φ renders in the absence of body forces and external tractions, which can be expressed as: where the Cauchy stress tensor is expressed as: Application of the Gauss theorem, the corresponding Eulerian equations can be written (strong form) for or any arbitrary value of the kinematic variables δu and δφ: To ensure the irreversibility when the elastomer is pressed or unloaded, the so-called history variable field H can be explained as [31]: where H t is the energy determined at any specific time increment t.Equations ( 10) and ( 11) can be combined to obtain the phase-field damage evolution equation which satisfies the Kuhn-Tucker conditions:

Governing Equation of Hydrogen Diffusion
The hydrogen diffusion analysis under conditions of stress is defined as an extension of Fick's law [32].The diffusion equation is developed by using the law of mass conservation for diffusion: where C is the diffusible hydrogen concentration in a solid body Ω, J is the hydrogen flux and n is the outward unit vector normal to the surface ∂Ω.The strong form of the balanced equation can be obtained by using the divergence theorem for any arbitrary volume: For an arbitrary, suitably continuous, scalar field, δC, the variational statement Equation ( 14) can be written as: Again, using the divergence theorem, the weak form of the equilibrium equation can be expressed as: where q = J • n is the concentration flux exiting the body across ∂Ω.
The chemical potential µ gradient instigates the diffusion process, and it can be expressed as [29]: where µ 0 is the chemical potential at a reference state, V H is the partial molar volume of hydrogen (2000 mm 3 /mol for iron [33]), T is the temperature, R is the ideal constant of gas (8314 N • mm/mol • K), θ L is the occupancy of the lattice sites, and σ H is the hydrostatic stress defined as by σ H = (σ 11 + σ 22 + σ 33 )/3.
A linear Onsager relationship is used to relate the flow to the chemical potential: where D is the ideal diffusivity.By substituting Equation (17) into Equation (18), and assuming the low occupancy (θ L 1) for most practical purposes, the diffusion flux can be expressed as: Substituting Equation (19) into Equation ( 16), the hydrogen diffusion equation can be obtained as:

Hydrogen Degradation Function
Atomic hydrogen dissolved in metal reduces the binding energy between metal atoms and leads to a reduction in fracture resistance.Alvaro et al. [34] proposed a HEDE model in which the surface energy decreases linearly with an increase in hydrogen concentration.According to this model, a quantum-mechanically informed hydrogen-dependent degradation function can be defined as: where G c (θ) represents the critical energy release rate with the HEDE effect and χ denotes the hydrogen damage factor, G c (0) is the critical energy release rate with no hydrogen influence.
Considering the hydrogen coverage θ as a function of the bulk hydrogen concentration C the framework of the Langmuir-McLean isotherm is defined as [12,35]: where C is the impurity mole fraction in units and ∆G 0 g is the Gibbs free energy.As defined in Equations ( 21) and ( 22), the hydrogen embrittlement parameter G C (θ)/G C (0) is dependent on the selected value of Gibbs free energy and hydrogen damage factor.The curves of hydrogen coverage with hydrogen concentration for various levels of Gibbs free energy range of 10 to 60 k/mol are shown in Figure 1a.Clearly, hydrogen concentration covers about four orders of magnitude for two parameters.The lower boundary and the upper boundary represent the hydrogen concentration at fracture initiation and the ultimate saturation level, respectively.Therefore, the value of the Gibbs free energy is very important and cannot be determined arbitrarily.The Gibbs free energy ∆G 0 g is assigned a value of 30 kJ/mol, representing the initiation hydrogen concentration and an ultimate saturation level ranging between 0.01 and 100 ppm in the α-Fe grain boundary [13].
pendent degradation function can be defined as: where   c G  represents the critical energy release rate with the HEDE effect and χ de- notes the hydrogen damage factor,   is the critical energy release rate with no hydrogen influence.
Considering the hydrogen coverage  as a function of the bulk hydrogen concentra- tion C the framework of the Langmuir-McLean isotherm is defined as [12,35]: where C is the impurity mole fraction in units and is the Gibbs free energy.As defined in Equations ( 21) and ( 22), the hydrogen embrittlement parameter is dependent on the selected value of Gibbs free energy and hydrogen damage factor.The curves of hydrogen coverage with hydrogen concentration for various levels of Gibbs free energy range of 10 to 60 k/mol are shown in Figure 1a.Clearly, hydrogen concentration covers about four orders of magnitude for two parameters.The lower boundary and the upper boundary represent the hydrogen concentration at fracture initiation and the ultimate saturation level, respectively.Therefore, the value of the Gibbs free energy is very important and cannot be determined arbitrarily.The Gibbs free energy is assigned a value of 30 kJ/mol, representing the initiation hydrogen concentration and an ultimate saturation level ranging between 0.01 and 100 ppm in the α-Fe grain boundary [13].The Gibbs free energy ∆G 0 g is fixed at a value of 30 kJ/mol.The curves of hydrogen embrittlement G C (θ)/G C (0) with hydrogen concentration for various levels of the hydrogen damage factor χ range of 0.30 to 1.00 are shown in Figure 1b.Considering the initiation hydrogen concentration ranging between 0.01 and 100 ppm, the hydrogen embrittlement parameter G C (θ)/G C (0) values decrease with the value of χ increasing.That is to say, the larger the χ value, the lower the critical energy release rate of hydrogen influence, and the hydrogen-induced fracture will easily occur.According to Jiang and Carter's DFT (Periodic Density Functional) calculations [12], for iron, χ value is taken as 0.89.

Finite Element Implementation
To implement the solution in ABAQUS, the displacements, phase-field, and hydrogen concentration nodal values are interpolated using Voigt notation: where m denotes the number of nodes, and N i is the shape function associated with node i, u i = u x , u y T is the displacement vector at the node i.The shape function matrix is expressed as: The gradient quantities associated with displacements, phase field, and hydrogen concentration are as follows: where B u i stands for the strain-displacement matrices and B i represents the spatial derivatives of the shape functions, which can be described as follows: where N i,x and N i,y are the derivatives of the corresponding shape function with respect to x and y, respectively.

Finite Element Discretization of the Deformation Phase-Field
The residual vectors for the displacement and phase-fields can be considered that Equation (9) must hold for arbitrary values of δu and δφ, which can be expressed as [29]: The value of the parameter k is set as the minimum to help in converging the solution.Normally, a value of k = 1 × 10 −7 is used [31].
By differentiating the residuals for the incremental nodal variables, the components of the consistent stiffness matrices for the displacement field and phase field can be developed as: where C ep in Equation ( 28) is the linear elastic stiffness matrix.

Finite Element Discretization of the Hydrogen Diffusion
Equation (20) can be discretized to generate the corresponding residual vector for any arbitrary virtual variation of the hydrogen concentration δC: From which can be used to create a diffusivity matrix: To solve the diffusion matrix, it is necessary to obtain the hydrostatic stress gradient ∇σ H , which is computed from integration points σ H values, extrapolated to the nodes using each element shape function, and subsequently multiplied by B i to compute ∇σ H .
According to the standard finite element formulation using the Galerkin method, we can obtain the discretized hydrogen transport equation with an unknown vector for the hydrogen concentration C, as follows [36]: The concentration capacity matrix and the diffusion flux vector are given by: where N i is the shape function associated with node i, q the hydrogen flux vector on the surface with an outward normal unit vector J • n.
We can perform the transient hydrogen diffusion analysis by solving Equation ( 33) by means of discretizing the time-derivative of C using the backward Euler method.

Finite Element Implementation in ABAQUS
To solve the coupled mechanical diffusion phase field problem (Equations ( 30) and ( 33)), the finite element (FE) approach is implemented in ABAQUS via user subroutine UEL.The element's characteristics include four integration points and isoparametric 2D quadrilateral quadratic with four degrees of freedom per node, namely u x , u y , φ and C. In ABAQUS/VIEWER, the user element subroutine has the drawback that integration point variables cannot be visualized, considering this limitation, we adopted an auxiliary dummy mesh taking the same number of nodes at integration points that are used in the user-defined element, (i.e., CPE8R).A user material subroutine (UMAT) allows the user to create the constitutive matrix and estimate stresses from the strain data.In this auxiliary mesh, the stress components and the constitutive matrix are made equal to zero, (i.e., they have no influence on the solution of the global system).The corresponding material response is calculated at each integration point in the auxiliary mesh.The data from the developed user subroutine (UEL) is saved in a Fortran module and exported in the UMAT subroutine.The equivalence between model variables and SDVs (Solution-Dependent state variables) are shown in Table 1.

Numerical Modeling
In Section, we have considered the compact tension CT specimens of 45CrNiMoVA steel that were hydrogen-charged to track the experimental load-displacement curves and characterize failure microstructural characteristics.Afterward, a hydrogen-assisted CT specimen failure model is developed to compare the simulation result with experimental data and to demonstrate the model's predictive capability for hydrogen-assisted brittle fracture.

Material and Experiment
A round bar of 45CrNiMoVA high-strength steel (∅ = 60 mm) is used for experimentation.Before fabricating specimens, the material is heated (860 • C, 1 h) and quenched in oil, then tempered (460 • C, 1 h) and gradually cooled with air.The stress-strain curves were performed on standard cylindrical specimens (∅ = 6 mm) by laboratory tensile testing.The Young's modulus and Poisson's ratio of this steel are used as 202 GPa and 0.3, respectively.Tensile testing is performed on a compact tensile CT specimen with an initial notch a = 15 mm, the width W = 30 mm, and the thickness B = 7.5 mm, as illustrated in Figure 2a.CT specimens are designed with their orientations parallel to the radial direction of the round bar.All specimens were undergone under fatigue loading at a frequency of 5 Hz and a load ratio R = 0.1 until the desired crack length versus width ratio of a/W = 0.6 was obtained, as specified in (American Society of Testing Materials) ASTM E1820 [37].The electrochemical approach was used to pre-charge all of the specimens in 0.1 mol• l −1 NaOH solution with a constant current density of 2 mA•cm −2 .The hydrogen atom diffusion distance ( x ) must be viewed as a function of the ideal diffusivity ( D ) and hydro- gen charging time ( t ).This can be expressed as [38]: The electrochemical approach was used to pre-charge all of the specimens in 0.1 mol• l −1 NaOH solution with a constant current density of 2 mA•cm −2 .The hydrogen atom diffu-sion distance (x) must be viewed as a function of the ideal diffusivity (D) and hydrogen charging time (t).This can be expressed as [38]: Low alloy martensitic steels have a hydrogen diffusion coefficient (D) around 10 −5 mm 2 •s −1 [39].Equation (34) uses x = 3.75 mm and D = 10 −5 mm 2 •s −1 with hydrogen charging time 24 h.In this research, the hydrogen pre-charging times were used as 12, 24, and 36 h.Prior to charging, the pre-charged duration was chosen to achieve a variable starting hydrogen concentration in order to study the mechanical properties of materials.In order to evaluate the effect of the extended diffusion period, hydrogencharged specimens were tested at displacement rates of 0.002 mm•s −1 .Additional testing was performed on an uncharged specimen to serve as a baseline for comparison.The MTS 632.02F-20COD extensometer (MTS Systems Corporation, Eden Prairie, MN, USA) was used to gauge the crack growth with the compliance method.Measurement of the crack length, a, was obtained by compliance method based on the following equations [40]: where a is the crack length, W is the specimen width, E is Young's modulus, V 0 is the crack mouth opening displacement, and P is the applied load.The calculation of load-line displacement (V ll ) was based on the linear elastic relationship with the crack mouth opening displacement, which is represented as follows [41]: Figure 2b depicts the load-line displacement curves corresponding to the specimens.The black dotted vertical lines show the maximum load and the accompanying displacement of the curves' load-lines.The 45CrNiMoVA alloy steel shows a tensile strength greater than 1500 MPa.The initiation of a crack causes the load to drop rapidly, and an unstable propagation eventually causes a rupture.In a hydrogen charging time of 36 h, the maximum load and the associated load-line displacement are reduced by 37% and 36%, respectively, due to the influence of HE.As a result, when the specimen is charged with hydrogen for 36 h in excess of the theoretical calculation, hydrogen atoms gain enough time to transmit into the core region, and the hydrogen in the specimen might be totally saturated, leading to the degradation of steel mechanical properties.
Using a scanning electron microscope (SEM) with an acceleration voltage of 20 kV, the specimen microstructures were examined.Figure 3a shows the fracture surface of a hydrogen-free specimen.Region 1 is the pre-crack region, while Region 2 is the area where the crack begins to form.The fracture micromechanics are shown in Figure 3b.The totally ductile material contains numerous fine shallow microvoid coalescence (MVC) dimples.
Hydrogen-charged specimen fracture surfaces were investigated using the same method.Specimen charged with hydrogen for 36 h shows typical fracture surfaces, as shown in Figure 4. Figure 4a illustrates a cross-sectional view of the fracture along the whole length and breadth of the specimen.Pre-crack, crack initiation, and crack growth are three zones of the fracture surfaces.As shown in Figure 4b,c, a high-magnification SEM picture of the fracture edge and the central crack initiation region is shown (a).Quasi-cleavages (QC), deep secondary cracks (SCs), and small-scale micro-void coalescence (MVC) are some of the fracture surfaces that make up the fracture mode.A detailed perspective is shown in Figure 4c,d.The fracture surfaces display an intergranular (IG) fracture with clear crystal sugar block morphology and large "chicken claw" traces on the crystal interface, as well as minor tearing edges, and lengthy and wide secondary cracking (SCs) on the grain boundaries.The HEDE mechanism leads to QC or IG fracture when the local stress accentuation and required hydrogen accumulation at the cleavage planes or grain boundaries (GB) exceed a critical value by weakening the cohesive strength of GB for their separations.This has been elaborated in the published literature as IG and QC modes [42,43].The high-strength steel specimens made of 45CrNiMoVA showed that the HEDE mechanism was dominant in the HE.
The black dotted vertical lines show the maximum load and the accompanying displacement of the curves' load-lines.The 45CrNiMoVA alloy steel shows a tensile strength greater than 1500 MPa.The initiation of a crack causes the load to drop rapidly, and an unstable propagation eventually causes a rupture.In a hydrogen charging time of 36 h, the maximum load and the associated load-line displacement are reduced by 37% and 36%, respectively, due to the influence of HE.As a result, when the specimen is charged with hydrogen for 36 h in excess of the theoretical calculation, hydrogen atoms gain enough time to transmit into the core region, and the hydrogen in the specimen might be totally saturated, leading to the degradation of steel mechanical properties.
Using a scanning electron microscope (SEM) with an acceleration voltage of 20 kV, the specimen microstructures were examined.Figure 3a shows the fracture surface of a hydrogen-free specimen.Region 1 is the pre-crack region, while Region 2 is the area where the crack begins to form.The fracture micromechanics are shown in Figure 3b.The totally ductile material contains numerous fine shallow microvoid coalescence (MVC) dimples.Hydrogen-charged specimen fracture surfaces were investigated using the same method.Specimen charged with hydrogen for 36 h shows typical fracture surfaces, as shown in Figure 4. Figure 4a illustrates a cross-sectional view of the fracture along the whole length and breadth of the specimen.Pre-crack, crack initiation, and crack growth are three zones of the fracture surfaces.As shown in Figure 4b,c, a high-magnification SEM picture of the fracture edge and the central crack initiation region is shown (a).Quasicleavages (QC), deep secondary cracks (SCs), and small-scale micro-void coalescence (MVC) are some of the fracture surfaces that make up the fracture mode.A detailed perspective is shown in Figure 4c,d.The fracture surfaces display an intergranular (IG) fracture with clear crystal sugar block morphology and large "chicken claw" traces on the crystal interface, as well as minor tearing edges, and lengthy and wide secondary cracking (SCs) on the grain boundaries.The HEDE mechanism leads to QC or IG fracture when the local stress accentuation and required hydrogen accumulation at the cleavage planes or grain boundaries (GB) exceed a critical value by weakening the cohesive strength of GB for their separations.This has been elaborated in the published literature as IG and QC modes [42,43].The high-strength steel specimens made of 45CrNiMoVA showed that the HEDE mechanism was dominant in the HE.

The Finite Element Model
The two-dimensional (2D) numerical model of the compact tension (CT) specimen is developed in ABAQUS.Due to symmetry, only the upper half was modeled with 5788 quadratic quadrilateral plane strain elements (CPE8R), as illustrated in Figure 5.The mesh

The Finite Element Model
The two-dimensional (2D) numerical model of the compact tension (CT) specimen is developed in ABAQUS.Due to symmetry, only the upper half was modeled with 5788 quadratic quadrilateral plane strain elements (CPE8R), as illustrated in Figure 5.The mesh size surrounding the fracture tip and region of crack extension is refined in the expected crack propagation domain, where the characteristic element size is 0.0625 mm.The length scale parameter l is normally taken two times the characteristic element length around the fracture route [30].We adopt a length scale 10 times larger than the characteristic element length to ensure enough elements within the process zone, l = 0.625 mm.The experimental tests are performed to verify the material parameters applied in the numerical model: Young's modulus 202 GPa and Poisson's ratio 0.3.The hydrogen damage coefficient x = 0.89 defined by Jiang and Carter [12] is applied.The diffusion coefficient D = 4 × 10 −5 mm 2 •s −1 is examined, as reported by Wu et al. [44].A velocity load F was applied at the loading hole and a data collection point was placed on the model loading line V ll .We assumed the initial hydrogen concentration as C(t = 0) = C 0 evenly distributed throughout the specimen.All exterior borders, including the crack surface, are prescribed a constant concentration, as specified by C(t = 0) = C 0 .The displacement rate was set to 10 −9 mm•s −1 for t = 1 × 10 7 s.This time step was selected for hydrogen redistribution inside the fracture process zone.With a magnitude G c (0) = 78 MPa•mm, the critical energy release rate without hydrogen offered the best fit to the data.

Result and Discussion
By using the presented model, the influential analyses of the hydrogen-assisted brittle fracture through varying phase-field parameters were conducted, and the results were compared with the experimental curves.As a result, for hydrogen-charged times of 12 h, 24 h, and 36 h, the initial hydrogen content values are 2.0 ppm, 7.0 ppm, and 11 ppm, respectively.Based on simulation and experimental data, the load-line displacements (Vll) are depicted in Figure 6.Maximum response force and major curve characteristics are nearly identical in this figure and corroborate with experimental results, i.e., specimen load capacity declines with increasing hydrogen concentration.The HEDE mechanism adopted in this study to diffuse hydrogen into the phase-field model can be used numerically to calculate brittle fracture.

Result and Discussion
By using the presented model, the influential analyses of the hydrogen-assisted brittle fracture through varying phase-field parameters were conducted, and the results were compared with the experimental curves.As a result, for hydrogen-charged times of 12 h, 24 h, and 36 h, the initial hydrogen content values are 2.0 ppm, 7.0 ppm, and 11 ppm, respectively.Based on simulation and experimental data, the load-line displacements (V ll ) are depicted in Figure 6.Maximum response force and major curve characteristics are nearly identical in this figure and corroborate with experimental results, i.e., specimen load capacity declines with increasing hydrogen concentration.The HEDE mechanism adopted in this study to diffuse hydrogen into the phase-field model can be used numerically to calculate brittle fracture.
respectively.Based on simulation and experimental data, the load-line displacem are depicted in Figure 6.Maximum response force and major curve character nearly identical in this figure and corroborate with experimental results, i.e., s load capacity declines with increasing hydrogen concentration.The HEDE me adopted in this study to diffuse hydrogen into the phase-field model can be used ically to calculate brittle fracture.Fracture phase-field values and hydrogen concentrations near the crack tip before crack initiation are shown in Figure 10.As expected, hydrogen concentration is high on the positive hydrostatic stress, i.e., in the vicinity of the crack tip, while the crack begins to propagate when the phase-field values approach one.The maximum load and displacement are noted as 4810 N and 0.29 mm, respectively, which are less than the 6915 N and 0.39 mm values obtained in the absence of hydrogen.Fracture phase-field values and hydrogen concentrations near the crack tip after a certain degree of crack propagation are shown in Figure 11. Figure 11a shows the new crack surfaces as compared to the results shown in Figure (red contour level).The high hydrostatic stress zone moves with the crack tip, resulting in hydrogen concentration, and the peak moves to the newly formed crack surfaces during the crack propagation process.Figure 11b shows that hydrogen is transported ahead of the crack tip.When comparing crack tip hydrogen concentrations in Figures 10b and 11b, the hydrogen concentration declines with fracture propagation.The hydrogen content drops from 14.30 ppm to 11.88 ppm during the crack propagation.The decrease in hydrogen concentration in the crack tip region is mainly related to the change of hydrostatic stress gradient in the region from crack initiation to crack propagation.In the process, the peak position of hydrostatic stress at the front of the crack tip moves with the crack propagation, and the strain gradient at the front area of the crack tip will be slightly lower than that at crack initiation.Therefore, the corresponding hydrostatic stress peak will also decrease slightly, thus resulting in a slight decrease in the hydrogen concentration at the front edge of the crack tip.Fracture phase-field values and hydrogen concentrations near the crack tip before crack initiation are shown in Figure 10.As expected, hydrogen concentration is high on the positive hydrostatic stress, i.e., in the vicinity of the crack tip, while the crack begins to propagate when the phase-field values approach one.The maximum load and displacement are noted as 4810 N and 0.29 mm, respectively, which are less than the 6915 N and 0.39 mm values obtained in the absence of hydrogen.Fracture phase-field values and hydrogen concentrations near the crack tip after a certain degree of crack propagation are shown in Figure 11. Figure 11a shows the new crack surfaces as compared to the results shown in Figure (red contour level).The high hydrostatic stress zone moves with the crack tip, resulting in hydrogen concentration, and the peak moves to the newly formed crack surfaces during the crack propagation process.Figure 11b shows that hydrogen is transported ahead of the crack tip.When comparing crack tip hydrogen concentrations in Figures 10b and 11b, the hydrogen concentration declines with fracture propagation.The hydrogen content drops from 14.30 ppm to 11.88 ppm during the crack propagation.The decrease in hydrogen concentration in the crack tip region is mainly related to the change of hydrostatic stress gradient in the region from crack initiation to crack propagation.In the process, the peak position of hydrostatic stress at the front of the crack tip moves with the crack propagation, and the strain gradient at the front area of the crack tip will be slightly lower than that at crack initiation.Therefore, the corresponding hydrostatic stress peak will also decrease slightly, thus resulting in a slight decrease in the hydrogen concentration at the front edge of the crack tip.
decrease in hydrogen concentration in the crack tip region is mainly related to the change of hydrostatic stress gradient in the region from crack initiation to crack propagation.In the process, the peak position of hydrostatic stress at the front of the crack tip moves with the crack propagation, and the strain gradient at the front area of the crack tip will be slightly lower than that at crack initiation.Therefore, the corresponding hydrostatic stress peak will also decrease slightly, thus resulting in a slight decrease in the hydrogen concentration at the front edge of the crack tip.

Conclusions
This paper presents a hydrogen-assisted fracture methodology, which incorporates the HEDE mechanism into both hydrogen transport and a phase-field-based formulation The distribution of the phase-field value after some amount of fracture propagation is shown in Figure 12, and the associated displacements are 0.38, 0.48, and 0.80 mm, respectively.At this point, the phase-field value reaches one, signaling the start of the crack.The red color indicates newly formed crack surfaces.The crack begins at the tip of the notch and rapidly propagates along the straight route.Meanwhile, the specimen undergoes a rapid load reduction.

Conclusions
This paper presents a hydrogen-assisted fracture methodology, which incorporates the HEDE mechanism into both hydrogen transport and a phase-field-based formulation

Conclusions
This paper presents a hydrogen-assisted fracture methodology, which incorporates the HEDE mechanism into both hydrogen transport and a phase-field-based formulation to study hydrogen embrittlement.The numerical model is demonstrated through hydrogen-charged (CT) 45CrNiMoVA high-strength steel specimens and corroborates with the experimental results.The following are the main conclusions of this research: 1.
The simulation results show that hydrogen accumulates near the crack tip due to positive hydrostatic stress and the peak increases gradually with loading before crack initiation.

2.
The HEDE was implemented by determining the critical energy release rate drops when hydrogen concentration increases.In the presented simulations, the hydrogen concentration reaches a peak near the newly formed crack surfaces and gradually falls as the crack propagates.For load-line displacement curves, the maximum loadcarrying ability decreases as the hydrogen content increases.

3.
The microstructural fracture mechanism of the hydrogen-charged 45CrNiMoVA highstrength steel compact-tension (CT) specimens demonstrate a brittle mixed fracture mode of QC and IG fracture that is consistent with the HEDE mechanism in the suggested model.

4.
The simulated load-line displacement curves show good agreement with computational and experimental curves.The model quantitatively estimates the initial hydrogen level.The proposed model provides a numerical tool for predicting the mechanical reaction of materials that are subjected to hydrogen-assisted brittle fracture, provided that the mechanical properties and phase-field model parameters are properly calibrated in advance.
The hydrogen-enhanced localized plasticity (HELP) is not accounted for in this model.To accurately forecast the hydrogen-assisted fracture ductile to brittle transition, this will be considered in future research.
Fracture Theory Based on Phase-Field Model 2.1.Phase Field Approximation of Diffusive Crack Topology Miehe et al. suggested the phase-field model of the diffusive crack [21-23].The model was developed by considering an infinite one-directional bar having a constant cross-sectional area Γ along the x-axis and occupying a solid body Ω.Where Γ denotes the completely fractured crack surface.The sharp crack topology can be characterized by a field variable φ(x) ∈ [0, 1].

Figure 1 .
Figure 1.(a) The hydrogen coverage  as a function of hydrogen concentration for a range of Gibbs free energy values of 10 kJ/mol to 60 kJ/mol; (b) hydrogen embrittlement parameter

Figure 1 .
Figure 1.(a) The hydrogen coverage θ as a function of hydrogen concentration for a range of Gibbs free energy values of 10 kJ/mol to 60 kJ/mol; (b) hydrogen embrittlement parameter G C (θ)/G C (0) as a function of hydrogen concentration for a range of hydrogen damage factor χ values of 0.30 to 1.00 (∆G 0 g = 30 kJ/mol ).

Figure 3 .
Figure 3. (a) Scanning electron microscopy (SEM) images for 45CrNiMoVA low alloy steel CT specimen, and (b) SEM high magnification image of the selected region in (a).

Figure 3 .
Figure 3. (a) Scanning electron microscopy (SEM) images for 45CrNiMoVA low alloy steel CT specimen, and (b) SEM high magnification image of the selected region in (a).

Figure 4 .
Figure 4. SEM images of the fracture surface of 45CrNiMoVA low alloy steel CT specimen Hcharged for 36 h: (a) overview of the fracture surface; (b,c) high magnification SEM images of the indicated region in image (a); (d) high magnification SEM image of the indicated region in image (c); (red arrowheads indicate secondary cracks (SCs) and white arrowheads indicate tear ridges).QC: quasi-cleavage; IG: intergranular fracture; MVC: microvoid coalescence.

Figure 4 .
Figure 4. SEM images of the fracture surface of 45CrNiMoVA low alloy steel CT specimen Hcharged for 36 h: (a) overview of the fracture surface; (b,c) high magnification SEM images of the indicated region in image (a); (d) high magnification SEM image of the indicated region in image (c); (red arrowheads indicate secondary cracks (SCs) and white arrowheads indicate tear ridges).QC: quasi-cleavage; IG: intergranular fracture; MVC: microvoid coalescence.

Figure 5 .
Figure 5. Finite element model for the (CT) specimen.

Figure 6 .
Figure 6.Load-line displacement curves for uncharged and hydrogen-charged CT specim

Figure 6 .
Figure 6.Load-line displacement curves for uncharged and hydrogen-charged CT specimens.

Figures 7 -
Figures 7-9 illustrates the distributions of hydrostatic stress and hydrogen concentration along the fracture path at various load steps to determine hydrogen concentration at 7.0 ppm.The two distributions are found consistent, exhibiting a peak value near the crack tip, and the peak gradually increases with loading before crack initiation.With the displacement loading from 0.10 to 0.214 mm, the hydrogen concentration at the crack tip increases from 9.189 to 10.76 ppm and the corresponding hydrostatic stress are found as 2568 to 6034 MPa, respectively.The modeling results indicate that hydrostatic stress is the main cause of hydrogen diffusion at the crack tip.Hydrogen atoms concentrate in places with a high local hydrostatic stress gradient, as projected by the hydrogen diffusion flux formula Equation (19).

Figure 9 .
Figure 9. (a) Hydrostatic stress and (b) hydrogen concentration contours near the crack tip at u = 0.214 mm for C 0 = 7.0 ppm.

Figure 10 . 18 Figure 10 .CFigure 11 .Figure 12 .
Figure 10.(a) Distributions of crack phase-field value and (b) hydrogen concentration at crack tip before the onset of crack initiation for C 0 = 7.0 ppm.

Figure 11 .
Figure 11.(a) Distributions of crack phase-field value and (b) hydrogen concentration at crack tip after some amount of crack propagation for C 0 = 7.0 ppm.

Metals 2022 , 18 Figure 10 .C = 7 Figure 11 .Figure 12 .
Figure 10.(a) Distributions of crack phase-field value and (b) hydrogen concentration at crack tip before the onset of crack initiation for 0 C = 7.0 ppm

Table 1 .
List of solution-dependent state variables in two-dimensional.