A simple and robust Abaqus implementation of the phase field fracture method

The phase field fracture method is attracting significant interest. Phase field approaches have enabled predicting - on arbitrary geometries and dimensions - complex fracture phenomena such as crack branching, coalescence, deflection and nucleation. In this work, we present a simple and robust implementation of the phase field fracture method in the commercial finite element package Abaqus. The implementation exploits the analogy between the phase field evolution law and the heat transfer equation, enabling the use of Abaqus' in-built features and circumventing the need for defining user elements. The framework is general, and is shown to accommodate different solution schemes (staggered and monolithic), as well as various constitutive choices for preventing damage under compression. The robustness and applicability of the numerical framework presented is demonstrated by addressing several 2D and 3D boundary value problems of particular interest. Focus is on the solution of paradigmatic case studies that are known to be particularly demanding from a convergence perspective. The results reveal that our phase field fracture implementation can be readily combined with other advanced computational features, such as contact, and deliver robust and precise solutions. The code developed can be downloaded from www.empaneda.com/codes.


Introduction
Modelling the morphology of an evolving interface is considered to be a longstanding mathematical and computational challenge. Tracking interface boundaries explicitly is hindered by the need of defining moving interfacial boundary conditions and manually adjusting the interface topology with arbitrary criteria when merging or division occurs (Biner, 2017). Phase field formulations have proven to offer a pathway for overcoming these challenges.
In the phase field modelling paradigm, the interface is smeared over a diffuse region using an auxiliary field variable φ, which takes a distinct value for each of the two phases (e.g., 0 and 1) and exhibits a smooth change between these values near the interface. The temporal evolution of the phase field variable φ is described by a partial differential equation (PDE) and thus the method enables the simulation of complex interface evolution phenomena by integrating a set of PDEs for the whole system, avoiding the explicit treatment of interface conditions. The phase field paradigm has quickly gained significant traction in the condensed matter and materials science communities, becoming the de facto tool for modelling microstructural evolution (Provatas and Elder, 2011). The change in shape and size of microstructural features such as grains can be predicted by defining the evolution of the phase field in terms of other fields (temperature, concentration, strain, etc.) through a thermodynamic free energy. This success has been extended to other interfacial problems, such as corrosion, where the phase field smoothens the metal-electrolyte interface (Cui et al., 2021), or fracture mechanics, where the phase field is used to implicitly track the evolution of the crack-solid boundary (Bourdin et al., 2000).
The success of phase field fracture methods has also triggered a notable interest for the development of robust solution algorithms to solve the coupled deformation-fracture problem (Miehe et al., 2010b;Gerasimov and De Lorenzis, 2016;Wu et al., 2020a;Kristensen and Martínez-Pañeda, 2020).
The total potential energy functional, including the contributions from the bulk and fracture energies, is minimised with respect to the two primary kinematic variables: the displacement field u and the phase field φ. Thus, the phase field φ, a damage-like variable, is solved for at the finite element nodes, as an additional degree of freedom. This requires performing the numerical implementation at the element level, as opposed to local damage models, which are implemented at the integration point level. In the context of commercial finite element packages, solving for the phase field as a degree-of-freedom requires the development of user element subroutines. The commercial finite element package Abaqus has received particular attention in the phase field fracture community, and a vast literature has emerged on the implementation of the phase field fracture method on this popular software suite (Liu et al., 2016;Molnár and Gravouil, 2017;Fang et al., 2019;Molnár et al., 2020b;Wu and Huang, 2020). These implementations require programming an ad hoc finite element, effectively using Abaqus as a solver and not being able to exploit most of its in-built features. In this work, we circumvent this issue by exploiting the analogy between the heat conduction equation and the phase field evolution law. This approach enables using the vast majority of Abaqus' in-built features, including the coupled temperature-displacement elements from its finite element library, which avoids coding user-defined elements and the associated complications in meshing and visualisation (e.g., Abaqus2Matlab is frequently used to preprocess input files, Papazafeiropoulos et al., 2017). Moreover, the phase field implementation presented can accommodate both staggered and monolithic solution schemes, ensuring convergence in all cases. We demonstrate the potential and robustness of the implementation presented by addressing several paradigmatic 2D and 3D boundary value problems. The framework provided is general and can be easily implemented in other finite element packages.
The remainder of this manuscript is organised as follows. In Section 2 we describe the theory underlying the phase field fracture method. The analogy with the heat transfer problem and the implementation details are given in Section 3. Representative results are shown in Section 4. First, unstable fracture is addressed with the paradigmatic benchmark of a cracked square plate under uniaxial tension. Secondly, convergence under stable crack propagation conditions is investigated using a cracked square plate subjected to shear. The performance of monolithic and staggered schemes is compared.
Thirdly, the screw tension tests presented by Wick et al. (2015) are examined. Finally, we simulate the so-called Brazilian laboratory test, which is widely used for measuring the tensile strength of rock-like materials. A comprehensive 3D analysis is conducted, including the modelling of the contact between the jaws and the specimen. The manuscript ends with concluding remarks in Section 5.

