A Unified Abaqus Implementation of the Phase Field Fracture Method Using Only a User Material Subroutine

Abstract

We present a simple and robust implementation of the phase field fracture method in Abaqus. Unlike previous works, only a user material (UMAT) subroutine is used. This is achieved by exploiting the analogy between the phase field balance equation and heat transfer, which avoids the need for a user element mesh and enables taking advantage of Abaqus’ in-built features. A unified theoretical framework and its implementation are presented, suitable for any arbitrary choice of crack density function and fracture driving force. Specifically, the framework is exemplified with the so-called AT1, AT2 and phase field-cohesive zone models (PF-CZM). Both staggered and monolithic solution schemes are handled. We demonstrate the potential and robustness of this new implementation by addressing several paradigmatic 2D and 3D boundary value problems. The numerical examples show how the current implementation can be used to reproduce numerical and experimental results from the literature, and efficiently capture advanced features such as complex crack trajectories, crack nucleation from arbitrary sites and contact problems. The code developed is made freely available.

1. Introduction

Variational phase field methods for fracture are enjoying a notable success [1,2]. Among many others, applications include shape memory alloys [3], glass laminates [4,5], hydrogen-embrittled alloys [6,7], dynamic fracture [8,9], fiber-reinforced composites [10,11,12,13], functionally graded materials [14,15,16], fatigue crack growth [17,18], and masonry structures [19]. The key to the success of the phase field paradigm in fracture mechanics is arguably three-fold. First, the phase field paradigm can override the computational challenges associated with direct tracking of the evolving solid-crack interface. The interface is made spatially diffuse by using an auxiliary variable, the phase field $\varphi$, which varies smoothly between the solid and crack phases and evolves based on a suitable governing equation. Such a paradigm has also opened new horizons in the modelling of other interfacial problems such as microstructural evolution [20] or corrosion [21]. Secondly, phase field modelling has provided a suitable platform for the simple yet rigorous fracture thermodynamics principles first presented by Griffith [22]. This energy-based approach enables overcoming the issues associated with local approaches based on stress intensity factors, such as the need for ad hoc criteria for determining the crack propagation direction [23,24]. Thirdly, phase field fracture modelling has shown to be very compelling and robust from a computational viewpoint. Advanced fracture features such as complex crack trajectories, crack branching, nucleation, and merging can be captured in arbitrary geometries and dimensions, and on the original finite element mesh (see, e.g., [25,26,27,28]). Also, computations can be conducted in a Backward Euler setting without the convergence issues observed when using other computational fracture methods. One reason behind this robustness is the flexibility introduced by solving the phase field, a damage-like variable, independently from the deformation problem. So-called staggered solution schemes have been presented to exploit this flexibility by computing sequentially the displacement and phase field solutions [29], avoiding computationally demanding phenomena such as snap-backs.

The success of phase field modelling has been, not surprisingly, accompanied by a vast literature devoted to the development of open-source codes and finite element implementations of variational phase field methods for fracture. These works have been aimed at both commercial finite element packages, such as COMSOL [30], and open-source platforms like FEniCS [31]. The development of phase field fracture implementations in the commercial package Abaqus has received particular attention [32,33,34,35,36,37,38], due to its popularity in the solid mechanics community. However, these works require the use of multiple user subroutines, most often including a user element (UEL) subroutine. Abaqus’ in-built elements cannot be employed due to the need for solving for the phase field $\varphi$ as a nodal degree-of-freedom. Having to adopt a user-defined finite element carries multiple limitations; namely post-processing requires the use of a dummy mesh or ad hoc scripts, and most in-built features of Abaqus cannot be exploited, as the software suite is effectively used as a solver. In this work, we overcome these limitations by presenting a new implementation that only requires the use of a user material (UMAT) subroutine. The simple yet robust implementation presented is achieved by taking advantage of the analogy between the phase field evolution equation and heat transfer. This not only greatly simplifies the use of Abaqus for conducting phase field fracture studies but also enables taking advantage of the many in-built features provided by this commercial package. In addition, we present a generalized theoretical and numerical framework that encapsulates what are arguably the three most popular phase field fracture models presented to date: (i) the so-called AT2 model [24], based on the Ambrosio and Tortorelli regularization of the Mumford-Shah functional [39], (ii) the AT1 model [40], which includes an elastic phase in the damage response, and (iii) the phase field-cohesive zone model PF-CZM [41,42], aimed at providing an explicit connection to the material strength. Our framework also includes two strain energy decompositions to prevent damage in compressive states: the spectral split [29] and the volumetric-deviatoric one [43]—both available in the context of anisotropic and hybrid formulations [44]. Moreover, the implementation can use both monolithic and staggered solution schemes, enhancing its robustness. Two example codes are provided with this work (www.empaneda.com/codes), both capable of handling 2D and 3D analyses without any modification. One is a simple 33-line code, which showcases the simplicity of this approach by adopting the most widely used constitutive choices (AT2, no split). The other one is an extended version, with all the features mentioned above, aimed at providing a unified implementation for phase field fracture. To the authors’ knowledge, the present work provides the simplest Abaqus implementation of the phase field fracture method.

