Phase-field modelling of hydrogen embrittlement in high-strength ferritic steels under cyclic loading

Abstract

Hydrogen-assisted fatigue cracking limits the use of high-strength ferritic steels in pressure vessels, pipelines, and structural components exposed to cathodic protection or hydrogen-bearing environments. We present a coupled phase-field and crystal-plasticity framework for predicting hydrogen-assisted crack advance in S690QL ferritic steel under cyclic loading. The model combines stress-assisted hydrogen diffusion, reversible trapping at dislocations and grain-boundary sites, a fatigue-degraded fracture energy, and an orientation-dependent slip resistance calibrated from electron backscatter diffraction and compact-tension fatigue tests. Electrochemical pre-charging and Devanathan-Stachurski permeation experiments were used to constrain the initial hydrogen inventory and apparent diffusivity, while thermal desorption analysis provided an independent check on trap occupancy. Across three nominal hydrogen contents (0.5, 1.0, and 2.0 wppm), two stress ratios, and loading frequencies from 0.01 to 1 Hz, the model predicted crack-growth rates with a mean absolute percentage error of 7.4% and reproduced the observed transition from transgranular cleavage-like facets at high frequency to intergranular crack advance at low frequency. Sensitivity analysis indicates that the apparent diffusivity, grain-boundary trap binding energy, and fatigue degradation length scale dominate uncertainty in predicted da/dN. The results show that phase-field fracture can provide a physically interpretable tool for hydrogen-service assessment when transport, trapping, microstructure, and cyclic damage are calibrated together rather than fitted independently.

Introduction

Hydrogen embrittlement remains a central reliability concern for steels deployed in hydrogen production, transmission, storage, and mechanically demanding transport applications [1,2]. High-strength ferritic steels are attractive because they combine moderate alloy cost with high yield strength, weldability, and established manufacturing routes. The same microstructural features that raise strength, including high dislocation density, fine lath packets, and segregation-prone boundaries, also provide trapping sites that can localise hydrogen near crack tips and interfaces [3,4]. Under monotonic loading, this localisation may reduce cohesive strength, enhance localised plasticity, or promote hydride-free decohesion depending on alloy, hydrogen fugacity, and loading history [5,6]. Under cyclic loading, the problem is more difficult: crack-tip plasticity, trap occupancy, and boundary damage evolve over many thousands of increments, while the hydrogen field may be out of equilibrium with the mechanical field at each point in the cycle.

Engineering assessment methods still rely heavily on empirical knock-down factors, conservative threshold stress-intensity estimates, or tests conducted under one prescribed charging condition. These approaches are useful for qualification, but they do not reveal why the same steel may exhibit transgranular crack advance at one frequency and intergranular advance at another. Mechanistic models must represent both hydrogen transport and damage evolution across atomistic, microstructural, and continuum scales [13]. Oriani-type trap equilibrium models are widely used for steels [7], and continuum descriptions of stress-assisted diffusion near cracks have a long history [8]. However, coupling these descriptions to a crack-growth model without imposing the crack path a priori remains challenging.

Phase-field fracture provides a natural route to this coupling because the crack is represented by a continuous damage variable rather than by an explicitly tracked discontinuity [9]. Recent extensions have incorporated hydrogen-dependent fracture resistance and have shown promising agreement with hydrogen-assisted cracking experiments [10,11]. The present study extends this class of models to cyclic loading of a high-strength ferritic steel. Our aim is not to introduce a new fitting law for fatigue crack growth, but to test whether independently constrained transport and microstructural parameters can reproduce the measured frequency dependence of da/dN and the observed fracture-mode transition.

Materials and experimental procedure

The investigated material was commercial S690QL plate supplied in the quenched-and-tempered condition. Optical emission spectroscopy confirmed that the composition was within the EN 10025-6 tolerance band; the measured carbon and manganese contents were 0.17 wt% and 1.42 wt%, respectively. Tensile coupons machined transverse to the rolling direction gave a 0.2% proof stress of 887 +/- 11 MPa and ultimate tensile strength of 943 +/- 9 MPa. Prior-austenite grain boundaries were reconstructed from EBSD maps collected on a JEOL JSM-7900F field-emission SEM with an Oxford Instruments Symmetry detector using a 0.18 um step size. The median reconstructed prior-austenite grain diameter was 11.6 um, with martensitic packet boundaries forming the dominant high-angle boundary population.

Compact-tension C(T) specimens with W = 50 mm and B = 10 mm were machined according to ASTM E647-23. Specimens were fatigue pre-cracked in laboratory air to a/W = 0.48 +/- 0.01 before hydrogen charging. Hydrogen was introduced cathodically in 0.1 M NaOH containing 1 g/L thiourea at current densities selected to produce nominal bulk hydrogen contents of 0.5, 1.0, and 2.0 wppm. These contents were measured by thermal desorption analysis immediately after charging and again after mechanical testing. Devanathan-Stachurski permeation tests on 1 mm membranes from the same plate gave an apparent diffusivity of (7.2 +/- 0.3) x 10^-10 m^2/s at 23 deg C, consistent with literature values for high-strength ferritic steels [12].