Phase field fracture model
The phase field fracture method builds upon Griffith's thermodynamics framework (Griffith, 1920). In agreement with the first law of thermodynamics, a crack can form (or grow) only if this process causes the total energy of the system to decrease or remain constant. Accordingly, a critical condition for fracture can be defined upon the assumption of equilibrium conditionsno net change in total energy. Consider an elastic solid containing a crack.
In the absence of external forces, the variation of the total energy E due to an incremental increase in the crack area dA is given by where W c is the work required to create new surfaces and ψ is the strain energy density, which is a function of the displacement field u and the strain field ε = ∇u T + ∇u /2. The last term in Eq. (1) is the so-called critical energy release rate G c = dW c /dA, a material property that characterises the fracture resistance. Thus, Griffth's premise is a local minimality principle for the sum of the elastic and fracture energies. For an arbitrary body Ω ⊂ IR n (n ∈ [1, 2, 3]) with internal discontinuity boundary Γ, this minimality principle can be expressed in a variational form as (Bourdin et al., 2008), Thus, within this framework, crack growth along any trajectory can be predicted without arbitrary criteria, driven by global minimality and the transformation of stored energy into fracture energy. However, minimisation of the variational Griffith energy functional (2) is hindered by the complexi-ties associated with tracking the propagating fracture surface Γ. The problem can be made computationally tractable by employing an auxiliary phase field φ that enables tracking the crack interface. The phase field φ can be interpreted as a damage-like variable that goes from 0 in intact regions to 1 inside of the crack. Accordingly, following continuum damage mechanics arguments, a degradation function g(φ) = (1 − φ) 2 can be defined to reduce the material stiffness with evolving damage. Hence, the regularised energy functional is given by, where is a length scale parameter that governs the size of the fracture process zone and γ is the crack density function. A common choice for γ reads, As rigorously proven using Gamma-convergence, the (u, φ) sequence that constitutes a global minimum for the regularised functional E converges to that of E for a fixed → 0 + . Thus, can be interpreted as a regularising parameter in its vanishing limit. However, for > 0 + a finite material strength is introduced and thus becomes a material property governing the strength (Tanné et al., 2018); e.g., for plane stress: where K Ic is the material fracture toughness. It has been shown that the consideration of a finite > 0 + enables to accurately predict crack nucleation, capturing its transition from strength-driven to fracture-driven (Tanné et al., 2018), and in agreement with the predictions from the coupled criterion in finite fracture mechanics (Molnár et al., 2020a).
We will restrict our analysis to the behaviour of linear elastic materials, such that the strain energy density of the intact material is given by, where C 0 is the (undamaged) linear elastic stiffness tensor. Accordingly, the Cauchy stress tensor is defined as where the undamaged Cauchy stress is given by σ 0 = C 0 : ε.
Considering the constitutive choices just described and taking the first variation of the E with respect to the primal kinematic variables u and φ renders, The local force balances can be readily derived by applying Gauss' divergence theorem and noting that (8) must hold for any kinematically admissible variations of the virtual quantities. Thus, the coupled field equations read, The discretised forms of the field equations can be solved using a monolithic scheme, where u and φ are solved simultaneously, or by means of a so-called staggered scheme, where an alternate minimisation strategy is used.