The remainder of this manuscript is organised as follows. In Section 2 we provide a generalised formulation for phase field fracture, which can accommodate a myriad of constitutive choices. This is exemplified with the AT2, AT1 and CZ-PFM models. Then, in Section 3, the details of the finite element implementation are presented, including the analogy with heat transfer and the particularities of the Abaqus usage. The potential of the implementation presented is showcased in Section 4, where several boundary value problems of particular interest are addressed. Specifically, (i) a three-point bending test, to compare with the results obtained with other numerical methods; (ii) a concrete single-edge notched beam, to compare with experimental data; (iii) a notched plate with a hole, to simulate complex crack paths, merging and nucleation; and (iv) a 3D gear, where cracking occurs due to contact between the teeth. Finally, concluding remarks are given in Section 5.

2. A Generalised Formulation for Phase Field Fracture

In this section, we formulate our generalised formulation, suitable for arbitrary constitutive choices of crack density function and fracture driving force. Consider an elastic body occupying an arbitrary domain $\mathsf{\Omega }\subset {\mathrm{I}\mathrm{R}}^n$$(n\in [1,2,3\left]\right)$, with an external boundary $\partial \mathsf{\Omega }\subset {\mathrm{I}\mathrm{R}}^{n-1}$ with outwards unit normal $\mathbf{n}$.

2.1. Kinematics

The primary kinematic variables are the displacement field vector $\mathbf{u}$ and the damage phase field $\varphi$. In this work, we limit our attention to small strains and isothermal conditions. Consequently, the strain tensor $\mathbf{\epsilon }$ reads

$$\mathbf{\epsilon }=\frac{1}{2}\left(\nabla {\mathbf{u}}^T+\nabla \mathbf{u}\right).$$

(1)

The nucleation and growth of cracks are described by using a smooth continuous scalar phase field $\varphi \in [0;1]$. The phase field describes the degree of damage, being $\varphi =0$ when the material point is in its intact state and $\varphi =1$ when the material point is fully broken. Since $\varphi$ is smooth and continuous, discrete cracks are represented in a diffuse manner. The smearing of cracks is controlled by a phase field length scale ℓ. The aim of this diffuse representation is to introduce, over a discontinuous surface $\mathsf{\Gamma }$, the following approximation of the fracture energy [24]:

$$\mathsf{\Phi }={\int }_{\mathsf{\Gamma }}{G}_c\mathrm{d}S\approx {\int }_{\mathsf{\Omega }}{G}_c\gamma (\varphi ,\nabla \varphi )\mathrm{d}V,\mathrm{for}\ell \to 0,$$

(2)

where $\gamma$ is the so-called crack surface density functional and $G_c$ is the material toughness [22,45]. This approximation circumvents the need to track discrete crack surfaces, a well-known challenge in computational fracture mechanics.

2.2. Principle of Virtual Work. Balance of Forces

Now, we shall derive the balance equations for the coupled deformation-fracture system using the principle of virtual work. The Cauchy stress $\mathbf{\sigma }$ is introduced, which is work conjugate to the strains $\mathbf{\epsilon }$. Also, a traction $\mathbf{T}$ is defined on the boundary of the solid $\partial \mathsf{\Omega }$, work conjugate to the displacements $\mathbf{u}$. Regarding fracture, we introduce a scalar stress-like quantity $\omega$, which is work conjugate to the phase field $\varphi$, and a phase field micro-stress vector $\mathsf{\xi }$ that is work conjugate to the gradient of the phase field $\nabla \varphi$. The phase field is assumed to be driven solely by the solution to the displacement problem. Thus, no external traction is associated with $\varphi$. In the absence of body forces, the principle of virtual work reads:

$${\int }_{\mathsf{\Omega }}\left\{\mathbf{\sigma }:\delta \mathbf{\epsilon }+\omega \delta \varphi +\mathsf{\xi }·\delta \nabla \varphi \right\}\mathrm{d}V={\int }_{\partial \mathsf{\Omega }}\left(\mathbf{T}·\delta \mathbf{u}\right)\mathrm{d}S$$

(3)

where $\delta$ denotes a virtual quantity. This equation must hold for an arbitrary domain $\mathsf{\Omega }$ and for any kinematically admissible variations of the virtual quantities. Thus, by application of the Gauss divergence theorem, the local force balances are given by:

