{'fraction_non_alphanumeric': 0.029886695162494065, 'fraction_numerical': 0.012076803039554922, 'mean_word_length': 5.058158651870119, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Applications of phase field fracture in modelling hydrogen assisted failures

The phase field fracture method has emerged as a promising computational tool for modelling a variety of problems including, since recently, hydrogen embrittlement and stress corrosion cracking. In this work, we demonstrate the potential of phase field-based multi-physics models in transforming the engineering assessment and design of structural components in hydrogen-containing environments. First, we present a theoretical and numerical framework coupling deformation, diffusion and fracture, which accounts for inertia effects.Several constitutive choices are considered for the crack density function, including choices with and without an elastic phase in the damage response. The material toughness is defined as a function of the hydrogen content using an atomistically-informed hydrogen degradation law. The model is numerically implemented in 2D and 3D using the finite element method. The resulting computational framework is used to address a number of case studies of particular engineering interest. These are intended to showcase the model capabilities in: (i) capturing complex fracture phenomena, such as dynamic crack branching or void-crack interactions, (ii) simulating standardised tests for critical components, such as bolts, and (iii) enabling simulation-based paradigms such as Virtual Testing or Digital Twins by coupling model predictions with inspection data of large-scale engineering components. The evolution of defects under in-service conditions can be predicted, up to the ultimate failure. By reproducing the precise geometry of the defects, as opposed to re-characterising them as sharp cracks, phase field modelling enables more realistic and effective structural integrity assessments.


Introduction
Hydrogen has been known for decades to notably reduce the toughness, ductility and fatigue life of engineering components [1][2][3]. Hydrogen ingress into a metallic sample can happen during its initial forming, during the coating or plating of a protective layer, through exposure to hydrogen or hydrogen-containing molecules in the air, soil or water, or through corrosion processes. Hydrogen absorbed from gaseous or aqueous environments diffuses within the metal and is attracted to areas of high hydrostatic stress, where damage takes place by means of mechanisms that are still being debated [4][5][6][7].
This so-called hydrogen embrittlement phenomenon is now pervasive across applications in the construction, defence, transport and energy sectors, due to the ubiquity of hydrogen and the higher susceptibility of modern, highstrength alloys [8,9].
The use of fracture mechanics-based models for hydrogen-sensitive applications could be a game-changer in preventing catastrophic failures and optimising material performance. For example, reliable modelling of hydrogen assisted fracture could enable a controlled use of high strength alloys, accelerate material certification, and govern inspection planning and fitnessfor-service assessment. Yet, the development of models capable of predicting crack initiation and growth as a function of material, loading and environmental variables has not been an easy task. Two main challenges hold back the use of predictive models in engineering assessment. The first one is the physical complexity of the problem at hand. Hydrogen embrittlement is a complicated chemical and micro-mechanical phenomenon that involves multiple hydrogen-metal interactions at several scales. However, several mechanistic models have been proposed that show good agreement with experiments with little or no calibration [10][11][12]. Predictions based on nominal material properties and parameters that can be independently determined are now possible. The second challenge lies in developing a computational framework capable of capturing, in arbitrary geometries and dimensions, the multi-physics elements of the problem and their interaction with the complex cracking phenomena occurring in engineering applications. The phase field fracture method [13] appears to provide a suitable framework for overcoming this obstacle.
Martínez-Pañeda et al. [14] have recently extended the phase field fracture method to predict hydrogen assisted failures and the approach has quickly gained popularity as a framework for incorporating various hydrogen embrittlement models [15][16][17][18][19][20]. The main experimental trends have been captured and advanced fracture features such as crack merging, nucleation from arbitrary sites and branching are predicted without convergence problems or the need for remeshing. However, the vast majority of the analyses are restricted to 2D boundary value problems of mostly academic relevance.
In this work, we demonstrate the potential of the phase field fracture method in predicting large scale hydrogen assisted failures of practical engineering interest. This includes, for the first time, (i) the modelling of hydrogen assisted fractures resulting from dynamic loading, where inertia is relevant, and (ii) the consideration of both AT1 [21] and AT2 [13,22] constitutive choices for the dissipation function of the phase field fracture method. Overall, the goal is to showcase the capabilities of phase field-based hydrogen assisted cracking formulations in enabling Virtual Testing in hydrogen-sensitive applications.
To this end, case studies will be addressed involving (i) crack branching in a hydrogen embrittled plate subjected to dynamic loading, (ii) crack-void interactions in a hydrogen-containing 3D bar, (iii) failure of a screw anchor exposed to an aggressive solution, simulating a standardised experiment, and (iv) cracking evolution in a pipeline with internal defects, as measured by in-line inspection. The last example showcases the possibility of combining phase field modelling with inspection data to create (so-called) Digital Twins of critical infrastructure, minimising expensive testing and monitoring. In this regard, there is a further motivation for the use of phase field for engineering assessment. Unlike other computational approaches, such as discrete fracture methods [23][24][25], predictions are not restricted to the evolution of sharp cracks but the growth of defects of any arbitrary shape can be simulated, and without any prior knowledge regarding the extent of growth or the growth direction. This opens the possibility of conducting defect mechanics-based assessments of notch-like defects, significantly reducing the conservatism associated with re-characterising all detected defects into sharp cracks. The concept is illustrated in Fig. 1 for a given defect length (larger than the transition flaw size); simulating the real defect geometry, as opposed to an equivalent crack, can provide more realistic and sustainable criteria for engineering assessment.

Strength
Defect bluntness/radius Cracks Grooves, pits, etc. Local thin areas

Fracture assessment
(Defect re-characterised as a sharp crack)

Plastic collapse assessment
Values used in ¡ tness-for-service assessment Component behaviour Figure 1: Using phase field to enable defect mechanics-based assessments that can reduce overly conservative fitness-for-service assessment. Schematic of current assessment procedures for a defect of length larger than the transition flaw size; adapted from standardised engineering assessment of burst pressure, see Ref. [26].
The remainder of this paper is organised as follows. A generalised phase field framework for chemo-mechanical fracture is given in Section 2, including inertia effects. Details of the finite element implementation are given in Section 3. The four case studies described above are presented in Section 4.
Finally, the manuscript ends with concluding remarks in Section 5.

Potential energy of the solid
Consider a domain Ω ∈ IR n (n ∈ [1, 2, 3]), with outer boundary δΩ ∈ IR n−1 . The domain contains a deformable solid with displacement field u and internal crack surface Γ ∈ IR n−1 . An absorbed species of concentration c might diffuse through the solid and interact with the mechanical behaviour.
The variational energy functional for the solid can be postulated as [27,28]: Here, the first term of the volume integral represents the kinetic energy density of the solid, withu denoting the velocity field and ρ the mass density.
Also within the volume integral, the second and third terms respectively denote the strain energy density and the chemical energy density related to the transport of solute species. The surface integral represents the fracture energy as proposed by Griffith [29], with G c being the material toughness.
The crack geometry Γ is unknown, hindering minimisation of (1). This can be overcome by introducing a continuous phase field variable ϕ ∈ [0, 1]. The phase field resembles a damage variable, representing the intact state of the material when ϕ = 0 and the completely fractured state when ϕ = 1. A regularised potential energy can be defined considering this auxiliary phase field and a regularisation length scale ℓ, such that (1) can be approximated as [13,28]: In the regularized formulation, γ represents a crack surface density function.
In the following, we introduce a number of specific choices for the above formulation, specialising it for brittle fracture of hydrogen-containing metals.
As discussed (e.g.) in Refs. [30,31], the phase field length scale can have different interpretations. As rigorously proven using Γ-convergence, the regularised functional Π ℓ (2) converges to that of Π (1) for a fixed ℓ → 0 + and thus minimising (2) provides the solution for the Griffith variational problem.
Accordingly, ℓ can be interpreted as a regularising parameter in its vanishing limit. However, a finite material strength is introduced for ℓ > 0 + and the phase field length scale becomes a material property governing the material strength, σ c . For example, under plane stress conditions, Thus, the fracture behaviour of the solid is governed by the material toughness G c and strength σ c , as defined by the choice of G c and ℓ. Therefore, phase field models provide a suitable framework to link fracture (in a Griffith sense) and damage. The phase field length scale also regularises the numerical solution; results are mesh-independent if the finite element mesh is sufficiently fine to resolve the fracture process zone, whose size is governed by ℓ.

Linear elasticity
We start by introducing the assumption of linear elastic behaviour. Even within the realm of brittle fracture, this is an appreciable simplification as plasticity will develop locally at the tip of sharp defects. However, the crack tip stresses predicted by conventional plasticity are insufficient to accurately predict crack tip hydrogen concentrations [32,33] and the use of strain gradient plasticity models reveals: (i) the existence of an elastic core surrounding the crack tip [34,35], and (ii) a crack tip stress distribution over the fracture process zone that is closer to that of linear elasticity [36,37]. Thus, linear elasticity provides a conservative, less computationally demanding alternative to multi-scale plasticity models. In addition, the formulation is made under the assumptions of small strains such that the strain tensor is given by, For an (undamaged) linear elastic stiffness tensor, C 0 , the strain energy density of the intact material is defined as Thus, it is assumed that the diffusive species has no influence on the strain energy, which is a common assumption for hydrogen in metals [38].

Phase field fracture
We proceed to make constitutive choices for the phase field fracture formulation. Two models will be considered, which are known to provide accurate descriptions of fracture phenomena in a regularized setting. The strain energy density and the crack surface density function are given by: Where the degradation function g(ϕ) is continuous and monotonic and takes the values g(0) = 1 and g(1) = 0. The function w(ϕ) must fulfill w(0) = 0 and w(1) = 1. Finally, This choice of phase field formulation has been shown to Γ-converge to the Griffith solution [39]. Common choices for g(ϕ) and w(ϕ) include g(ϕ) = (1 − ϕ) 2 and w(ϕ) = ϕ 2 , which produces the so-called standard or AT2 phase field model [22]. Another common choice is the same quadratic degradation function and w(ϕ) = ϕ, which is often referred to as the AT1 phase field model [21]. The latter formulation introduces an elastic regime prior to any damage in the solid.

Modified Fickian diffusion
Finally, for the absorbed diffusive species, we adopt a modified version of Fickian diffusion. In addition to concentration gradients, diffusion is assumed to be driven by gradients of hydrostatic stress, such that atomic hydrogen accumulates in areas where the lattice is being expanded. Mass conservation requirements relate the rate of change of the hydrogen concentration c with the hydrogen flux J through the external surface, Diffusion is driven by the gradient of the chemical potential ∇µ. Thus, for a diffusion coefficient D, the flux is related to ∇µ through a linear Onsager relationship, The chemical potential includes a hydrostatic stress σ H -dependent term to account for the role of volumetric strains in driving diffusion. For an occupancy of lattice sites θ L and partial molar volume of hydrogenV H , the chemical potential of hydrogen in lattice sites is given by, Here, µ 0 denotes the chemical potential in the reference case, R is the gas constant and T is the absolute temperature. Consider now the relation between the occupancy and the number of sites θ L = c/N , and make the common assumptions of low occupancy (θ L ≪ 1) and constant interstitial sites concentration (∇N = 0); inserting (11) into (10) then renders,

Hydrogen degradation of the toughness
The material toughness is defined to be sensitive to the hydrogen content, G c (c), as consistently observed experimentally. First, the Langmuir-McLean isotherm is used to estimate, from the bulk concentration, the hydrogen coverage θ at decohering interfaces, Here, g 0 b is the binding energy for the impurity at the site of interest. A value of 30 kJ/mol is assumed throughout this study, which is representative of hydrogen trapped at grain boundaries [40]. Following atomistic studies of surface energy sensitivity to hydrogen coverage [41], a linear degradation of the fracture energy with θ is assumed, where χ is a hydrogen damage coefficient, a material parameter that can be estimated by calibrating with experiments [18] or inferred from atomistic calculations [14]. We adopt the latter approach and assume a value of χ = 0.89 throughout our calculations, this magnitude provides the best fit to atomistic calculations on iron [14,42]. A failure model based on grain boundary decohesion is implicitly assumed for the above choices: a degradation law based on atomistic calculations of surface energy sensitivity to hydrogen content, and the consideration of a grain boundary binding energy in (13). However, we emphasise that the framework is universal; the degradation law can be adapted to accommodate any other mechanistic interpretation (e.g., including a dependence on parameters such as the dislocation density or the void volume fraction) or chosen to be phenomenological, relating G c to c and to a hydrogen damage coefficient χ to be calibrated.

Coupled force balances
We proceed to present the weak and strong form of the problem considering the constitutive choices above. First, following Refs. [15,43], a scalar field η is defined to determine the kinematics of composition changes, such Thus, from a kinematic viewpoint, the domain Ω can be described by the displacement u, phase field parameter ϕ, and chemical displacement η. Taking the first derivative of (2) and incorporating the constitutive prescriptions adopted, the coupled weak form reads, in the absence of body forces and external tractions and fluxes, as follows: The local force balances can then be readily derived by applying Gauss' divergence theorem and noting that (16) must hold for any kinematically admissible variations of the virtual quantities. Accordingly, The coupling between the different physical elements of the problem is evident in (17). First, as damage increases, the phase field reduces the stiffness of the solid in the linear momentum equation (17a). As observed in (17b), the phase field evolves driven by the competition between strain energy density ψ 0 and toughness G c , the latter diminishing with increasing hydrogen concentration. Finally, (17c), diffusion of atomic hydrogen is governed by concentration gradients and the lattice dilation, as characterised by the hydrostatic stress.

Finite element implementation
We shall now describe the details of the numerical implementation in the context of the finite element method. First, some numerical considerations are presented for the phase field problem in Section 3.1, to guarantee damage irreversibility and prevent crack growth from compressive stresses. Secondly, a threshold is defined in Section 3.2 to address the implementation peculiarities inherent to the AT1 model. Thirdly, the discretisation of the problem and the formulation of residuals and stiffness matrices is described in Sec-

Addressing irreversibility and crack growth in compression
First, a history variable field H is introduced to ensure damage irreversibility. Thus, for a time t, Secondly, we introduce a strain energy decomposition to prevent cracking in compression. The volumetric-deviatoric split by Amor et al. [45] is adopted, by which the compressive volumetric strain energy does not contribute to damage. Thus, the strain energy density is decomposed into the following terms: where K is the bulk modulus, ⟨a⟩ ± = (a ± |a|)/2 and ε ′ = ε − tr(ε)I/3.
In addition, we follow the hybrid implementation of Ambati et. al [46] in considering only ψ + 0 in the evaluation of the history variable field H, therefore referring to it as H + henceforth, while considering ψ 0 in the displacement problem.

Implementing the AT1 phase field formulation
Unlike the AT2 phase field model, the AT1 formulation does not inherently ensure that the lower bound on the phase field is enforced. If no measures are taken, the phase field can become negative for all strains below the critical strain, which is given by To overcome this, we introduce a lower bound by re-defining the history field as: The minimum threshold employed for the history field corresponds to the strain energy density magnitude that yields a vanishing phase field in the homogeneous 1D case [47]. Other methods exist for implementing the AT1 formulation, such as using constrained optimization solvers.

Finite element discretisation of variational principles
We proceed to discretise the linearised problem and present the associated residuals and stiffness matrices. First, making use of Voigt notation, the displacement field u, phase field ϕ and hydrogen concentration c can be discretised as where N i denotes the shape function associated with node i, for a total number of nodes m. Here, N u i is a diagonal interpolation matrix with the nodal shape functions N i as components. Similarly, using the standard straindisplacement B matrices, the associated derivatives are discretised as, Considering this finite element discretisation and the weak form balances (16), the resulting discrete equations of the balances for the displacement, phase field and concentration can be expressed as the following residuals: where k is a numerical parameter introduced to keep the system of equations well-conditioned; a value of k = 1 × 10 −7 is adopted throughout this study based on previous studies. This choice is grounded on previous studies [48,49]; the use of smaller values has no influence in the results. Subsequently, the tangent stiffness matrices are calculated as: Unless otherwise stated, the global system of equations for the linearised problem will be solved using a staggered, alternative minimisation scheme [14,50].

Results
The potential of the formulation in simulating complex fracture phenomena and transforming engineering assessment is demonstrated by addressing four case studies of particular interest. First, we model for the first time dynamic failure of a hydrogen pre-charged steel plate, using the AT1 model; see Section 4.1. Next, in Section 4.2, we address the failure of a tensile bar due to the interaction between a tilted crack and a neighboring void.
Thirdly, in Section 4.3, we model the brittle fracture of an anchor in a concrete element subjected to a corrosive environment, following the ASTM E488 standard [51]. Finally, the progressive failure of a pipeline is simulated; coupling modelling with in-line inspection data, the model incorporates the numerous defects that typically arise due to pitting corrosion and captures their growth and coalescence under in-service conditions.

Crack branching in an embrittled steel plate due to dynamic loading
For the first case study, we consider the paradigmatic boundary value problem of dynamic crack branching in a rectangular plate [28,52,53]. The geometry and boundary conditions are given in Fig. 2. This well-known case study is based on dynamic experiments on brittle materials such as glassy polymers [54]. Here, we aim to illustrate the influence of hydrogen on the dynamic fracture pattern of a martensitic steel of type 440C, which is known to exhibit very brittle fracture in the presence of hydrogen [55].  [56], which corresponds to a Griffith energy of G c = 2.33 J/mm 2 ; recall that K Ic = G c E/(1 − ν 2 ). The material is assumed to be pre-charged with either 5 or 0.1 wppm of hydrogen. As a reference case, a specimen without hydrogen has also been included. The fracture length scale has been chosen as two times the characteristic element length ℓ = 2h e = 0.4 mm and the specimen is subjected to a tensile impact load of 60 MPa. The domain is meshed with quadratic rectangular (square) elements and the initial crack has been introduced through the phase field rather than as a discontinuity in the mesh. A total of 50,000 elements are used, taking advantage of symmetry. In this dynamic case, crack propagation is expected to be faster than hydrogen diffusion by at least an order of magnitude, allowing us to neglect the influence of hydrogen transport. Without the diffusion equations the global stiffness matrix is symmetric, which enables us to make use of a monolithic scheme in conjunction with the quasi-Newton method [57,58] to solve the problem more efficiently. For this specific problem, the AT1 phase field model is used. Thus, this case study also constitutes the first example involving AT1 and quasi-Newton solution methods, confirming the successful performance observed with other phase field models [57,58].
The increment size has been chosen as ∆t < h e /v n to completely resolve the stress waves in the material. The finite element results obtained are reported in Fig. 3 in terms of the phase field contours; blue colour represents intact material (ϕ ≈ 0) while red colour is used to denote the crack (ϕ ≈ 1). The results reveal that the crack pattern is very sensitive to the hydrogen content. A noticeable influ-ence is observed for the case of a c 0 = 0.1 wppm hydrogen concentration and extensive branching is seen for the case of c 0 = 5 wppm, which might be interpreted as mild shattering. We emphasise that branching is intrinsically related to inertia effects, which are significant in this boundary value problem, and the material brittleness, which is sensitive to the hydrogen content.
The results showcase the capability of the proposed framework for capturing inertia effects and the complex crack patterns that might occur in hydrogenassisted fractures. The remaining case studies will deal with quasi-static loading conditions and accordingly the kinetic energy terms will be dropped from the model.

Void-crack interaction in a 3D tensile bar
The second case study deals with the prediction of crack initiation, growth and unstable failure in a rectangular prismatic bar containing a circular void and a tilted circular edge crack. The bar is shown in Fig. 4, with length

Virtual experiments: design of screw anchors against brittle fracture
Our third case study models the failure of screw anchors in aggressive environments, mimicking the testing conditions of the ASTM E488/E488M standard [51]. The aim is to showcase the capabilities of the model in optimising experimental campaigns and certification, while addressing the phe-nomenon of bolt failure due to hydrogen ingress, a significant concern in offshore engineering [59].
Screw anchors for use in concrete are often made out of high-strength galvanized steel. The zinc coating provides excellent corrosion protection but it can potentially increase the risk of hydrogen embrittlement. If the coating is damaged, the corrosion potential can be lowered sufficiently for hydrogen evolution to occur even in the highly alkaline conditions observed in concrete. To simulate the brittle failure of screws exposed to hydrogencontaining environments, we chose to replicate the test provided in the ASTM E488/E488M standard [51] for assessing the susceptibility of screw anchors to hydrogen embrittlement. The setup for the test is illustrated in Fig. 7.
The general approach for the standardised test is to pre-charge the anchor with hydrogen and then carry out a tensile test up to ultimate failure. The pre-charging occurs by exposing the sample, for an extended period of time, to a solution representative of the one found in concrete pores. The precharging is carried out with the anchor in tension and potentiostatic control of the potential. The applied potential is kept sufficiently low for hydrogen evolution to occur. To estimate the hydrogen content of the bolt we consider the work by Recio et al. [60], where a similar pre-charging protocol was employed. Based on the results therein and the differences between protocols, the initial hydrogen concentration is expected to be equal to 2 wppm or higher. A magnitude of 2 wppm is assumed. Both parts of the modeling domain are meshed using tetrahedral elements with quadratic shape functions. The mesh is illustrated in Fig. 8. The concrete domain is discretised with 117,456 while a total of 155,278 elements are used for the screw. As the problem considers crack initiation from a stress concentration, it is deemed suitable to adopt the AT1 phase field model [30].

Coupling with NDE: Failure of a pipeline with internal defects
The final case study showcases how phase field methods can be coupled to non-destructive evaluation (NDE) to develop high-fidelity models for realtime assessment of critical infrastructure. This is exemplified with the modelling of the progressive failure of a pipeline with numerous internal defects.
The structural integrity of pipelines subjected to aggressive environments is a major concern in the energy industry [61]. Corrosion damage nucleates pits that act as stress concentrators, attracting hydrogen that triggers early cracking. Thus, we consider the structural failure of a pipeline as a paradigmatic case study and show how our multi-physics model can predict cracking in large scale components as a function of the material properties, the environmental hydrogen concentration and the loading conditions. In addition, we aim at enabling a Virtual Testing paradigm by coupling the present phase field predictions with defect characterisation from NDE inspections. The distribution of defects in a pipeline is shown in Fig. 11, as measured from in-line inspection (ILI) and reported by Larrosa et al. [26].
A three-dimensional finite element model is developed based on this data. Figure 11: Initial distribution of defects in a corroded pipeline, as measured by in-line inspection (ILI). Adapted from [26].
A total of 1,750 defects were characterised by Larrosa et al. [26] in the pipeline shown in Fig. 11, of length 11 km, outer radius 162 mm, and thickness 40 mm. We mimic the outer radius and thickness of the pipeline but we restrict our attention to a critical span of 2 meters length, containing a total of 112 defects. Taking  The phase field AT2 model is considered. We assume that the pipeline has been pre-charged with an initial hydrogen concentration of c 0 = 1 wppm and is continuously exposed to the same environment, which corresponds to a 3% NaCl aqueous solution. We note that a fixed concentration is prescribed for simplicity but a constant chemical potential should instead be used as boundary condition when aiming at quantitative results [32,62,63]. Mimicking in-service loading conditions in a riser, we subject the pipeline to an internal pressure of 152 MPa and to axial tension, with a remote stress of 105 MPa. The internal pressure is increased linearly in time while the remote stress is held constant throughout the analysis. The evolution of defects and cracks predicted is shown in Fig. 12 by plotting the completely cracked regions (ϕ > 0.8).

Conclusions
We have presented a generalised phase field formulation for predicting hydrogen assisted fracture. The coupled deformation-diffusion-damage framework developed considers, for the first time, both AT1 and AT2 phase field models, the role of inertia, and a 3D finite element implementation.