Finite element implementation
We shall describe the numerical framework proposed. First, we introduce a history field to ensure damage irreversibility. Secondly, the analogy with heat transfer is presented. Thirdly, the particularities of the Abaqus implementation are described. Finally, we show how our implementation can accommodate different solution schemes, and discuss the advantages and limitations of the options available. For the sake of brevity, we limit our description to the constitutive and implementation choices inherent to the code provided, and describe in Appendix A other potential extensions, which are considered in the numerical examples.

Damage irreversibility
A history variable field H is introduced to prevent crack healing, ensuring that the following condition is always met where φ t+∆t is the phase field variable in the current time increment while φ t denotes the value of the phase field on the previous increment. For both loading and unloading scenarios, the history field must satisfy the Kuhn-Tucker conditions Accordingly, the history field for a current time t can be written as:

Heat Transfer Analogy
For a solid with thermal conductivity k, specific heat c p and density ρ, the field equation for heat transfer in the presence of a heat source r reads: where T is the temperature field. Under steady-state conditions the rate term vanishes and Eq. (13) is reduced to, The analogy of this elliptic partial differential equation (PDE) with the phase field evolution law is evident, with the temperature field acting as the phase field T ≡ φ. Making use of the history field described above, one can reformulate the phase field local force balance, Eq. (9)b, as And thus (14) and (15) are equivalent upon assigning the value of unity to the thermal conductivity (k = 1) and defining the following heat flux due to internal heat generation, Finally, for the computation of the Jacobian matrix, one should also define the rate of change of heat flux (r) with temperature (T ≡ φ), We have restricted ourselves to the steady-state scenario, treating the phase field evolution law as rate-independent. This is, by far, the most common formulation for phase field fracture. However, one can also introduce a viscous regularisation term in the phase field equation by exploiting instead the transient problem -Eq. (13). In such scenario, the quantity ρc p is analogous to a viscosity parameter (Miehe et al., 2010a). The heat capacity terms help stabilising the solution and thus one might wish to address a rate-independent (steady-state) problem by conducting instead a transient analysis over a long time. However, as demonstrated in the numerical examples below, we do not see the need to consider viscous regularisation to achieve convergence.

Abaqus particularities
The heat transfer analogy described can be readily implemented in Abaqus  To avoid editing the user subroutine, mechanical and fracture properties are defined in the input file only, as user material properties, and are then transferred between subroutines using solution dependent variables. Consistent with the heat transfer analogy outlined above, one must activate the heat generation option and define as material properties the thermal conductivity k, the density ρ and the specific heat c p , all of them with a value of unity. Also, one should assign an initial temperature distribution of No additional pre-processing or post-processing steps are needed, all actions can be conducted within the Abaqus/CAE graphical user interface and the phase field solution can be visualised by plotting the nodal solution temperature (NT11).

Solution schemes
The global system of equations, shown in Fig. 1, can be solved in either a monolithic or a staggered manner. In a monolithic approach, the displacement sub-system K u u = R u and the phase field sub-system K φ φ = R φ are solved simultaneously. On the other hand, a staggered solution scheme entails an alternative minimisation approach, by which the sub-systems are solved sequentially. Monolithic solution strategies are unconditionally stable and, therefore, more efficient (in principle). However, the total potential energy functional (3) is non-convex with respect to u and φ. As a consequence, the Jacobian matrix in Newton's method becomes indefinite, hindering convergence when solving for the displacement and the phase field at the same time. It has been recently shown that the use of quasi-Newton methods such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm enables the implementation of robust monolithic schemes that are very efficient and do not exhibit convergence issues (Kristensen and Martínez-Pañeda, 2020;Wu et al., 2020a) -see also (Kristensen et al., 2020b;Wu et al., 2021) for application examples. Unfortunately, the quasi-Newton solution scheme is not available in Abaqus for thermo-mechanical problems. Accordingly, we implement a conventional monolithic scheme, based on Newton's method, and a staggered scheme of the single-pass type. The flowchart associated with each of these solution schemes is presented in Fig. 2. In the staggered case, the residual and the stiffness matrix for the phase field sub-system are built considering the history field of the previous increment H t ; i.e., the history field is frozen during the iterative procedure, facilitating convergence in demanding problems at the cost of scarifying unconditional stability. A recursive iteration or multi-pass staggered scheme can be implemented by using a Fortran module to transfer the history field between the UMAT and the HETVAL. Thus, we provide a general framework that provides flexibility to enhance robustness or efficiency, as required for the problem at hand. This trade-off between efficiency and robustness, and the differences in performance between solution schemes, are addressed in the numerical examples below.