$$\begin{array}{cc}\hfill & \nabla ·\mathbf{\sigma }=0\hfill \\ \hfill & \nabla ·\mathsf{\xi }-\omega =0\hfill \end{array}\mathrm{in}\mathsf{\Omega },$$

(4)

with natural boundary conditions:

$$\begin{array}{c}\hfill \mathbf{\sigma }·\mathbf{n}=\mathbf{T}\\ \hfill \mathsf{\xi }·\mathbf{n}=0\end{array}\mathrm{on}\partial \mathsf{\Omega }.$$

(5)

2.3. Constitutive Theory

The constitutive theory is presented in a generalised fashion, and the AT1 [40], AT2 [24] and PF-CZM [41,42] models are then derived as special cases. The total potential energy of the solid reads,

$$W\left(\mathbf{\epsilon }\left(\mathbf{u}\right),\varphi ,\nabla \varphi \right)=\psi \left(\mathbf{\epsilon }\left(\mathbf{u}\right),g\left(\varphi \right)\right)+\phi \left(\varphi ,\nabla \varphi \right)$$

(6)

where $\psi$ is the elastic strain energy density and $\phi$ is the fracture energy density. The former diminishes with increasing damage through the degradation function $g\left(\varphi \right)$, which must fulfill the following conditions:

$$g\left(0\right)=1,g\left(1\right)=0,g^{\prime }\left(\varphi \right)\le 0\mathrm{for}0\le \varphi \le 1.$$

(7)

We proceed to formulate the fracture energy density as,

$$\begin{array}{c}\hfill \phi \left(\varphi ,\nabla \varphi \right)=G_c\gamma (\varphi ,\nabla \varphi )=G_c{\displaystyle \frac{1}{4c_w\ell }}\left(w\left(\varphi \right)+{\ell }^2{|\nabla \varphi |}^2\right).\end{array}$$

(8)

where ℓ is the phase field length scale and $w\left(\varphi \right)$ is the geometric crack function. The latter must fulfill the following conditions:

$$w\left(0\right)=0,w\left(1\right)=1,w^{\prime }\left(\varphi \right)\ge 0\mathrm{for}0\le \varphi \le 1.$$

(9)

Also, $c_w$ is a scaling constant, related to the so-called geometric crack function:

$$c_w={\int }_0^1\sqrt{w\left(\zeta \right)}\mathrm{d}\zeta .$$

(10)

Damage is driven by the elastic energy stored in the solid, as characterized by the undamaged elastic strain energy density ${\psi }_0$. To prevent cracking under compressive strain states, the driving force for fracture can be decomposed into active ${\psi }_0^{+}$ and inactive ${\psi }_0^{-}$ parts. Accordingly, the elastic strain energy density can be defined as [46]:

$$\psi \left(\mathbf{\epsilon }\left(\mathbf{u}\right),g\left(\varphi \right)\right)={\psi }^{+}\left(\mathbf{\epsilon }\left(\mathbf{u}\right),\varphi \right)+{\psi }_0^{-}\left(\mathbf{\epsilon }\left(\mathbf{u}\right)\right)=g\left(\varphi \right){\psi }_0^{+}\left(\mathbf{\epsilon }\left(\mathbf{u}\right)\right)+{\psi }_0^{-}\left(\mathbf{\epsilon }\left(\mathbf{u}\right)\right)$$

(11)

Also, damage is an irreversible process: $\dot{\varphi }\ge 0$. To enforce irreversibility, a history field variable $\mathcal{H}$ is introduced, which must satisfy the Karush–Kuhn–Tucker (KKT) conditions:

$${\psi }_0^{+}-\mathcal{H}\le 0,\dot{\mathcal{H}}\ge 0,\dot{\mathcal{H}}({\psi }_0^{+}-\mathcal{H})=0.$$

(12)

Accordingly, for a current time t, over a total time $\tau$, the history field can be defined as,

$$\mathcal{H}={\max }_{t\in [0,\tau ]}{\psi }_0^{+}\left(t\right).$$

(13)

Consequently, the total potential energy of the solid (6) can be re-formulated as,

$$W=g\left(\varphi \right)\mathcal{H}+\frac{G_c}{4c_w}\left(\frac{1}{\ell }w\left(\varphi \right)+\ell {|\nabla \varphi |}^2\right)$$

(14)

Now we proceed to derive, in a generalised fashion, the fracture micro-stress variables $\omega$ and $\mathsf{\xi }$. The scalar micro-stress $\omega$ is defined as:

$$\omega ={\displaystyle \frac{\partial W}{\partial \varphi }}=g^{\prime }\left(\varphi \right)\mathcal{H}+\frac{G_c}{4c_w\ell }{w}^{\prime }\left(\varphi \right),$$

(15)

while the phase field micro-stress vector $\mathsf{\xi }$ reads,