Fatigue tests were conducted on an Instron 8801 servo-hydraulic frame at stress ratios R = 0.1 and R = 0.5. Loading frequencies were 1, 0.1, and 0.01 Hz. Crack length was monitored by direct-current potential drop and verified by optical measurements after unloading. Fracture surfaces were imaged by SEM at 5 kV, and selected regions were cross-sectioned by focused ion beam lift-out for transmission electron microscopy. Because hydrogen cannot be quantified by conventional EDS, local hydrogen redistribution was assessed indirectly by thermal desorption spectra, crack-growth kinetics, trap-informed diffusion modelling, and, for two representative samples, deuterium charging followed by time-of-flight secondary ion mass spectrometry.

Phase-field formulation

The mechanical problem was formulated in finite strain crystal plasticity with a scalar phase-field damage variable d bounded between zero and one. The undamaged elastic response used cubic elasticity averaged over local lath-packet orientations, while plastic deformation was resolved on the ferritic {110}<111> and {112}<111> slip families. The crack density functional followed the regularised variational fracture framework of Miehe and co-workers [9]. The degradation function was chosen as g(d) = (1 - d)^2 + k, with residual stiffness k = 10^-7 to maintain numerical conditioning.

Hydrogen transport was described by a lattice concentration c_L and a set of reversible traps representing dislocations, lath boundaries, and reconstructed prior-austenite grain boundaries. Lattice diffusion included the hydrostatic-stress contribution to chemical potential. Trap occupancy followed Oriani equilibrium for the baseline simulations [7], while a kinetic trapping variant was used as a robustness check at the lowest loading frequency. The local critical energy release rate was degraded as a function of hydrogen coverage, G_c(theta) = G_c0(1 - chi theta), where theta is the normalised occupancy of the dominant boundary trap population. The parameter chi was calibrated only once, using the 1.0 wppm, R = 0.1, 0.1 Hz dataset, and then held fixed for all remaining predictions.

Cyclic degradation was introduced through a history variable that accumulates positive elastic energy density over repeated loading. The fatigue term reduces the effective fracture toughness once the accumulated cyclic energy exceeds a threshold calibrated from hydrogen-free fatigue crack-growth data. This separation is important: the Paris-law-like fatigue response is calibrated without hydrogen, and the hydrogen terms then alter crack-tip resistance through transport and trapping rather than through a direct empirical multiplication of da/dN.

Numerical implementation and calibration

The coupled problem was implemented as an ABAQUS/Standard 2023 user element with monolithic solution of displacement, phase-field damage, lattice hydrogen concentration, and accumulated fatigue history. A staggered scheme was also tested, but the monolithic solver was retained because it reduced cycle-to-cycle drift in the hydrogen field near the crack tip. The smallest element size in the crack-process zone was 0.8 um, approximately one fifth of the selected phase-field length scale. Mesh refinement studies using 0.6, 0.8, 1.2, and 1.6 um minimum element sizes changed predicted da/dN by less than 4% once the minimum element size was below 1.0 um.

Boundary conditions reproduced the measured load histories from the test frame rather than idealised sinusoidal amplitudes. The initial hydrogen concentration field was taken from the charging simulations and allowed to redistribute during mechanical cycling. Boundary hydrogen exchange during testing was neglected for the baseline case because the free-surface flux estimated from interrupted TDA was smaller than the crack-tip redistribution flux over the simulated cycle window. This assumption was relaxed in a supplementary calculation, which changed predicted da/dN by less than 3% for the 1 Hz cases and by 6% for the 0.01 Hz cases.

Model calibration used a hierarchical procedure. First, elastic-plastic parameters were fitted to monotonic tension and cyclic stabilisation loops in air. Second, the apparent diffusivity and trap densities were constrained by permeation and desorption measurements. Third, hydrogen-free fatigue data fixed the fatigue history threshold and cyclic degradation exponent. Finally, the hydrogen degradation coefficient chi and grain-boundary trap binding energy were adjusted within literature-supported bounds to match one intermediate-frequency dataset. No parameters were re-fitted for the remaining hydrogen contents, stress ratios, or loading frequencies.

Results

The coupled model reproduced the measured crack-growth trends across all tested hydrogen contents. For the baseline calibration case, predicted da/dN agreed with experiment to within 5.9%. Across the full matrix, the mean absolute percentage error was 7.4%, with R^2 = 0.94 for log(da/dN). The largest error occurred in the 2.0 wppm, 0.01 Hz condition, where the model under-predicted crack growth by 12.8%. This discrepancy coincided with the highest fraction of intergranular facets and may indicate local irreversible trap saturation or oxide-assisted boundary weakening not represented in the present formulation.