Results
We shall show the robustness and capabilities of the present implementation by simulating fracture in several paradigmatic boundary value problems.
First, crack initiation and growth in a notched square plate is addressed under both uniaxial tension (Section 4.1) and shear (Section 4.2). Then, the failure of screws subjected to tension, with and without initial cracks, is sim- N/mm (Kristensen and Martínez-Pañeda, 2020). We discretise the model using linear quadrilateral elements for coupled displacement-thermal analyses, CPE4T in Abaqus terminology. A total of 8,532 elements are used. As shown in Fig. 3b, the mesh is refined along the expected crack path, such that the characteristic element size is at least five times smaller than the phase field length scale . For this case study, the monolithic implementation is used and no strain energy decomposition is assumed. The predicted crack path is showcased in Fig. 3c by plotting the contours of the phase field variable φ. The force versus displacement response predicted is shown in Fig. 4. The result agrees with that of Kristensen and Martínez-Pañeda (2020), which was obtained using a quasi-Newton solution scheme. Cracking is unstable, with the crack extending through the ligament instantaneously. This leads to a dramatic drop in the load carrying capacity, as shown in Fig. 4. However, despite this drastic change in the structural response, convergence can be attained and the fracture event is captured in one single load increment. Fig.   4 also shows the number of iterations required to achieve convergence in each increment, superimposed to the force versus displacement response. We use time increments of constant size and resolve the analysis with a total of 100 load increments. Convergence throughout can be achieved with as few as 10 increments, but using a larger number facilitates capturing the sudden load drop with greater fidelity. An adaptive time stepping scheme, such as the one developed by Kristensen and Martínez-Pañeda (2020), can be easily incorporated. This will allow for the increment size to increase or decrease as needed, enabling accurate results at an even smaller computational cost.
In any case, it can be observed that the problem can be solved efficiently, with most time increments requiring a small number of iterations to achieve convergence (10 or fewer). However, resolving the fracture event requires a load increment with over 400 iterations. Unlike other computational fracture methods, the Newton-Raphson algorithm can converge after hundreds of iterations in phase field models (Gerasimov and De Lorenzis, 2016). The solution controls of Abaqus have to be edited to increase the maximum number of iterations that are allowed before convergence is deemed unlikely and the load increment is aborted (see the accompanying input file, to be downloaded from www.empaneda.com/codes). It must be noted that, despite the good performance observed, this boundary value problem can be resolved more efficiently using quasi-Newton solution schemes (see Kristensen and Martínez-Pañeda, 2020).

Notched square plate under shear
We shall now address the case of stable crack growth by simulating the fracture of the notched square plate considered in Section 4.1, but subjected to shear loading. As shown in Fig. 5a The force versus displacement response is shown in Fig. 6, along with the size of each increment and the number of iterations that were needed to achieve convergence. The crack propagates in a stable manner, leading to a progressive reduction in the reaction force. Again, the results agree with those obtained by Kristensen and Martínez-Pañeda (2020) using a monolithic quasi-Newton solution scheme. This boundary value problem is known to be particularly challenging from a convergence viewpoint and is thus used to compare the monolithic and staggered solution schemes. Consider first the monolithic analysis, Fig. 6a. While the entire crack propagation process can be captured, many increments require a very significant number of iterations to achieve convergence -unlike in the uniaxial tension case where cracking is unstable. It is clear that, for this boundary value problem, the monolithic implementation struggles to converge and becomes inefficient. Now let us examine the output of the staggered case. The results obtained with the single-pass staggered implementation also make use of a uniform increment size, with the entire analysis being conducted using 10 4 load steps. This is a sufficiently large number of increments such that the solution is similar to that obtained with the unconditionally stable monolithic model -see Fig. 6b.
In the staggered case, all load steps converge after two increments. Notwithstanding, as discussed before, this solution scheme is not unconditionally stable and results can be sensitive to the number of time increments. We also conduct the analysis using 10 3 load steps; the crack trajectory and the maximum force attained agree with those predicted with the monolithic scheme but the force versus displacement result differs in the softening region (not shown). The staggered implementation appears to be more robust and efficient than the monolithic one for this specific case study; as quantified in Fig.   6c, the total number of iterations is larger in the monolithic case. However, one should note that both implementations are significantly outperformed by a monolithic approach based on the quasi-Newton solution method. As shown in (Kristensen and Martínez-Pañeda, 2020), a precise solution to this specific boundary value problem can be obtained with a number of iterations that is one order of magnitude smaller than the accurate staggered solution.