$$\mathsf{\xi }={\displaystyle \frac{\partial W}{\partial \nabla \varphi }}=\frac{\ell }{2c_w}{G}_{\mathrm{c}}\nabla \varphi .$$

(16)

Inserting these into the phase field balance Equation (4b), one reaches the following phase field evolution law:

$$\frac{G_c}{2c_w}\left(\frac{w^{\prime }\left(\varphi \right)}{2\ell }-\ell {\nabla }^2\varphi \right)+g^{\prime }\left(\varphi \right)\mathcal{H}=0$$

(17)

We shall now make specific constitutive choices, particularising the framework to the so-called AT2, AT1 and PF-CZM models.

Degradation function$g\left(\varphi \right)$. Both AT2 and AT1 models were originally formulated using a quadratic degradation function:

$$g\left(\varphi \right)={\left(1-\varphi \right)}^2+\kappa$$

(18)

where $\kappa$ is a small, positive-valued constant that is introduced to prevent ill-conditioning when $\varphi =1$. A value of $\kappa =1\times {10}^{-7}$ is adopted throughout this work. Alternatively, the PF-CZM model typically uses the following degradation function,

$$g\left(\varphi \right)=\frac{{(1-\varphi )}^d}{{(1-\varphi )}^d+a\varphi (1+b\varphi )},$$

(19)

with,

$$a=\frac{4E{G}_c}{\pi \ell f_t^2},$$

(20)

where E denotes Young’s modulus and $f_t$ is the tensile strength of the material. The choices of b and d depend on the softening law employed. Two commonly used softening laws are the linear one, with $b=-0.5$ and $d=2$, and the exponential one, with $b=2^{(5/3)}-3$ and $d=2.5$.

Dissipation function. The dissipation function is governed by the magnitude of w and, consequently, $c_w$. For the AT2 model: $w\left(\varphi \right)={\varphi }^2$ and $c=1/2$. Since $w^{\prime }\left(0\right)=0$, this choice implies a vanishing threshold for damage. An initial, damage-free linear elastic branch is introduced in the AT1 model, with the choices $w\left(\varphi \right)=\varphi$ and $c=2/3$. Finally, in the PF-CZM case we have $w\left(\varphi \right)=2\varphi -{\varphi }^2$ and $c=\pi /4$.

Fracture driving force${\psi }_0^{+}$. The variationally consistent approach, as proposed in the original AT2 model, is often referred to as the isotropic formulation:

$${\psi }_0^{+}\left(\mathbf{\epsilon }\right)=\frac{1}{2}\mathbf{\epsilon }:{\mathit{C}}_0:\mathbf{\epsilon }=\frac{1}{2}\lambda {\mathrm{tr}}^2\left(\mathbf{\epsilon }\right)+\mu \mathrm{tr}\left({\mathbf{\epsilon }}^2\right),{\psi }_0^{-}\left(\mathbf{\epsilon }\right)=0.$$

(21)

where ${\mathit{C}}_0$ is the undamaged elastic stiffness tensor and $\lambda$ and $\mu$ are the Lamé parameters. In the context of the AT1 and AT2 models, damage under compression is prevented by decomposing the strain energy density following typically two approaches. One is the so-called volumetric-deviatoric split, proposed by Amor et al. [43], which reads

$${\psi }_0^{+}\left(\mathbf{\epsilon }\right)=\frac{1}{2}K{\langle \mathrm{tr}\left(\mathbf{\epsilon }\right)\rangle }_{+}^2+\mu \left({\mathbf{\epsilon }}^{\prime }:{\mathbf{\epsilon }}^{\prime }\right),{\psi }_0^{-}\left(\mathbf{\epsilon }\right)=\frac{1}{2}K{\langle \mathrm{tr}\left(\mathbf{\epsilon }\right)\rangle }_{-}^2.$$

(22)

Here, K is the bulk modulus, ${\langle a\rangle }_{\pm }=\left(a\pm \left|a\right|\right)/2$, and ${\mathbf{\epsilon }}^{\prime }=\mathbf{\epsilon }-\mathrm{tr}\left(\mathbf{\epsilon }\right)\mathit{I}/3$. The second one is the so-called spectral decomposition, proposed by Miehe et al. [46], which builds upon the spectral decomposition of the strain tensor ${\mathbf{\epsilon }}^{\pm }={\sum }_{a=1}^3{\langle {\epsilon }_I\rangle }_{\pm }{\mathbf{n}}_I\otimes {\mathbf{n}}_I$, with ${\epsilon }_I$ and ${\mathbf{n}}_I$ being, respectively, the strain principal strains and principal strain directions (with $I=1,2,3$). The strain energy decomposition then reads [46]:

$${\psi }_0^{\pm }\left(\mathbf{\epsilon }\right)=\frac{1}{2}\lambda {\langle \mathrm{tr}\left(\mathbf{\epsilon }\right)\rangle }_{\pm }^2+\mu \mathrm{tr}\left[{\left({\mathbf{\epsilon }}^{\pm }\right)}^2\right]$$

(23)

The split can be applied not only to the phase field balance law but also to the balance of linear momentum. Considering the split only in the phase field balance (17) is typically referred to as the hybrid approach [44]. Alternatively, an anisotropic formulation can be used, such that the damaged version of the stress tensor $\mathbf{\sigma }$ is computed as,

$$\mathbf{\sigma }\left(\mathbf{u},\varphi \right)=g\left(\varphi \right)\frac{\partial {\psi }_0^{+}\left(\mathbf{\epsilon }\right)}{\partial \mathbf{\epsilon }}+\frac{\partial {\psi }_0^{-}\left(\mathbf{\epsilon }\right)}{\partial \mathbf{\epsilon }}.$$

(24)

On the other hand, in the PF-CZM model the driving force for fracture is defined as [41]:

$${\psi }_0^{+}=\frac{{\langle {\sigma }_1\rangle }_{+}^2}{2E},$$

(25)

with the other term of the split being given by,

$${\psi }_0^{-}=\frac{1}{2E}\left[{\sigma }_1{\langle {\sigma }_1\rangle }_{-}+{\sigma }_2^2+{\sigma }_3^2-2\nu \left({\sigma }_2{\sigma }_3+{\sigma }_1{\sigma }_3+{\sigma }_1{\sigma }_2\right)\right],$$

(26)

where $\nu$ is Poisson’s ratio and ${\sigma }_i$ are the principal stresses, with ${\sigma }_1$ being the maximum principal (undamaged) stress. The variational consistency is lost but the failure surface of concrete under dominant tension can be well captured [41]. This formulation is only used with the hybrid approach.

In addition, it is important to note that for the AT1 and PF-CZM models there is a minimum value of the fracture driving force, which we denote as ${\mathcal{H}}_{min}$. This is needed as otherwise $\varphi \le 0$, as can be observed by setting $\varphi =0$ and solving the balance Equation (17). The magnitude of ${\mathcal{H}}_{min}$ is then given by the solution of (17) for $\mathcal{H}$ under $\varphi =0$. For the AT1 case: ${\mathcal{H}}_{min}=3G_c/(16\ell )$; while for the PF-CZM model: ${\mathcal{H}}_{min}=2G_c/(\pi a\ell )=f_t^2/\left(2E\right)$.

3. Finite Element Implementation

We proceed to present our finite element model. The unified phase field fracture theory presented in Section 2 is numerically implemented in Abaqus using only a user material (UMAT) subroutine; i.e., at the integration point level. This is achieved by taking advantage of the similitude between the heat transfer law and the Helmholtz-type phase field balance equation. The analogy between heat transfer and phase field fracture is described in Section 3.1, while the specific details of the Abaqus implementation are given in Section 3.2. The present implementation does not require the coding of residual and stiffness matrix terms; however, these are provided in Appendix A for completeness.

3.1. Heat Transfer Analogy

Consider a solid with thermal conductivity k, specific heat $c_p$ and density $\rho$. In the presence of a heat source r, the evolution of the temperature field T in time t is given by the following balance law:

$$k{\nabla }^{2}T-\rho c_p\frac{\partial T}{\partial t}=-r,$$

(27)

Under steady-state conditions the $\partial T/\partial t$ term vanishes and Equation (27) is reduced to,

$$k{\nabla }^{2}T=-r$$

(28)

Now, rearrange the phase field evolution law (17) as,

$${\nabla }^2\varphi =\frac{g^{\prime }\left(\varphi \right)\mathcal{H}2c_w}{\ell G_c}+\frac{w^{\prime }\left(\varphi \right)}{2{\ell }^2}.$$

(29)

Equations (28) and (29) are analogous upon considering the temperature to be equivalent to the phase field $T\equiv \varphi$, assuming a unit thermal conductivity $k=1$, and defining the following heat flux due to internal heat generation,

$$r=-\frac{g^{\prime }\left(\varphi \right)\mathcal{H}2c_w}{\ell G_c}-\frac{w^{\prime }\left(\varphi \right)}{2{\ell }^2}.$$

(30)

Finally, we also define the rate of change of heat flux (r) with temperature ($T\equiv \varphi$),

$$\frac{\partial r}{\partial \varphi }=-\frac{g^{″}\left(\varphi \right)\mathcal{H}2c_w}{\ell G_c}-\frac{w^{″}\left(\varphi \right)}{2{\ell }^2},$$

(31)

as required for the computation of the Jacobian matrix.