Frequency had a strong effect on fracture morphology. At 1 Hz, fracture surfaces were dominated by transgranular cleavage-like facets and shallow fatigue striations. The predicted crack path followed high resolved-shear regions within packets and did not dwell at reconstructed prior-austenite boundaries. At 0.01 Hz, SEM fractography showed a continuous intergranular network over 43 +/- 6% of the crack extension region in the 2.0 wppm condition. The simulation reproduced this transition by predicting elevated boundary trap occupancy and reduced local fracture resistance at prior-austenite boundaries intersecting the cyclic plastic zone.

The corrected transport timescale helps explain why the transition is not controlled by lattice diffusion over a single 5 um packet dimension alone. Using the measured apparent diffusivity, L^2/D for L = 5 um is approximately 0.035 s. Hydrogen can therefore redistribute over packet-scale distances within each cycle even at 1 Hz. The rate-limiting process in the simulations is instead the coupled evolution of trap occupancy, cyclic plastic strain, and boundary damage over the larger crack-tip process zone. For an effective redistribution length of 50-80 um, L^2/D is 3.5-8.9 s, which places the 0.1 and 0.01 Hz tests in a regime where boundary occupancy can approach the values required for decohesion during the tensile part of the cycle.

Sobol sensitivity analysis using 10,000 quasi-Monte Carlo samples identified the apparent diffusivity, grain-boundary trap binding energy, and phase-field length scale as the dominant contributors to variance in predicted crack-growth rate [14]. First-order indices were 0.34, 0.27, and 0.16, respectively, with a total interaction contribution of 0.18. The hydrogen degradation coefficient chi was influential only at the highest hydrogen content, whereas the residual stiffness parameter and time-step size had negligible effects within the tested numerical ranges.

Discussion

The results support a picture in which hydrogen-assisted fatigue crack growth is governed by the competition between crack-tip cyclic plasticity, trap-mediated hydrogen redistribution, and boundary-specific fracture resistance. This is consistent with atomistic and continuum studies showing that hydrogen can alter local cohesion and dislocation processes without requiring a single universal embrittlement mechanism [2,6]. The present model does not attempt to distinguish experimentally between HELP and HEDE at the atomic scale. Instead, it uses measurable continuum consequences of trapping and boundary weakening to predict crack advance at the specimen scale.

A useful outcome is that the model captures a fracture-mode transition without prescribing a switch from transgranular to intergranular failure. The phase-field crack follows the weakest available path as the local resistance evolves. At high frequency, the crack advances before extensive boundary trap occupancy accumulates in the process zone. At low frequency, cyclic plasticity and stress-assisted diffusion give boundary traps time to approach high occupancy over a process-zone length, lowering the intergranular resistance. This interpretation also explains why changing R from 0.1 to 0.5 increases da/dN less dramatically than decreasing frequency: the higher mean load raises hydrostatic stress, but frequency controls the time available for coupled trap and damage evolution.

Several limitations remain. The grain-boundary trap population is represented by a reconstructed boundary map and a single effective binding energy, even though real boundaries differ in character, segregation state, carbide decoration, and local chemistry. The model uses isothermal transport and does not include environmental reactions at newly created crack surfaces. The low-frequency under-prediction at the highest hydrogen content suggests that irreversible trapping or local chemical effects may become important near saturation. Finally, the computational cost is high: one 10,000-cycle simulation with the finest mesh required 31 h on 64 CPU cores. Reduced-order surrogates will be needed before this framework can be used routinely in component-scale hydrogen-service assessment.

Conclusion

We developed and tested a coupled phase-field, crystal-plasticity, and hydrogen-transport model for hydrogen-assisted fatigue cracking in S690QL ferritic steel. With transport and fatigue parameters constrained by independent experiments, the model predicted crack-growth rates within 7.4% on average and reproduced the observed shift from transgranular to intergranular crack advance as loading frequency decreased. The analysis indicates that the controlling timescale is not simple lattice diffusion over a single packet dimension, but the coupled evolution of trap occupancy and cyclic boundary damage over the crack-tip process zone.

The framework provides a route toward microstructure-aware assessment of high-strength steels in hydrogen service. Future work will incorporate boundary-character-dependent trap energies, explicit surface reaction kinetics, and validation against welded heat-affected-zone microstructures. The code and calibration datasets are provided as supplementary material to support reproduction of the benchmark calculations.

Data and code availability

All finite-element input decks, post-processing scripts, calibration tables, and raw experimental time histories are included in the supplementary archive. The archive also contains a README describing compiler versions, ABAQUS environment variables, mesh-generation settings, and the exact scripts used to produce the crack-growth and sensitivity-analysis results reported here.