Screw tension tests
We proceed now to simulate the fracture of a screw subjected to tension, following the work by Wick et al. (2015). The geometry, dimensions and boundary conditions mimic those of (Wick et al., 2015) and are shown in Fig. 7. Three different cases are considered. First, we model a screw with no initial damage; i.e., without the initial crack displayed in Fig. 7. Secondly, we will assume that the screw contains an initial short crack, with size a = 3 mm. Thirdly, a screw with a long crack will be modelled, where a = 6 mm.
In all cases, the initial cracks are introduced by defining as initial condition φ = 1. Moreover, the initial crack is vertical, as shown in Fig. 7, has a thickness of 0.16 mm, and its bottom tip is located at a distance of 7 mm to the bottom of the screw. Following Wick et al. (2015), the material properties are taken to be E = 210 GPa, ν = 0.3, = 0.2 mm, and G c = 2.7 N/mm. The crack growth trajectories predicted for the three cases described above are shown in Fig. 8, by plotting the phase field contours. The results agree qualitatively with those obtained by Wick et al. (2015). In the absence of an initial defect, crack nucleation takes place near the head of the screw. This is in agreement with expectations, as the first winding of the thread carries the highest load (see Kristensen et al., 2020b). However, when an initial defect is present, two cracks branch from it and propagate until reaching the sides of the screw. The force versus displacement response is shown in Fig. 9a. In agreement with expectations, the sample without an initial defect is able to carry a larger load. In regard to the screws with an existing defect, the stiffness of the solid is degraded faster in the case of a long crack, relative to the sample with a smaller crack, but the magnitude of the maximum force attained is similar in both cases. The number of iterations required to achieve convergence is shown for every load increment in Figs. 9b-9d for, respectively, the case without an initial defect, the case with an initial long crack and the case with an initial short crack. In all three cases convergence can be readily attained.
The crack grows in an unstable fashion and the situation thus resembles that of Section 4.1; convergence can be readily attained but one specific increment requires more than 100 iterations to do so.

3D Brazilian test
Finally, we showcase the potential of the framework presented in capturing structural failure in 3D solids. We do so by simulating the fracture of a brittle solid subjected to the Brazilian test. The Brazilian test is a laboratory experiment widely used in the rock mechanics community to indirectly measure the tensile strength of brittle materials. As shown in Fig. 11a, a circular disk is compressed between two jaws until fracture occurs. Upon the assumption that failure occurs at the centre of the disk, closed form expressions can be used to determine the material tensile strength from the remote load (Garcia-Fernandez et al., 2018). As shown in Fig. 11b, The material properties are defined as follows. On the one side, the jaws are typically made of steel, for which E = 210 GPa and ν = 0.3 are assumed.
For the disk we consider a brittle solid with elastic properties E = 25 GPa and ν = 0.2 and fracture properties = 0.5 mm and G c = 0.16 N/mm. The jaws radius to disk radius ratio is chosen to be R j /R d = 1.5. The contact between the jaws and the disk is defined as surface to surface contact with a finite sliding formulation. The normal behaviour is based on a hard contact formulation, where the contact constraint is enforced with a Lagrange multiplier representing the contact pressure in a mixed formulation. The tangential contact behaviour is assumed to be frictionless. To prevent damage under compression, the spectral tension-compression decomposition by Miehe et al. (2010a) is adopted -see Appendix A. Also, an anisotropic formulation is used, such that the strain energy density split is accounted for in the balance equation for the displacement problem (see Appendix A for details).
The results obtained are shown in Fig. 12 in terms of the phase field contours for the different loading stages. The evolution of the phase field is also shown in Video 1, provided in the online version of this manuscript. Subfigures 12 (a)-(c) show in red colour the phase field contours where φ > 0.9.
The crack appears to initiate at the centre of the disk and propagates towards the jaws very fast. Also, smaller cracks nucleate near the loading region.
These calculations have been obtained using 345 load increments and using the monolithic implementation, no convergence issues were observed.