3.2. Abaqus Particularities

The analogy between heat transfer and phase field fracture lays the grounds for a straightforward implementation of variational phase field fracture models in Abaqus. Only a user material (UMAT) subroutine is needed, as it is possible to define within the UMAT a volumetric heat generation source (30) and its variation with respect to the temperature (31). It must be noted that a recent version of Abaqus should be used, as the UMAT volumetric heat generation option does not function properly for versions older than 2020. The alternative for versions 2019 or older is to combine the UMAT with a heat flux (HETVAL) subroutine [38].

Abaqus’ in-built displacement-temperature elements can be used, significantly facilitating model development. The same process as for a standard Abaqus model can be followed, with a few exceptions. The user should employ an analysis step of the type coupled temperature-displacement, with a steady-state response. Also, one should define as material properties the thermal conductivity k, the density $\rho$ and the specific heat $c_p$, all of them with a value of unity. To avoid editing the UMAT subroutine, the mechanical and fracture properties are provided as mechanical constants in the user material definition. Also, one should define a zero-temperature initial condition $T(t=0)=0\forall \mathit{x}$. No other pre-processing or post-processing steps are needed, everything can be done within the Abaqus/CAE graphical user interface, and the phase field solution can be visualized by plotting the nodal solution temperature (NT11). Inside of the UMAT, the material Jacobian ${\mathit{C}}_0$ and the Cauchy stress ${\mathbf{\sigma }}_0$ are computed from the strain tensor. The current (undamaged) stress-strain state is used to determine the driving force for fracture, $\mathcal{H}$. Both ${\mathit{C}}_0$ and ${\mathbf{\sigma }}_0$ are degraded using the current value of the phase field $\varphi$ (temperature), which is passed to the subroutine by Abaqus, such that $\mathit{C}=g\left(\varphi \right){\mathit{C}}_0$ and $\mathbf{\sigma }=g\left(\varphi \right){\mathbf{\sigma }}_0$. Finally, $\mathcal{H}$ and $\varphi$ are used to compute r (30) and $\partial r/\partial \varphi$ (31), defined as the volumetric heat generation and its derivative with respect to the temperature. In its simplest form, the code requires only 33 lines.

The implementation also accommodates both monolithic and staggered schemes, enabling convergence even in computationally demanding problems. We choose not to define the non-diagonal, coupling terms of the displacement-phase field stiffness matrix; i.e., ${\mathit{K}}_{\mathit{u}\varphi }={\mathit{K}}_{\varphi \mathit{u}}=\mathbf{0}$. This makes the stiffness matrix symmetric. By default, Abaqus assumes a non-symmetric system for coupled displacement-temperature analyses but one can configure the solver to deal with a symmetric system by using the separated solution technique. The current values of the phase field (temperature) and displacement solutions are provided to the subroutine, so they can used to update the relevant variables (${\mathit{C}}_0$, $\mathbf{\sigma }$, r and $\partial r/\partial \varphi$), such that the deformation and fracture problems are solved in a simultaneous (monolithic) manner. Conversely, one can use solution dependent state variables (SDVs) to store and use the history field of the previous increment ${\mathcal{H}}_t$, effectively freezing its value during the iterative procedure taking place for the current load increment. This is known as a single-pass staggered solution scheme. Although single-pass staggered schemes are very robust, unconditional stability no longer holds and one should conduct a sensitivity analysis to ensure that the load increments employed are sufficiently small. Robustness and unconditional stability can be achieved by using quasi-Newton methods [47,48], but such option is not currently available in Abaqus for coupled temperature-displacement analyses. Independently of the solution scheme, it is known that phase field fracture analyses can achieve convergence after many iterations [48,49]. Thus, the solution controls are modified to enable this (see the example input file provided in www.empaneda.com/codes).

4. Results

We address several paradigmatic boundary value problems to showcase the various features of the implementation, as well as its robustness and potential. First, we use the PF-CZM model to simulate fracture in a three-point bending experiment and compare the results with those obtained by Wells and Sluys [50] using an enriched cohesive zone model. Secondly, we model mixed-mode fracture in a concrete beam to compare the crack trajectories predicted by the AT2 model to those observed experimentally [51]. Thirdly, cracking in a mortar plate with an eccentric hole is simulated to benchmark our predictions with the numerical and experimental results of Ambati et al. [44]. Finally, the AT1 model is used in a 3D analysis of crack nucleation and growth resulting from the interaction between two gears.

4.1. Three-Point Bending Test

First, we follow the work by Wells and Sluys [50] and model the failure of a beam subjected to three-point bending. In their analysis, Wells and Sluys combined the concepts of cohesive zone modelling and partition of unity, using an exponential traction-separation law [50]. To establish a direct comparison, we choose to adopt the so-called phase field-cohesive zone model (PF-CZM) [41,42] using the exponential degradation function.

The geometry, dimensions and boundary conditions are shown in Figure 1a. A vertical displacement of 1.5 mm is applied at the top of the beam, at a horizontal distance of 5 mm to each of the supports. No initial crack is defined in the beam. Following Ref. [50], the mechanical behaviour of the beam is characterized by a Young’s modulus of $E=100$ MPa and a Poisson’s ratio of $\nu =0$, while the fracture behaviour is characterized by a tensile strength of $f_t=1$ MPa and a toughness of $G_c=0.1$ N/mm. Recall that in the PF-CZM model the material strength is explicitly incorporated into the constitutive response and, as a consequence, results become largely insensitive to the choice of phase field length scale, which is here assumed to be $\ell =0.1$ mm. The model is discretised using 4-node coupled temperature-displacement plane strain elements (CPE4T in Abaqus notation). As shown in Figure 1b, the mesh is refined in the center of the beam, where the crack is expected to nucleate and grow. The characteristic element is at least five times smaller than the phase field length scale and the total number of elements equals 5820. Results are computed using the monolithic scheme.

In agreement with expectations and with the results by Wells and Sluys [50], a crack nucleates at the bottom of the beam, in the center of the beam axis. The crack then propagates in a straight manner until reaching the top, as shown in Figure 1c. The resulting force versus displacement response reveals a quantitative agreement with the predictions by Wells and Sluys [50]—see Figure 2.

4.2. Mixed-Mode Fracture of a Single-Edge Notched Concrete Beam

We proceed to model the failure of a concrete beam containing a notch. The aim is to compare the predictions obtained with the AT2 model with the experimental observations by Schalangen [51]. Schalangen subjected a concrete beam to the loading configuration shown in Figure 3. The beam is supported at four locations, and each support is connected to a girder beam through a rod. The cross-sections of the outer rods are smaller than those of the inner rods, to ensure an equal elongation. The load is applied to the center of the girder beams and then transferred through the rods to the concrete beam. The resulting fracture is stable and mixed-mode.

The geometry and boundary conditions of our finite element model aim at mimicking the experimental configuration, see Figure 4a. Two rigid beams are defined, tied to the reference points RP1 and RP2, where the boundary conditions are applied. Both girder beams can rotate around their reference points. The steel rods and supports are modelled and assigned a Young’s modulus $E=210$ GPa and a Poisson’s ratio equal to $\nu =0.3$. The cross-section of the inner rods equals 1000 mm${}^2$ while the cross-section of the outer rods is taken to be ten times smaller, in agreement with the experimental configuration. As shown in Figure 4a, both horizontal and vertical displacements are constrained at the reference point RP1, while RP2 has its horizontal displacement constrained but is subjected to a vertical displacement of 0.5 mm.

Fracture is simulated using the AT2 model. To prevent failure of elements under compression, the strain energy density is divided into tensile and compressive parts employing the strain spectral decomposition proposed by Miehe et al. [29], using the anisotropic formulation (24). The material properties of the concrete beam are taken to be: Young’s modulus $E=35$ GPa, Poisson’s ratio $\nu =0.2$, and toughness $G_c=0.1$ N/mm. The phase field length scale is assumed to be equal to $\ell =2$ mm and, consequently, the characteristic size of the elements along the potential crack propagation region equals 0.5 mm (see Figure 4b). The rods are modelled using truss elements, while the concrete beam is discretised with a total of 28,265 linear quadrilateral coupled temperature-displacement plane strain elements. The results obtained are presented in Figure 5. Both experimental (Figure 5a) and numerical (Figure 5b) results are shown. A very good agreement can be observed, with the crack initiating in both cases at the right corner of the notch and deflecting, following a very similar trajectory, towards the right side of the bottom support.

4.3. Notched Plate with an Eccentric Hole

In this case study, we demonstrate the capabilities of the framework in capturing the interaction of cracks with other defects, and in predicting crack nucleation from arbitrary sites. This is achieved by using the monolithic scheme and without observing convergence issues. Specifically, we chose to model the failure of a mortar plate, which has been experimentally and numerically investigated by Ambati et al. [44]. As shown in Figure 6a, the plate contains a 10 mm notch and an eccentric hole of 10 mm radius. Mimicking the experimental setup, the plate contains two loading pin holes; the bottom one is fixed in both vertical and horizontal directions, while a vertical displacement of 2 mm is applied to the top one. The material properties are $E=5982$ MPa, $\nu =0.22$, $\ell =0.25$ mm and $G_c=2.28$ N/mm. The AT2 phase field model is considered, with no split applied to the strain energy density. We discretise the plate with 56,252 linear plane stress coupled displacement-thermal elements (CPS4T, in Abaqus notation). The characteristic element length in the regions surrounding the notch and the hole is five times smaller than the phase field length scale.