Conclusions
We have presented a simple and robust implementation of the phase field The potential of the framework is demonstrated by addressing four 2D and 3D paradigmatic boundary value problems. First, unstable fracture is examined using a notched square plate subjected to tension. Secondly, stable crack growth is investigated by subjecting the square plate to shear loading. Thirdly, the fracture of screws with and without internal cracks is investigated. Finally, the Brazilian test is simulated, including the modelling of the contact between the jaws and the disk. We observe that the monolithic standard Newton implementation provided is able to reach convergence in all cases. However, a single-pass staggered scheme appears to be more efficient in convergence-wise demanding problems. Computations are efficient but both schemes seem to perform worse than quasi-Newton methods (Kristensen and Martínez-Pañeda, 2020;Wu et al., 2020a). We also find that the use of interpolation schemes might not lead to efficiency improvements in phase field fracture. The framework can be very easily extended to other material models (e.g., plasticity) and damage mechanisms, such as fatigue.

Acknowledgements
The authors would like to acknowledge financial support from the Ministry of Science, Innovation and Universities of Spain through grant PGC2018-099695-B-I00. E. Martínez-Pañeda additionally acknowledges financial support from the EPSRC (grants EP/R010161/1 and EP/R017727/1) and from the Royal Commission for the 1851 Exhibition (RF496/2018).

Appendix A. Additional details of numerical implementation
The framework can be easily extended to incorporate other constitutive choices. Specifically, as shown in the results section, a tension-compression split of the driving force for fracture should be considered to prevent damage from developing under compressive stresses. Alternative strain energy splits are described below, together with an anisotropic phase field formulation where the split is incorporated into the linear momentum equation. All these extensions are implemented in the user material (UMAT) subroutine.
For simplicity, the code accompanying this manuscript (to be downloaded from www.empaneda.com/codes) does not include these additional features, but an extended version can be provided upon request. (2009) volumetric-deviatoric split. In both cases, the strain energy density is decomposed as follows, and only ψ + is considered in the evaluation of the history field H, Eq. (12). In regard to the specific constitutive definition of ψ + 0 , the volumetric-deviatoric split assumes that the compressive part of the volumetric strain energy does not contribute to the fracture process. Accordingly, where K is the bulk modulus, µ is the shear modulus, denote the Macaulay brackets, such that a ± = (a ± |a|)/2, and ε is the deviatoric part of the strain tensor, such that ε = ε − tr(ε)1/3. Here, 1 is the second-order unit tensor.
On the other hand, the spectral decomposition considers, ψ + 0 = 1 2 λ tr ε + 2 + µ tr ε + 2 (A.4) where λ is the first Lamé constant and a spectral decomposition is applied to the strain tensor, such that: where ε I and n I are the principal strains and principal strain directions (with I = 1, 2, 3). The components ε + and ε − are obtained by considering in (A.6) the tensile and compressive principal strains, respectively.

Appendix A.2. Anisotropic formulation
While the majority of the representative results presented are obtained using the hybrid approach proposed by Ambati et al. (2015), we have also extended our implementation to incorporate the so-called anisotropic approach (Miehe et al., 2010a). Thus, the decomposition into tension and compression components is also considered in the field equation for the displacement problem, such that the Cauchy stress (7) would instead read, From an implementation perspective, this translates into a more elaborate computation of the material Jacobian, C = ∂σ/∂ε. Thus, the material behaviour is characterised by the following 4th order elasticity tensor: where H ε is the Heaviside function, such that H ε (x) = 1 for x ≥ 0 or H ε (x) = 0 for x < 0, and J ≡ J ijkl = δ ij δ kl , with δ ij being the Kronecker delta. Also, the projection tensor P + = ∂ ε [ε + (ε)] is computed as (Miehe, 1998) H ε (ε a ) δ ab n ai n aj n bk n bl If ε a = ε b then P + ijkl (A.9) cannot be evaluated. Under such circumstances we replace the term ε a + − ε b + / (ε a − ε b ) with H ε (ε a ).