The results obtained, in terms of the crack trajectory, are shown in Figure 6. A very good agreement with the experimental observations is attained (Figure 6b). As shown in Figure 6c, the crack starts from the notch tip and deflects towards the hole. The location of the point of interaction between the hole and the crack originating from the notch appears to be the same for experiments and simulations. Upon increasing the applied load, a new crack eventually nucleates from the right side of the hole, and propagates until reaching the end of the plate. The resulting force versus displacement response is shown in Figure 7, where various images of the crack path have been superimposed to facilitate interpretation. The curve exhibits a linear behaviour until crack nucleation occurs ($u\approx 0.28$ mm), when a sudden drop in the load carrying capacity is observed. The interaction between the crack and the hole induces mixed-mode conditions and crack deflection, which is reflected in the force versus displacement curve. Once the crack has reached the hole, the applied displacement can be further increased without a drop in the load. This is observed until the nucleation of the second crack, which leads to the complete failure of the plate.

4.4. 3D Analysis of Cracking Due to the Contact Interaction between Two Gears

Finally, we proceed to showcase the abilities of the model in simulating complex 3D boundary value problems, involving advanced features such as contact. It should be emphasized that the same subroutine is used for both 2D and 3D analyses as the implementation is conducted at the integration point level. We chose to simulate the nucleation and growth of cracks in the teeth of two interacting gears, a problem of important technological relevance. The geometries of the two gears are shown in Figure 8, with dimensions given in mm. The circular pitch equals 8 mm, the pressure angle is 20${}^{\circ }$ and both the clearance and the backlash equal 0.05 mm. Both gears have a thickness of 3 mm. The boundary conditions are also depicted in Figure 8. The inner hole of each gear is tied to the gear center point. The center of the small, right gear is subjected to a rotation of 1 radian, while a linear rotational spring is considered at the center of the large, left gear. The stiffness of the rotational spring is $7\times {10}^6$ N·mm/rad.

The modelling requires a non-linear geometrical analysis and the use of a contact algorithm to simulate the interaction between the gear teeth. Frictionless contact is assumed for the tangential contact behaviour, which is enforced by making the Lagrangian multiplier equal to zero. The normal contact behaviour is considered to be a hard contact with a surface-to-surface interaction. The penetration of the slave surface into the master surface is minimised under hard contact conditions. The normal contact constraint is enforced through a Lagrangian multiplier. The material properties read $E=210$ GPa, $\nu =0.3$, $\ell =0.25$ mm, and $G_c=2.7$ N/mm. Fracture is predicted using the AT1 model and no split is used for the strain energy density. The model is discretised with more than 120,000 three-dimensional coupled temperature-displacement brick elements. The results obtained are shown in Figure 9, in terms of phase field $\varphi$ contours. Cracking initiates from the root of one of the teeth from the smaller gear and propagates towards the opposite root until the rupture of the gear teeth.

5. Conclusions

We have presented a unified Abaqus implementation of the phase field fracture method. Unlike previous works, our implementation requires only one user subroutine, of the user material type (UMAT). This enables avoiding the use of user elements, with the associated complications in pre- and post-processing, as well as exploiting most Abaqus’ in-built features. The implementation is compact, requiring only 33 lines of code in its simpler form, and can be used indistinctly for 2D and 3D problems. It is also robust, as both staggered and monolithic solution schemes have been incorporated. Moreover, the implementation can accommodate any constitutive choice of phase field model. We present a unified theoretical framework that resembles the code, and particularize it to three of the most widely used phase field models: AT1, AT2 and PF-CZM. In addition, several strain energy splits are considered, in the framework of both hybrid and anisotropic formulations.

We have demonstrated the robustness and capabilities of the framework presented by addressing several boundary value problems of particular interest. First, we showed that the PF-CZM version leads to an excellent agreement with the enriched cohesive zone model analysis by Wells and Sluys [50] of crack nucleation and growth in a beam subjected to three-point bending. Secondly, we validated the crack trajectories predicted by the AT2 model with the experimental observations by Schalangen [51] on a concrete beam exhibiting mixed-mode fracture. Thirdly, we simulated the failure of a mortar plate with an eccentric hole to showcase the capabilities of the framework in capturing the interaction between cracks and other defects, as well as the nucleation of secondary cracks. The simulations agree qualitatively and quantitatively with the results obtained by Ambati et al. [44]. Finally, we used the AT1 version to model cracking due to the interaction between gears to showcase the capabilities of the model in dealing with 3D problems incorporating complex computational features, such as contact and geometric non-linearity. The codes developed have been made freely available, with examples and documentation at www.empaneda.com/codes.