# An open-source Abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation

Gergely Molnar, Anthony Gravouil, Rian Seghir, Julien Réthoré

► To cite this version:

Gergely Molnar, Anthony Gravouil, Rian Seghir, Julien Réthoré. An open-source Abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation. Computer Methods in Applied Mechanics and Engineering, 2020, 10.1016/j.cma.2020.113004. hal-02541808

HAL Id: hal-02541808

https://hal.science/hal-02541808

Submitted on 14 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
An open-source Abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation.

Gergely Molnár<sup>a,*</sup>, Anthony Gravouil<sup>a</sup>, Rian Seghir<sup>b</sup>, Julien Réthoré<sup>b</sup>

**Abstract**

Brittle and ductile dynamic fracture in solids is a complex mechanical phenomenon which attracted much attention from both engineers and scientists due to its technological interests. Modeling cracks in dynamic cases based on a discontinuous description is difficult because it needs additional criteria for branching and widening. Due to the very short time scales and the spatial complexity of the dynamic problem, its experimental analysis is still very difficult. Therefore, researchers still have to rely on numerical simulations to find a deeper explanation for many observed phenomena. This study is set out to investigate the effect of plasticity on dynamic fracture propagation. On the other hand, the diffuse phase-field formulation makes it possible to initiate, propagate, arrest or even branch cracks while satisfying the basic principles of thermodynamics. An implicit, staggered elastoplastic version of the phase-field approach was implemented in the commercial finite element code Abaqus through the UEL option. By means of simple examples we show that localized ductile deformations first increase both resistance and toughness. Then, after a maximum value, the resistance starts to decrease with a significant increment in energy dissipation. By favoring shear deformation over tensile failure the fracture pattern changes. First the branching disappears, then the crack propagation angle changes and becomes a shear band. Finally, we observed the increment of the instantaneous dynamic stress intensity factor during the acceleration stage of the crack without introducing a rate dependent critical fracture energy. We explained this phenomenon with the increasing roughness of the fracture surface.

**Keywords:** Dynamic fracture, Ductile fracture, Phase-field, Abaqus UEL, Dynamic stress intensity factor

**1. Introduction**

Fracture is the most feared failure mode in solids and structures. Both brittle and ductile crack propagation can lead to life threatening and devastating financial consequences. Therefore, a predictive description of the phenomenon is of great importance to both scientists and engineers equally. One of the earliest successful descriptions were proposed by Griffith [1]. The criterion was based on an energy balance equation. Griffith postulated that if the release rate of the strain energy is higher than a critical value, the preexisting crack propagates. The theory became very popular because the critical energy release rate was found to be a material constant for measuring the fracture toughness

<footnote><sup>a</sup> Univ Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France</footnote>

<footnote><sup>b</sup> EC Nantes, CNRS UMR6183, GeM, F-44300, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email address: gergely.molnar@insa-lyon.fr (Gergely Molnár)</footnote>
of brittle materials. The original description did not consider the actual motion and the dynamics of these events. Since then, considerable experimental [2, 3, 4, 5] and theoretical studies [6, 7, 8, 9] were dedicated to understanding the underlying physics and mechanics. However, due to the complexity of the problem, dynamic fracture is still an active area of research.

Beyond a few theoretical cases, practical engineering examples are usually too complex to be solved using analytical techniques. Therefore, numerical methods play a key role in fracture analysis. To model cracks both discrete and continuum methods can be utilized. Molecular dynamics [10, 11] provide a valuable insight into the physics of crack initiation, however neither the length nor the time scales is sufficient for practical applications. Discrete element simulations [12] solve the scaling problem, however their use is mostly limited to highly inhomogeneous materials (e.g. granular mater or rocks). Nevertheless, even if discrete techniques propose some advantages, the continuum description is the most commonly used approach. The finite element method was shown to be able to model dynamic fracture with the help of various enrichments. The cohesive zones method [13, 14] is based on the separation of predefined surfaces assuming a finite process zone located at the front of the fracture. The so-called strong discontinuity approach [15] adds a finite displacement step at the crack front. In more recent years the eXtended finite element method [16] was successfully used to follow dynamic crack propagation in elastic [17] and ductile [18] cases. Generally, branching was identified as an instability and the loss of hyperbolicity [19, 20] of the acoustic matrix. The above mentioned methods represent fracture as discrete discontinuities. Therefore, initiation or branching of cracks were proven to be a complicated task and often needed a special criteria.

Nonlocal damage models [21, 22, 23, 24] provide an alternative description. The thick level set approach determines the crack path based on the damage evolution at the crack tip, and the discontinuity is introduced explicitly. While in other smeared models, such as the phase-field approach, discontinuities are not used directly. Instead, a smooth transition characterizes the effect of the crack. A single damage variable connects the fully damaged region to the intact material. The application of nonlocal damage models to Griffith's theory is based on the variational formulation of Francfort and Marigo [25, 26] and the phase-field implementation of Bourdin [27].

The reversibility of the damage field [28] was first enforced by prohibiting the crack from healing. This formulation leads to a minimization principle with inequality constraints, which is challenging to implement in most finite element methods [29]. Departing from the original idea, Miehe et al. [30] proposed a prevalent approach by introducing a new field variable, the history of the maxima of the elastic (undamaged) energy.

Two main numerical techniques are used in practice to solve the mechanical problem and find the phase-field topology: the monolithic (fully coupled) [31, 32, 33, 34, 35, 36] and the staggered (weakly coupled) [27, 30, 37, 38, 39] scheme. Monolithic algorithms usually converge when the initial estimate is close to the equilibrium solution. However, in the case of e.g., unstable crack initiation, they fail to find a definitive solution. This problem is due to the convexity of the global energy and the jumps in the solution [40]. Various modifications [40, 41, 34] were proposed to circumvent these issues.

The staggered scheme is based on the alternating minimization algorithm of Bourdin et al. [27]. Two well-posed, convex problems are solved independently. While the damage and the elastic strain energy is passed between the two equations only at the beginning of the time step. This iterative process was proven to be very robust, even in unstable cases [38]. The only disadvantage of the method that it requires a relatively small time step [39]. In this work, a common variant is used [30, 42, 43], when the iteration around the two decoupled problems is not performed. In quasi-static cases, for unstable propagation, this could violate the energy conservation of the system. Therefore, we controlled the time step automatically to detect unstable propagation and reduce the change in the external loading. Since the first introduction of the method, numerous open-source implementations
were proposed [32, 44, 39, 33]. For a detailed overview on the history the reader is referred to the work of Bourdin [45].

Since the method was originally introduced, it became surprisingly popular in the fracture science community. Several implementations were proposed to follow dynamic [46, 47, 48, 49, 50, 51, 52] and ductile [53, 54, 55, 56, 57, 58, 59, 60] fracture. Surprisingly, formulations concentrating on the effect of plasticity on the dynamic crack propagation is still scarce. Only two implementations [61, 62] exist with details on the development of the theory and some examples. Despite the benefits, the use of phase-fields is still limited in practice because none of the commonly available finite element software provide it as a built-in option. Recently Liu et al. [63] published a dynamic phase-field model using Abaqus, however their subroutine is still reserved for private use. A smart implementation was proposed recently by Azinpour et al. [64], who used the built-in thermal module of the software and used the temperature field as damage. Furthermore, hydrogen induced cracking [65] and a coupled cohesive zone model [35] was implemented in Abaqus using either UMATs or UELs. Both authors provide an open-source version of their codes.

Crack propagates in an unstable or dynamic manner if the energy release rate increases with crack growth. Therefore, the surplus of released energy can be transformed into kinetic energy. The difference between the fracture toughness and the additional energy release determines how much energy becomes available as kinetic energy and consequently governs the speed at which the crack propagates. Finally, if the energy release exceeds twice the critical value, the crack branches. It is commonly assumed that very close to the crack tip a so-called process zone forms. In this region the material starts dissipating energy (e.g. by plasticity) to regularize the stress singularity. As a result, the macroscopic fracture resistance is higher for more ductile materials. These properties are known to depend upon strain rate. For example at higher deformation rates, the yield strength is higher as well. It was postulated theoretically [2] and observed experimentally [66], that the instantaneous fracture toughness (stress intensity factor in this case) of most materials (even in brittle ones) increases when the crack accelerates. The magnitude of this increment is of course different for brittle and ductile materials. In general, strain rate sensitive materials, such as ferritic steels will show larger variations in the measured stress intensity factor than rate insensitive ones such as high strength steel of aluminum alloys [67]. The uniqueness of the relation between crack velocity and macroscopic fracture toughness is questioned by several authors [68, 66]. However, there is still no theoretical study which explains if this phenomenon is a geometrical effect, or if there should be a unique law to describe the rate dependent fracture surface energy for each material. It would be crucial to clarify this question, as most engineering guidelines [67] and numerical methods (e.g. XFEM [69]) use the law of Kanninen [2] to decide if the crack can propagate and at which speed it runs.

Motivated by the popularity [60] of our recent open-source Abaqus phase-field implementation for static cases [70], we wish to extend the subroutine to dynamic, ductile and coupled cases in both 2D and 3D. The new implementation contains all options and proposes switch variables to the user to decide which option should be used. We implemented the popular staggered time integration scheme [30] for rectangular and triangular (tetrahedral in 3D) elements. Presently, a von Mises [71] type plastic criterion is implemented, which is gradually degraded with damage. This allows us to directly take the plastic energy into account when calculating the fracture topology. Optionally, a penalty function is included if only an ultimate strain value is available for the user.

Through simple examples the validity of the implementation is proved by comparison to analytic solutions and results from literature. Finally, with the help of the newly formulated element we study the impact of plasticity on the fragmentation patterns, the crack velocity and the macroscopic fracture toughness in dynamic crack propagation.

As a summary, we propose a robust and versatile open-source subroutine for practical engineers
as well as research scientists in a widely available commercial finite element code Abaqus. The implementation can be used to study both brittle and ductile static and dynamic fracture.

The paper is structured as follows: Section 2 presents the theory of the elasto-plastic dynamic phase-field formulation. Section 3 asses the effect of the yield parameters in the quasi-static case. While section 4 shows 2D and 3D dynamic crack topologies, and section 5 presents the instantaneous fracture toughness obtained using the Abaqus implementation. Finally, Section 6 concludes the paper.

## 2. Methods

### 2.1. Phase-field approximation of the fracture surface

The principal idea of the phase-field formulation for fracture is the introduction of an internal length scale ($l_c$). Thus, this method diffuses the sharp crack into the volume of the elastoplastic solid as illustrated in Fig. 1.

By introducing damage variable ($d$), a smooth transition is created from undamaged ($d = 0$) to broken material ($d = 1$). Miehe [31] introduced the crack surface density:

$$ \gamma(d, \nabla d) = \frac{1}{2l_c} d^2 + \frac{l_c}{2} |\nabla d|^2. \qquad (1) $$

This approximation uses the spatial gradient of the damage field to describe the crack topology. The volume integral of eq. (2) on the whole body gives the theoretical fracture surface:

$$ \Gamma = \int_{\Omega} \gamma(d, \nabla d) d\Omega. \qquad (2) $$

Our previous work [70] provides further details about the theory of the basic phase-field formulations.

### 2.2. Energy functional

The energy functional of the elastoplastic dynamic problem involves the following Lagrangian function:

$$ \mathcal{L} = D(\dot{\mathbf{u}}) - \Pi(\mathbf{u}, d), \qquad (3) $$

where $D(\dot{\mathbf{u}})$ is the kinetic energy:

$$ D(\dot{\mathbf{u}}) = \frac{1}{2} \int_{\Omega} \dot{\mathbf{u}}^T \dot{\mathbf{u}} \rho d\Omega, \qquad (4) $$

and $\Pi(\mathbf{u}, d)$ is the potential energy:

$$ \Pi(\mathbf{u}, d) = E(\mathbf{u}, d) + P(\mathbf{u}, d) + W(d). \qquad (5) $$

In eq. (4) $\dot{\mathbf{u}}$ contains components of the velocity vector, while $\rho$ is the mass density. The potential energy is constructed from three components: the (i) elastic strain energy ($E$), (ii) plastic strain energy ($P$) and the (iii) fracture energy ($W$). All components depend on either the phase-field ($d$) or the displacement field ($\mathbf{u}$). The following sections are dedicated to discuss the details of theory behind each energy component.
### 2.2.1. Elastic strain energy

Classically the elastic strain energy density is composed of a tensile ($\psi_0^+$) and a compression components ($\psi_0^-$):

$$E(\mathbf{u}, d) = \int_{\Omega} \psi^{el}(\mathbf{u}, d) d\Omega, \qquad (6)$$

$$\psi^{el}(\mathbf{u}, d) = g(d) \psi_0^+(\boldsymbol{\varepsilon}^{el}(\mathbf{u})) + \psi_0^-(\boldsymbol{\varepsilon}^{el}(\mathbf{u})). \qquad (7)$$

The decomposition of the elastic strain energy [24] is meant to account for degradation upon tension, while the material remains intact upon compression. $g(d)$ is the degradation function:

$$g(d) = (1-d)^2 + k, \qquad (8)$$

where $k$ is a small $10^{-8}$ number responsible for numerical stability.

Depending on the implementation, the definition of the positive and negative energy parts can vary. In our case for brittle materials, we consider:

$$\psi^{el}(\varepsilon, d) = \mu \sum_{i=1}^{3} \left[ \langle \varepsilon_i \rangle_{-}^2 + g(d) \langle \varepsilon_i \rangle_{+}^2 \right] + \frac{\lambda}{2} \left[ g(d) \langle tr(\varepsilon) \rangle_{+}^2 + \langle tr(\varepsilon) \rangle_{-}^2 \right]. \qquad (9)$$

While for the ductile case, we assume that the material is damaged equally under the shear in compression and extension [37, 72]:

$$\psi^{el}(\varepsilon, d) = g(d) \left[ \mu \sum_{i=1}^{3} \varepsilon_i^2 + \frac{\lambda}{2} \langle tr(\varepsilon) \rangle_{+}^2 \right] + \frac{\lambda}{2} \langle tr(\varepsilon) \rangle_{-}^2. \qquad (10)$$

The overall strain field is divided into elastic ($\boldsymbol{\varepsilon}^{el}$) and plastic parts ($\boldsymbol{\varepsilon}^{pl}$):

$$\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^{el} + \boldsymbol{\varepsilon}^{pl}. \qquad (11)$$

While assuming small deformations, where the strains are defined by the symmetric part of the deformation gradient:

$$\boldsymbol{\varepsilon} = \nabla^S \mathbf{u}. \qquad (12)$$

Details considering the elastic stress and stiffness calculation can be found in Appendix A.

We note that in the ductile case, if the shear is applied under severe hydrostatic compression, volumetric locking can spoil the final results [73]. As in this paper, we use fracture cases mostly in extension we did not have to deal with the phenomenon published recently for purely shear driven cracks. However, if the end-user wants to apply the program for compression fracture cases, this should be tackled appropriately [74, 75, 56].

### 2.2.2. Plastic strain energy

The plastic energy can be expressed as:

$$P(\mathbf{u}, d) = \int_{\Omega} g(d) \psi_{0}^{pl}(\boldsymbol{\varepsilon}^{pl}(\mathbf{u})) d\Omega = \int_{\Omega} \int_{t} (\dot{\boldsymbol{\varepsilon}}^{pl} : \boldsymbol{\sigma}) dt d\Omega, \qquad (13)$$
where $\sigma$ is the damage degraded Cauchy stress tensor (see eq. (A.3)). Present implementation supports the von Mises [71] yield criterion:

$$f(\sigma) = \sigma_{eq} - g(d) [\sigma^{\text{lim}} + H\varepsilon_{eq}^{\text{pl}}], \quad (14)$$

with $\sigma_{eq}$ as the von Mises stress (calculated equally from the degraded Cauchy stress). As there is no permanent volume change, the ductile part of the energy history can be expressed as a function of the yield stress and the energy equivalent plastic shear strain [76]:

$$\psi_0^{pl}(\varepsilon_{eq}^{pl}(\mathbf{u})) = \varepsilon_{eq}^{pl}(\mathbf{u}) \left[ \sigma^{\text{lim}} + \frac{1}{2} H \varepsilon_{eq}^{pl}(\mathbf{u}) \right] + \frac{1}{2} I_{\varepsilon} \langle \varepsilon_{eq}^{pl} - \varepsilon_{eq}^{cr} \rangle^2, \quad (15)$$

where $\sigma^{\text{lim}}$ is the von Mises yield strength, $H$ is the hardening modulus and $\varepsilon_{eq}^{pl}$ is the energy equivalent plastic shear strain. The last term in eq. (15) with the positive part function<fnref n="1" /> is a penalty component which assures that the material breaks at a given critical strain. $I_{\varepsilon}$ is chosen according to the considerations discussed in Appendix C.

It can be noticed that with the evolution of the damage the yield strength of the material gradually becomes zero, allowing plasticity to initiate and guide fracture. Additionally, the UEL option allows the user to implement any desired yield criterion. Current formulation allows us to treat the plastic problem as a classical von Mises type model, for which details can be found in various textbooks [77, 76].

If the plastic option is active, plastic deformation under axial compression and tension contributes the same way to damage. Thus, the switch variables in the energy decomposition are $\alpha_{i=1,2,3} = 1$ independently of the sign of the principal strains. Details can be found in Appendix A. An anisotropic Hill yield criterion [78] could be an alternative treatment for this problem. Furthermore, we consider the effect of stress triaxiality [62] on the yield function essential to reproduce the experimentally observed crack nucleation in plastic materials. However, the implementation of such approach is beyond the scope of present paper.

### 2.2.3. Fracture energy

The main idea of the phase-field approach is that the discontinuity is smeared and treated as a continuous field between damaged and intact materials. Therefore, the fracture energy can be expressed from eq. (2) as a function of the damage variable (d) itself:

$$W^0(d) = \int_{\Gamma} g_c d\Gamma \approx \int_{\Omega} g_c \gamma(d, \nabla d) d\Omega = \int_{\Omega} \frac{g_c}{2l_c} [d^2 + l_c^2 |\nabla d|^2] d\Omega, \quad (16)$$

where $\gamma$ is the fracture energy density, $g_c$ is the surface energy needed to create a unit fracture surface and $l_c$ is the length scale parameter which measures the scale of the damage diffusion.

Interestingly, a slightly different fracture energy density was proposed by Fremond [79] which provides an elastic threshold:

$$W(d) = \int_{\Omega} \psi_c [2d + l_c^2 |\nabla d|^2] d\Omega. \quad (17)$$

Note that compared to eq. (16), in eq. (17) d enters the formulation by a linear term.

<footnote n="1">¹$\langle \cdot \rangle = (\cdot + |\cdot|)/2$</footnote>
After certain algebraic manipulations Miehe [55] achieved a very similar form to their original formulation:

$$W(d) = \int_{\Omega} t(d) + w^{0}(d) d\Omega, \qquad (18)$$

where $t(d)$ will act as a threshold energy in the staggered formulation (see later):

$$t(d) = \psi_c - g(d) \psi_c. \qquad (19)$$

$\psi_c = g_c/2l_c$ was chosen, as a result $w^0$ became identical to eq. (16):

$$w^0(d) = 2l_c \psi_c \gamma(d, \nabla d) = \frac{g_c}{2l_c} [d^2 + l_c^2 |\nabla d|^2] \qquad (20)$$

This equation remains fixed. However, to switch between the classic solution [30] and the elastic threshold the user can set $\psi_c$ to zero in $t(d)$. Thus omitting the initial change in the energy history.

### 2.3. Staggered time-integration algorithm

There is an ongoing discussion whether the monolithic or the staggered scheme should be used in practice. In the case of unstable crack propagation, the monolithic solution [31] tends to become numerically unstable as the algorithm needs to find the fully formed crack path in only one iteration, which is of course impossible. There are various techniques, such as viscous regularization [40], path following algorithms [41], or even globalization [34], which can improve the solution. However, none of them are suitable for the code implemented in Abaqus. Therefore, present the paper focuses on the staggered scheme with the history variable proposed by Miehe et al. [30]. This method solved the topological problem and provided incredible robustness to the technique.

The basic idea is that the displacement and the phase-field problems are coupled only weakly. In each iteration they are independent and interact only through a so-called history field. In Fig. 1 the schematic illustration of the staggered algorithm is shown. The sharp crack is regularized by the phase-field which is calculated based on the energy history. Then, the damage field is used to recalculate the displacement distribution.

The energy of the displacement problem is formulated as follows:

$$\Pi^u = D(\dot{\mathbf{u}}) - E(\mathbf{u}, d) - P(\mathbf{u}, d) + \Pi^{ext} \qquad (21)$$

where $\Pi^{ext}$ is the external work done by the body ($\bar{\gamma}$) and boundary ($\bar{\mathbf{t}}$) forces:

$$\Pi^{ext} = \int_{\Omega} \bar{\gamma} \cdot \mathbf{u} dV - \int_{\partial \Omega} \bar{\mathbf{t}} \cdot \mathbf{u} dA. \qquad (22)$$

By taking the variation of eq. (21), the corresponding strong form can be obtained:

$$\begin{align}
\delta\Pi^u &= 0 \quad \forall \delta\mathbf{u} \rightarrow \nabla\sigma - \bar{\gamma} = \rho\bar{\mathbf{u}} \quad \text{in } \Omega \nonumber \\
\sigma \cdot \mathbf{n} &= \bar{\mathbf{t}} \quad \text{on } \Gamma_N \nonumber \\
\mathbf{u} &= \bar{\mathbf{u}} \quad \text{on } \Gamma_D. \nonumber
\end{align} \qquad (23)$$

This is then solved for $\mathbf{u}$, assuming $d$ is constant.

Similarly, the Lagrangian equation of the phase-field problem is written as:

$$\Pi^d = \int_{\Omega} [g_c \gamma(d, \nabla d) + g(d) \mathcal{H}] d\Omega, \qquad (24)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of the staggered scheme for solving the phase-field problem in elastoplastic solids. $l_c$ controls the approximate width of the damaged zone.</figcaption></figure>

where the potential energy from the displacement problem is replaced by the history field<fnref n="2" />:

$$ \mathcal{H}_0 = 0 \\ \mathcal{H}_{n+1} = \max \left\{ \psi_0^+ + \psi_0^{pl} - \psi_c \quad \mathcal{H}_n \right\}. \qquad (25) $$

Furthermore, this formulation enforces the irreversibility of the damage ($\dot{d} \ge 0$). Thus, the history field satisfies the Karush-Kuhn-Tucker conditions. Finally, the corresponding strong form can be expressed as:

$$ \begin{array}{l} \delta\Pi^d = 0 \quad \forall \delta d \rightarrow \frac{g_c}{l_c} (d - l_c^2 \Delta d) = 2(1-d) \mathcal{H} \quad \text{in } \Omega \\ \nabla d \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_d. \end{array} \qquad (26) $$

Due to the limitations of the UEL option of Abaqus, the two problems are solved at the same time independently based on the variables obtained in the previous iteration. In Fig. 2 a flowchart shows the basic iteration process.

The staggered scheme is implemented so that the two elements are connected through only the common block. For dynamic cases Abaqus uses the Hilber-Hughes-Taylor (HHT) [80] method to obtain equilibrium. HHT solves the linearized equilibrium problem applying the following Newton-Raphson iteration:

$$ \left[ \begin{array}{cc} \mathbf{S}_n^u & \mathbf{0} \\ \mathbf{0} & (1+\alpha) \mathbf{K}_n^d \end{array} \right] \left[ \begin{array}{c} \mathbf{u}_{n+\Delta t} \\ \mathbf{d}_{n+\Delta t} \end{array} \right] = - \left[ \begin{array}{c} \alpha \mathbf{r}_n^u \\ \alpha \mathbf{r}_{n-1}^d - (1+\alpha) \mathbf{r}_n^d \end{array} \right], \qquad (27) $$

where $\mathbf{r}_n^u$ and $\mathbf{r}_n^d$ are the residues, $\mathbf{K}_n^u$ and $\mathbf{K}_n^d$ are the elementary stiffness matrices of the displacement and phase-field problems at time $t_n$. $\mathbf{u}_{n+\Delta t}$ and $\mathbf{d}_{n+\Delta t}$ are the new nodal solutions at

<footnote n="2">²If the energy decomposition is not used, both tension ($\psi_0^+$) and compression ($\psi_0^-$) energy densities contribute to $\mathcal{H}$.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Flowchart of the staggered solution used to implement the coupled displacement phase-field solution in Abaqus.</figcaption></figure>

$t_n + \Delta t$.

For details about the static calculations, the reader is referred to our previous paper [70].

In the dynamic case the displacement residue contains an inertial component as well ($\mathbf{f}_n^{\text{ine}}$) and can be expressed as:

$$ \mathbf{r}_n^u = (1 + \alpha) \mathbf{f}_n^{\text{int}} - \alpha \mathbf{f}_{n-1}^{\text{int}} + \mathbf{f}_n^{\text{ine}} - \mathbf{f}_n^{\text{ext}}, \quad (28) $$

where $\mathbf{f}_n^{\text{ext}}$ is the external and $\mathbf{f}_n^{\text{int}}$ is the internal force vector. Parameter $\alpha$ is a damping coefficient, which by default is set to $\alpha = -0.05$.

The tangent matrix of the displacement problem contains both materials stiffness ($\mathbf{K}_n^u$) and mass ($\mathbf{M}$) matrices:

$$ \mathbf{S}_n^u = \mathbf{M} \frac{d\ddot{\mathbf{u}}}{du} + (1 + \alpha) \mathbf{K}_n^u, \quad (29) $$

where the acceleration derivative is assumed to take the following form with $\beta = (1 - \alpha)^2/4$:

$$ \frac{d\ddot{\mathbf{u}}}{du} = \frac{1}{\beta\Delta t^2}. \quad (30) $$

All the corresponding residue vectors and stiffness matrices can be found either in Appendix B or in our previous paper [70].

### 3. Effect of yield parameters

In this section, through simple static examples the effect of each yield parameter will be explained. Our aim is to demonstrate, how plasticity can dominate the fracture process and sometimes alter the crack path.

#### 3.1. Yield strength

A single 2D plane strain element is the simplest case, where the phase-field model can be understood. A rectangular plate with dimensions of 1 mm × 1 mm is subjected to simple shear by moving its top side in the x direction while constraining its bottom. The shear modulus was set to $\mu = 80.77$ GPa with $H = 5$ GPa, $g_c = 5 \cdot 10^{-3}$ J/mm², $l_c = 0.1$ mm and $I_e = 0$. Fig. 3 shows the corresponding von Mises stress as a function of the applied strain for different yield strength values.
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Equivalent stress-strain curves for elastic and plastic solutions for different yield strength values. (b) Phase-field damage as a function of equivalent strain for elastic solution and different yield strength values.</figcaption></figure>

It can be seen, that the brittle case (blue line) is elastic until the maximum stress. This is due to the elastic threshold ($\psi_c$) introduced in the fracture energy description.

When the yield strength is reduced, the material plastifies and the stress increases with a reduced slope ($\mu_t$) related to the hardening constant $H$:

$$H = \frac{\mu_t}{1 - \mu_t/\mu}, \qquad (31)$$

where $\mu$ is the elastic shear modulus. When reversing the deformation and relaxing the stresses the stiffness of the material is unaffected in the plastic stage, but starts to decrease when the material sustains damage.

## 3.2. Damage enhanced plastic deformation

By applying the damage function directly on the plastic energy, we reduce the yield strength gradually. This creates a smooth transition from elastoplastic to fully damaged response. However, during this transformation additional plastic deformation can occur. In order to demonstrate this phenomenon, consider the same homogeneous problem as shown in section 3.1 with a yield strength of $\sigma^{lim} = 5$ GPa ($H = 0$). In Fig. 4 the stress strain curves are shown for analytic (a) and staggered solutions (b). In part (a) it can be seen that the yield strength is initially higher than the analytic elastic limit stress ($\sigma_{eq} = 3.47$ GPa). However, when the material is damaged the relative order of the two curves change. If the damaged elastic stress is higher then the damaged yield strength plastic deformations is induced.

The theoretical reason behind the phenomenon can be explained by writing the yield criterion as:

$$\sigma_{eq} = g(d) \sigma_{eq}^{el} = g(d) 3\mu\epsilon_{eq} \leq g(d) \sigma^{lim}. \qquad (32)$$
<figure><img src="image_5.png" /><figcaption>Figure 4: Shear stress and yield strength as a function of applied deformation: (a) analytic solution, (b) staggered solution with relatively large strain increments.</figcaption></figure>

This equation can be simplified by $g(d)$, and it can be shown that plastic deformation becomes independent of damage. Thus, when the elastic (undamaged) stress crosses the initial yield strength, plasticity takes over. The macroscopic response will be affected slightly due to this phenomenon as the degradation function is quadratic and the absolute difference is negligible compared to the maximum stress value. Nevertheless, caution is advised when analyzing plastic strain fields around the crack path.

The increased plastic deformation originates from comparing a local (stress based plasticity) to a nonlocal quantity (damage). To improve the description Miehe [81] used a gradient plasticity model, where the equivalent plastic multiplier was treated as an independent degree of freedom. However, in our implementation, to remain general and to allow the end users to implement any desired yield criterion, we kept the local description of the plastic material model.

Fig. 4(b) shows the staggered solution of the same problem. The strain steps are intentionally enlarged to demonstrate the issue at hand. When they are decreased, the stress result converges to the analytic one. Due to the solution scheme explained in Fig. 2, the damage is delayed by 2 time steps. This way, the elastic stress can sometimes overcome the theoretical damage threshold. Evidently, to avoid this problem it is recommended to perform a convergence analysis on the time step size, however there are cases when the reduction of the macroscopic time step will not affect the local strain increments. One well known example is the departure of an unstable crack propagation.

### 3.3. *Complex static problem*

To demonstrate the effect of the yield strength on the static fracture topology the classic single edge notched specimen is used. In Fig. 5a the geometry of the model is depicted. The bottom side is constrained in both $x$ and $y$ directions, while the top is extended in the $y$ direction maintaining $u_x = 0$. The Young's modulus was set to 210 GPa with a Poisson's ratio of 0.3. No hardening was applied, and the penalty term was also omitted ($I_\epsilon = 0$). The fracture properties were set to
<figure><img src="image_6.png" /><figcaption>Figure 5: Single edge notched specimen: (a) geometry; (b) horizontal fracture for high yield strength; (c) shear band for low yield strength cases; (d) reaction force as a function of the top side displacement.</figcaption></figure>

$g_c = 2.7 \cdot 10^{-3}$ J/mm$^2$ and $l_c = 6 \cdot 10^{-3}$ mm. The right hand side of the model, where the fracture is expected, was meshed in an unstructured manner with $h_{FEM} = 3 \cdot 10^{-3}$ mm size elements.

To control the time step an automatic integration scheme is used. The increment of the history field was limited to $d\mathcal{H} = \psi_c$. However, in case of an unstable crack initiation the increment of $\mathcal{H}$ becomes independent of the time step. When this phenomenon is detected, the step is limited to a minimum of $\Delta u_y^{min} = 10^{-12}$ mm. This way, the unstable crack propagation is identified precisely. However, we note that this problem is not present in dynamic cases.

Fig. 5(b) and (c) show the fracture pattern as a function of the yield strength. It can be seen, that if $\sigma^{lim}$ is sufficiently low, the fracture favors a tilted angle and the propagation is led by plastic shear deformation. This kind of shear localization is called shear banding [82].

In Fig. 5(d) the reaction force is shown as a function of the top side's displacement. It can be seen, that for high yield strength cases (brittle, $\sigma^{lim} = 3,4$ GPa) the crack initiates and propagates in an unstable manner. While for lower $\sigma^{lim}$ values it remains stable and a plateau is present. These results show how the addition of plasticity increases the fracture toughness of the configuration. Interestingly, by adding a small amount of plastic yield ($\sigma^{lim} = 4$ GPa) not only the toughness but the resistance of the material becomes higher. However, by further reducing $\sigma^{lim}$ the resistance decreases significantly with the increment of the energy dissipation.

## 4. Fracture pattern

To determine the effect of plasticity on the dynamic crack propagation patterns, as well as to validate our implementation, two well-known benchmark examples are used. First, a 2D plate is loaded with an instantaneous traction stress on the top and bottom boundaries. Secondly, a 3D slab is pre-extended statically (without damage), then the dynamic fracture is followed using the phase-field method.

### 4.1. 2D dynamic branching

The geometry of the 2D test is shown in Fig. 6(a). Two different finite element tessellations were tested: (i) an irregular one shown in Fig. 6(b) with ~92 000 elements, where the right hand side of the model is densified randomly with $h_{FEM} = 0.125$ mm size rectangular elements; and (ii) a structured
<figure><img src="image_7.png" /><figcaption>Figure 6: Rectangular plate subjected to uniform unidirectional dynamic traction: (a) Geometry with dimensions in mm; (b) random finite element mesh with a refined zone ($h_{FEM} = 0.125$ mm) in the expected propagation zone.</figcaption></figure>

one, where all sides are partitioned equally with the same $h_{FEM}$, creating a 256 000 element model. The Young's modulus is set to 32 GPa with a Poisson's ratio of 0.2. The density is 2450 kg/m$^2$, with fracture properties $g_c = 3$ J/m$^2$ and $l_c = 0.5$ mm. Neither hardening ($H = 0$), nor strain penalty ($I_e = 0$) is applied. The time step was set to $\Delta t = 10^{-8}$ s.

In the following sections the obtained results are discussed as a function of the yield strength: (i) fracture topology, (ii) crack velocity and finally in section 5 (iii) the instantaneous stress intensity factor.

First, to verify our implementation, we compared our fracture pattern to ones obtained in Ref. [48]. When the elastic threshold was suppressed ($\psi_c = 0$) the fracture patterns were in agreement. However, when the elastic threshold option was used, the crack started to propagate later as well as the branching occurred after a longer crack path. This is due to the fact that the material became more resistant with $\psi_c$. Fig. 7(a) and (b) also shows that in the elastic case the crack is confined, while in Borden's work it can diffuse more into the material.

Our next comparison will focus on the difference between structured and unstructured mesh. It can be seen in Fig. 7(b) and (c) that the mesh has a significant effect on the branching. Unfortunately, in dynamic calculation the wave propagation is significantly affected by the uniformity of the mesh [83]. Thus, when an instability occurs (like a branching) it will be influenced by the symmetry of the velocity field. Note, that this phenomenon was already observed in static calculations [70].

Finally, the effect of the yield strength can be seen in Fig. 7(c-e). By reducing the yield strength the branching disappears as the local kinetic energy is being consumed by plasticity. Then, if $\sigma^{lim}$ is decreased further, similarly to the static case, the crack starts to propagate in the direction of shear rather than the tensile mode.

Fig. 7(f) shows the equivalent plastic shear strain along the shear band in the case of $\sigma^{lim} = 2$ GPa. It can be seen that the plastic deformation has localized around the damaged zone to a band-like
shape. However, as highlighted in section 3.2, the actual value of the plastic strain should be treated with caution because of the enhancing effect of the phase-field formulation. The leaking shown in part (e) is present in almost all dynamic phase-field models if the simulation is conducted long enough, as the kinetic energy is not entirely dissipated when the crack reaches the boundaries. Therefore, the remaining waves concentrate usually at the damaged outer region.

For quantitative comparison, the crack tip velocity is shown in Fig. 8. The position of the crack tip was calculated as the largest extension of the phase-field where its value is larger then $d > 0.9$ in the lower half of the model ($y < 0$). An agreement was found between our work and Borden’s [48]. It can be seen, that adding an elastic threshold the fracture initiates and branches later, but reaches the same maximum velocity. This velocity limit is in good correspondence with the theoretical limit [84] related to the Rayleigh wave speed ($v_R$).

It can be observed that the velocity at the initiation of the crack slows down by reducing the plastic yield strength which is mostly due to the loss of rigidity around the crack tip. For $\sigma^{lim} = 2$ GPa, when a shear localization occurs, the maximum velocity cannot really be defined as damage develops at the same time along the shear band.

## 4.2. Effect of plane stress state

All techniques which are used to distinguish between tensile and compression energies are based on the decomposition of a 3D stiffness tensor. While by omitting the third directional strain this can be simplified for plane strain cases, in the plane stress state $\epsilon_z$ is not zero. Therefore, in elasticity the components of the stiffness matrix is modified to reproduce a $\sigma_z = 0$ state.

Unfortunately, the energy decomposition routines used for phase-field simulations are not suited for plane stress cases. Therefore, to analyze the effects of $\sigma_z$ being zero, a simple 3D trick is applied. The same geometry is created as shown in Fig. 6, however the plate is extruded in the $z$ direction and the thickness is meshed with one element. To achieve the plane stress state all the nodes were left free in the perpendicular ($z$) direction.

In Fig. 9 the crack tip velocity is shown for different yield strength as a function of time. The same phenomenon is visible as in plane strain cases, the crack tip slows down as the yield strength decreases. For the brittle material we did not observe a significant difference between the plane strain and plane stress cases. Unfortunately, neither the temporal nor the spatial resolution of the simulation is high enough to capture the difference caused by the 2 % decrease in the dilatational wave speed.

On the other hand, the fracture pattern shows a significant difference for the ductile cases. In the plane stress case the straight fracture (as seen in Fig. 7(d)) never appears. After the branching, the shear band forms between the main branch and the side of the specimen. Finally, at $\sigma^{lim} = 1.5$ MPa the fracture pattern is identical to Fig. 7(e).

## 4.3. Pre-stressed 3D slab

In the previous two cases, when an instantaneous traction is applied on the two opposite surfaces, the shock wave traveled in a relatively homogeneous manner, therefore the 3D effect of the crack tip [67] was not visible.

Recently, Henry & Adda-Bedia [50] and Bleyer & Molinari [85] showed that initiation, propagation and branching of a crack is a truly 3D phenomenon. In this paper our aim is to partially reproduce their results in order to validate our implementation, however due to the demanding computational cost we concentrate on simple examples where the asymmetric energy degradation has a significant effect.
<figure><img src="image_8.png" /><figcaption>Figure 7: Fracture pattern as a function of the yield strength and mesh uniformity. Comparison of part (a) should be made with the work of Borden et al. [48].</figcaption></figure>
<figure><img src="image_9.png" /></figure>

Figure 8: Crack tip velocity as a function of time for different yield strength. For the purpose of verification with dashed line a threshold free ($\psi_c = 0$) data is compared to the work of Borden et al. [48]. The crack velocity was calculated by taking the time derivative of the positions of the crack tip. The position of the crack tip was calculated as the largest extension of the phase-field where its value is larger than $d > 0.9$.

<figure><img src="image_10.png" /></figure>

Figure 9: Crack tip velocity as a function of time for different yield strength under plane stress conditions. The figure shows the fully formed fracture patterns for each case, except for the brittle one, which was recorded at $t = 100$ µs.
<figure><img src="image_11.png" /><figcaption>Figure 10: Geometry and boundary conditions of the 3D tensile test.</figcaption></figure>

Two models with different $W/l_c$ ratios are considered: 5 and 50. We note here, that Ref. [50] showed that the undulating fracture pattern appears at a higher value with the addition of local imperfections to the crack tip.

Fig. 10 shows the geometry of the model. The slab was first extended ($\Delta\hat{u}_y$) in the $y$ direction without any constraints in the $x$ or the $z$ translations. This was done statically without damage (all damage degrees of freedom were set to zero in the mechanical simulation). Then we switched to a dynamic calculation and the fracture propagation was followed using the staggered phase-field approach. The Young's modulus was set to 3.09 GPa with a Poisson's ratio of 0.35. The penalty term in the plastic energy was omitted. The fracture properties were set to $g_c = 300 \cdot 10^{-6} \text{ J/mm}^2$ and $l_c$ was varied between 0.2 mm and 0.03 mm for $W/l_c = 5$ and 50 respectively. A finite element size of $h_{FEM} = l_c/2$ was used with a randomly generated mesh.

In Fig. 11(a), the fracture pattern is shown for the brittle case. The micro branching instabilities compare well to the results showed in Ref. [85]. However, when switching to plasticity (even with an infinite $\sigma^{lim}$), the pattern stabilizes. This is due to the difference in the energy degradation. Present implementation uses isotropic degradation ($\alpha_i = 1$) if the plastic option is used. This result is consistent with the fact that the instabilities were not present in the 2D simulation of Bleyer & Molinari [86] either. Furthermore, when the plastic yield strength is reduced a localization appears at the crack tip due to the pre-plastified state. This can be reduced by increasing the hardening coefficient.

These results show that the complete anisotropic degradation function is necessary to reproduce the local instabilities. This problem raises a question: is it valid to apply the degradation function on the plastic energy in eq. (14)?

By increasing the $W/l_c$ ratio, more details can be seen. Fig. 12 shows the evolution of the fracture pattern as a function of time for the brittle case. It is interesting to see that the crack starts to propagate with a microscopic branching which then coalesces creating a chip. Then, it propagates
<figure><img src="image_12.png" /><figcaption>Figure 11: Fracture pattern for the 3D large scale test with different formulations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Fracture pattern for the 3D small scale test.</figcaption></figure>

straight until it branches macroscopically. This phenomenon is indeed due to the effect of the local elastic contraction at the crack tip [67].

After the analysis of the fracture patterns, the instantaneous stress intensity factor is calculated to investigate the issue of the velocity toughening mechanism of the crack propagation.

## 5. Instantaneous fracture toughness

The correlation between dynamic (instantaneous) fracture toughness ($K_{ID}$) and crack tip velocity ($v_c$) is a well documented phenomenon [66]. It was shown that most materials exhibit an increment in $K_{ID}$ when the crack accelerates. Kanninen [2] proposed the following a generalized description:

$$K_{ID}(v_c) = \frac{K_{IC}}{1 - \left(\frac{v_{cr}}{v_l}\right)^m}, \qquad (33)$$

where $K_{IC}$ is the static fracture toughness in mode I opening, $v_l$ is the maximum propagation velocity and $m$ is a material constant. It was shown that not only the material but also the geometry of the experiment has a significant effect on the observed results [66]. Therefore, the relationship can hardly be considered a unique material property. This is a significant problem, because when using other simulation techniques (e.g. XFEM [69]) the kinetics of the crack is determined by the above relationship which is assumed to be a material property.
<figure><img src="image_14.png" /><figcaption>Figure 13: Effective fracture toughness measured using the dynamic stress intensity factor as a function of the crack tip velocity. The simulated data are fitted using eq. (33).</figcaption></figure>

In this section, we calculate $K_{ID}$ for the 2D surface traction case (see results in section 4.1). This toughness value then will be correlated with the crack tip velocity. Our aim was to show how the mechanism leading to branching effects the measured macroscopic toughness. $g_c$ is kept constant independently of the deformation rate.

The macroscopic fracture toughness is determined using the dynamic J integral [87, 88] for tensile opening mode:

$$J_I = \int_S \left( \Pi_{\varepsilon} n_x - \mathbf{t} \frac{\partial \mathbf{u}}{\partial x} \right) dS + \int_A \rho \ddot{\mathbf{u}} \frac{\partial \mathbf{u}}{\partial x} dA, \quad (34)$$

where $\Pi_{\varepsilon}$ is the overall strain energy: $\Pi_{\varepsilon} = g(d)(\psi_0^+ + \psi_0^{pl}) + \psi_0^-$. The surface traction vector $\mathbf{t}$ can be calculated as: $\mathbf{t} = \sigma \cdot \mathbf{n}$, where $\mathbf{n}$ is the external normal unit vector to path $S$ which is a closed contour of area ($A$) around the crack tip. The second surface integral term takes the dynamic effect into account. We assume no body forces in our case.

The stress intensity factor is then calculated using the following assumption [89]:

$$K_I = \sqrt{\frac{8\mu J_I}{1 + \kappa}}, \quad (35)$$

where $\mu$ is the shear modulus and $\kappa = 3 - 4\nu$ for plane strain cases and $\nu$ is Poisson's ratio.

Fig. 13 shows the instantaneous stress intensity factor calculated between the initiation and the branching for different yield strength values. The dispersion observed around 700-800 m/s is attributed to the fluctuation of the crack velocity at $t = 30 - 40 \ \mu s$. A clear increment can be observed as a function of the propagation velocity. With solid lines eq. (33) is fitted on the simulated data. Furthermore, by decreasing the yield strength from brittle to 5 GPa, a gradual increase in macroscopic fracture toughness was found. This is in good qualitative correspondence with the static results as well as with experiments conducted on ship steel at different temperatures [2].

The observed increment in the fracture toughness is approximately 50%. This is smaller than the experimentally measured value [90]. However, Kalthoff [68] showed for Araldite B and Ravi-Chandar [91] for Homalite that based on the geometry this increment can vary between 30-110% for the same material.
<figure><img src="image_15.png" /><figcaption>Figure 14: Dynamic stress intensity factor as a function of normalized dissipation rate for the elastic/brittle case. To obtain the dissipation rate the fracture energy increment was normalized by $g_c$.</figcaption></figure>

To find a theoretical explanation for the increment in $K_{ID}$ without the change in $g_c$ we looked at the dissipation rate:

$$ \Gamma' = \frac{1}{h} \frac{d\Gamma}{da}, \qquad (36) $$

where $\Gamma$ is the fracture surface (see eq. (2)), $h$ is the thickness of the 2D elements and $a$ represents the length of the crack. This quantity basically describes the relative amount of fracture surface opening with a unit increment of the crack length. Fig. 14 shows $K_{ID}$ as a function of the dissipation rate. A linear relationship was found between the two quantities, which signifies that the increment of the instantaneous stress intensity factor is a result of the crack's widening.

This explanation might be considered as a particular case for the 2D phase-field model. However, several authors [50, 85] showed using 3D phase-field simulations, that when the local dissipation rate exceeds $\Gamma' \approx 2$ micro-branching occurs even before the macroscopic branching phenomenon. Unfortunately, the mentioned 3D nature of the crack propagation cannot be captured using 2D simulations [86]. Nevertheless, our observation shows that the phase-field approximation, independently from its limitations, can explain a very important physical phenomenon even in 2D.

The findings shown in Fig. 14 are supported by the experimental observations [92, 93, 94] who established a similar correlation between the roughness of the cracks surface, fracture velocity and the instantaneously measured stress intensity factor.

The increment of the process zone size can distort the displacement field and create a larger shadow in the caustic measurement [95], which can explain the generally higher $K_{ID}$ in macroscopic experiments [66]. However, depending on the rate sensitivity of the material (e.g. martensitic steel [5]) a rate dependent $g_c$ would be necessary to recreate the experimental measurement. Fortunately, the phase-field formulation is capable of handling the task [96].

It is a long lasting debate if the internal length scale ($l_c$) has a physical interpretation. Initially $l_c$ was only a numerical trick to solve the fracture mechanics problem with partial differential equation. However, recently it was suggested [42, 97] that $l_c$ can be correlated to the materials microstructure. Therefore, it is an important question that simply taking an $l_c$ close to zero is the real solution or by doing so we introduce another phenomenon which distorts the mechanical interpretation. Reducing $l_c$ to zero suggest that the material has no micro-structure and that it can be considered homogeneous
at any scale. As soon as, real experiments are analyzed, it is clear that there is a lower limit for $l_c$
bellow which the material has to be considered heterogeneous. The analysis, presented in Ref. [42, 97],
relying on a heterogeneous materials could maybe estimate this lower bound.

Furthermore, based on atomic scale simulations [98, 99] and nanoscopic experiments [100, 101] we
question the formulation of the plastic energy. Due to the small separation of the atomic bonds at
the micro-scale they can easily reform, leading to response where plasticity does not induce neither
damaged nor cracks. The behavior of the material remains plastic and a continuous flow stress
is observed. However, the void growth mechanism observed in experiments [102], that leads to
macroscopic coalescence and fracture propagation, indicates that locally the material is in tension.
Therefore, we argue that the degradation function in eq. (14) should be different or it should be
omitted based on the analyzed scale and the triaxiality [48] of the stress state. This would respond to
the problem of the anisotropic versus only hydrostatic energy degradation mechanism implemented
in the present case. Another option would be to use porous plastic models such as the Gurson-
Tvergaard-Needleman criterion [56, 103].

**6. Conclusion**

This study is set out to investigate the effect of plasticity on dynamic fracture propagation. An
elastoplastic version of the phase-field approach was implemented in the framework of a staggered
implicit dynamic time integration scheme. We chose to use the commercial finite element code Abaqus
with the UEL option because it is widely available in both industrial workspaces and in academic
laboratories. Therefore, this 2D and 3D implementation enables practicing engineers and scientists
to simulate easily not only static but dynamic and ductile crack propagation with the same UEL
subroutine.

With the help of this novel approach we studied the fracture patterns, the crack tip velocity and
the instantaneous stress intensity factor ($K_{ID}$) during dynamic fracture for various elastic and ductile
cases.

We have observed an abrupt shift between tensile and shear crack initiation as a function of the
shear strength and fracture toughness ratio. For static cases the presence of local ductile deformation
was clearly observable in the reaction force. With an optimal shear strength both resistance and
toughness can be increased. In dynamic crack propagation a similar phenomenon was observed.
Using an optimal value no branching was found, the crack propagated straight.

In both static and dynamic tests a high yield strength resulted in a perpendicular (mode I) tensile
fracture, while for low values the crack favored a tilted angle (mode II).

Without changing the materials fracture toughness ($g_c$) we observed an increment in the dynamic
stress intensity factor as a function of the crack velocity. This phenomenon is in agreement with
experimental results [66].

We explained this phenomenon by correlating the macroscopic fracture toughness to the relative
dissipation rate ($\Gamma'$). We found a linear correlation between the additional crack surface created in
a unit advancement and $K_{ID}$. This is also in agreement with past experimental measurements [93].
In case of a brittle calculation, when the dissipation rate exceeded $\Gamma' = 2.0$, the crack branched. The
mechanism of the local surface roughness was explained formerly by Bleyer [85], who showed that
in 3D conditions $\Gamma'$ can reach 2.0 locally creating a micro-branch before the macroscopic instability.
This way increasing the dissipated energy locally.

However, this is still not enough to reproduce experimental results conducted on materials sensitive to strain rate. Therefore, it would be advantageous to include both viscoelasticity, a rate
dependent $g_c$ or delayed damage [104] to further extend the analysis to these rate dependent materials. Furthermore, to fully understand the effect of the length scale it would be advantageous to compare phase-field simulations to other techniques such as molecular dynamics, eXtended FEM [69] or finite fracture mechanics [105]. The Abaqus implementation should be extended in the future to treat the phase-field problem with a simpler UMAT model using finite strains [106].

## Acknowledgment

The authors thank N. Moës for the interesting discussions and for sharing the details of the asymmetric energy decomposition algorithm. This work was supported by the French Research National Agency program through the grant ANR-16-CE30-0007-01.

## References

[1] A. A. Griffith, The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 221 (582-593) (1921) 163-198.

[2] M. F. Kanninen, C. H. Popelar, Advanced fracture mechanics, no. 15, Oxford University Press, 1985.

[3] H. Maigre, D. Rittel, Mixed-mode quantification for dynamic fracture initiation: application to the compact compression specimen, International Journal of Solids and Structures 30 (23) (1993) 3233-3244.

[4] M. Zhou, A. J. Rosakis, G. Ravichandran, Dynamically propagating shear bands in impact-loaded prenotched plates—I. experimental investigations of temperature signatures and propagation speed, Journal of the Mechanics and Physics of Solids 44 (6) (1996) 981-1006.

[5] J. F. Kalthoff, Modes of dynamic shear failure in solids, International Journal of Fracture 101 (1-2) (2000) 1-31.

[6] H. D. Bui, Mécanique de la rupture fragile, Masson, 1978.

[7] M. Adda-Bedia, R. Arias, M. B. Amar, F. Lund, Generalized Griffith criterion for dynamic fracture and the stability of crack motion at high velocities, Physical Review E 60 (2) (1999) 2366.

[8] L. B. Freund, Dynamic fracture mechanics, Cambridge university press, 1998.

[9] J. R. Rice, Some studies of crack dynamics, in: Physical Aspects of Fracture, Springer, 2001, pp. 3-11.

[10] S. J. Zhou, P. S. Lomdahl, R. Thomson, B. L. Holian, Dynamic crack processes via molecular dynamics, Physical review letters 76 (13) (1996) 2318.

[11] C. L. Rountree, R. K. Kalia, E. Lidorikis, A. Nakano, L. van Brutzel, P. Vashishta, Atomic aspects of crack propagation in brittle materials: Multimillion atom molecular dynamics simulations, Annual Review of Materials Research 32 (1) (2002) 377-400.
[12] Y.-L. Gui, H. H. Bui, J. Kodikara, Q.-B. Zhang, J. Zhao, T. Rabczuk, Modelling the dynamic failure of brittle rocks using a hybrid continuum-discrete element method with a mixed-mode cohesive fracture model, International Journal of Impact Engineering 87 (2016) 146–155.

[13] F. Zhou, J. F. Molinari, Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency, International Journal for Numerical Methods in Engineering 59 (1) (2004) 1–24.

[14] M. L. Falk, A. Needleman, J. R. Rice, A critical evaluation of cohesive zone models of dynamic fracture, Le Journal de Physique IV 11 (PR5) (2001) 43–50.

[15] A. E. Huespe, J. Oliver, P. J. Sanchez, S. Blanco, V. Sonzogni, Strong discontinuity approach in dynamic fracture simulations, Mecánica Computacional 25 (2006) 1997–2018.

[16] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1) (1999) 131–150.

[17] J. Réthoré, A. Gravouil, A. Combescure, An energy-conserving scheme for dynamic crack growth using the extended finite element method, International Journal for Numerical Methods in Engineering 63 (5) (2005) 631–659.

[18] B. Prabel, A. Combescure, A. Gravouil, S. Marie, Level set X-FEM non-matching meshes: application to dynamic crack propagation in elastic-plastic media, International Journal for Numerical Methods in Engineering 69 (8) (2007) 1553–1569.

[19] E. Sharon, S. P. Gross, J. Fineberg, Local crack branching as a mechanism for instability in dynamic fracture, Physical Review Letters 74 (25) (1995) 5096.

[20] T. Belytschko, H. Chen, J. Xu, G. Zi, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment, International journal for numerical methods in engineering 58 (12) (2003) 1873–1905.

[21] G. Pijaudier-Cabot, Z. P. Bažant, Nonlocal damage theory, Journal of engineering mechanics 113 (10) (1987) 1512–1533.

[22] M. Jirasek, Nonlocal models for damage and fracture: comparison of approaches, International Journal of Solids and Structures 35 (31-32) (1998) 4133–4145.

[23] N. Moës, C. Stolz, P.-E. Bernard, N. Chevaugeon, A level set based model for damage growth: The thick level set approach, International Journal for Numerical Methods in Engineering 86 (3) (2011) 358–380.

[24] P.-E. Bernard, N. Moës, N. Chevaugeon, Damage growth modeling using the thick level set (TLS) approach: Efficient discretization for quasi-static loadings, Computer Methods in Applied Mechanics and Engineering 233 (2012) 11–27.

[25] G. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (8) (1998) 1319 – 1342.

[26] G. Francfort, J. J. Marigo, Vers une théorie énergétique de la rupture fragile, Comptes Rendus Mécanique 330 (4) (2002) 225–233.
[27] B. Bourdin, G. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (4) (2000) 797–826.

[28] J.-J. Marigo, Constitutive relations in plasticity, damage and fracture mechanics based on a work property, Nuclear Engineering and Design 114 (3) (1989) 249–272.

[29] T. Gerasimov, L. D. Lorenzis, On penalization in variational phase-field models of brittle fracture, Computer Methods in Applied Mechanics and Engineering 354 (2019) 990–1026.

[30] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (45–48) (2010b) 2765 – 2778.

[31] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, International Journal for Numerical Methods in Engineering 83 (10) (2010a) 1273–1311.

[32] T. Heister, M. F. Wheeler, T. Wick, A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach, Computer Methods in Applied Mechanics and Engineering 290 (2015) 466–495.

[33] L. Kaczmarczyk, Z. Ullah, K. Lewandowski, X. Meng, X.-Y. Zhou, I. Athanasiadis, H. Nguyen, C.-A. Chalons-Mouriesse, E. Richardson, E. Miur, MoFEM: An open source, parallel finite element library, Journal of Open Source Software 5 (45) (2020) 1441.

[34] T. Wick, Modified newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation, Computer Methods in Applied Mechanics and Engineering 325 (2017) 577–611.

[35] J.-Y. Wu, Y. Huang, Comprehensive implementations of phase-field damage models in Abaqus, Theoretical and Applied Fracture Mechanics 106 (2020) 102440.

[36] F. Aldakheel, B. Hudobivnik, A. Hussein, P. Wriggers, Phase-field modeling of brittle fracture using an efficient virtual element scheme, Computer Methods in Applied Mechanics and Engineering 341 (2018) 443–466.

[37] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids 57 (8) (2009) 1209–1229.

[38] P. Sicsic, J.-J. Marigo, C. Maurini, Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling, Journal of the Mechanics and Physics of Solids 63 (2014) 256–284.

[39] P. Farrell, C. Maurini, Linear and nonlinear solvers for variational phase-field models of brittle fracture, International Journal for Numerical Methods in Engineering 109 (5) (2017) 648–667.

[40] J. L. Chaboche, F. Feyel, Y. Monerie, Interface debonding models: a viscous regularization with a limited rate dependency, International journal of solids and structures 38 (18) (2001) 3127–3160.
[41] E. Lorentz, A. Benallal, Gradient constitutive relations: numerical aspects and application to gradient damage, Computer methods in applied mechanics and engineering 194 (50-52) (2005) 5191–5220.

[42] T. T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, C. Chateau, A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure, Engineering Fracture Mechanics 139 (2015) 18 – 39.

[43] T. T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, C. Chateau, A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography, Computer Methods in Applied Mechanics and Engineering 312 (2016) 567–595.

[44] J.-J. Marigo, C. Maurini, K. Pham, An overview of the modelling of fracture by gradient damage models, Meccanica 51 (12) (2016) 3107–3128.

[45] B. Bourdin, G. A. Francfort, J.-J. Marigo, The Variational Approach to Fracture, Springer Netherlands, 2008.

[46] A. Karma, A. E. Lobkovsky, Unsteady crack motion and branching in a phase-field model of brittle fracture, Physical Review Letters 92 (24) (2004) 245510.

[47] B. Bourdin, C. J. Larsen, C. L. Richardson, A time-discrete model for dynamic fracture based on crack regularization, International journal of fracture 168 (2) (2011) 133–143.

[48] M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. Hughes, C. M. Landis, A phase-field description of dynamic brittle fracture, Computer Methods in Applied Mechanics and Engineering 217-220 (2012) 77 – 95.

[49] M. Hofacker, C. Miehe, Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation, International Journal of Fracture 178 (1) (2012) 113–129.

[50] H. Henry, M. Adda-Bedia, Fractographic aspects of crack branching instability using a phase-field model, Physical Review E 88 (6) (2013) 060401.

[51] T. Y. Li, J. J. Marigo, D. Guilbaud, S. Potapov, Variational approach to dynamic brittle fracture via gradient damage models, in: Applied mechanics and materials, Vol. 784, 2015, pp. 334–341.

[52] J. Carlsson, P. Isaksson, Dynamic crack propagation in wood fibre composites analysed by high speed photography and a dynamic phase field model, International Journal of Solids and Structures 144 (2018) 78–85.

[53] R. Alessi, J. J. Marigo, S. Vidoli, Gradient damage models coupled with plasticity: variational formulation and main properties, Mechanics of Materials 80 (2015) 351–367.

[54] M. Ambati, T. Gerasimov, L. D. Lorenzis, Phase-field modeling of ductile fracture, Computational Mechanics 55 (5) (2015) 1017–1040.
[55] C. Miehe, F. Aldakheel, A. Raina, Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory, International Journal of Plasticity 84 (2016) 1–32.

[56] Y. Zhang, E. Lorentz, J. Besson, Ductile damage modelling with locking-free regularised gtn model, International Journal for Numerical Methods in Engineering 113 (13) (2018) 1871–1903.

[57] R. Alessi, J. J. Marigo, C. Maurini, S. Vidoli, Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: one-dimensional examples, International Journal of Mechanical Sciences 149 (2018) 559–576.

[58] M. Dittmann, F. Aldakheel, J. Schulte, P. Wriggers, C. Hesch, Variational phase-field formulation of non-linear ductile fracture, Computer Methods in Applied Mechanics and Engineering 342 (2018) 71–94.

[59] J. Fang, C. Wu, J. Li, Q. Liu, C. Wu, G. Sun, L. Qing, Phase field fracture in elasto-plastic solids: variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening, International Journal of Mechanical Sciences 156 (2019) 382–396.

[60] J. Fang, C. Wu, T. Rabczuk, C. Wu, C. Ma, G. Sun, Q. Li, Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies, Theoretical and Applied Fracture Mechanics 103 (2019) 102252.

[61] C. McAuliffe, H. Waisman, A coupled phase field shear band model for ductile–brittle transition in notched plate impacts, Computer Methods in Applied Mechanics and Engineering 305 (2016) 173–195.

[62] M. J. Borden, T. J. Hughes, C. M. Landis, A. Anvari, I. J. Lee, A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects, Computer Methods in Applied Mechanics and Engineering 312 (2016) 130 – 166.

[63] G. Liu, Q. Li, M. A. Msekh, Z. Zuo, Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model, Computational Materials Science 121 (2016) 35 – 47.

[64] E. Azinpour, J. P. S. Ferreira, M. P. L. Parente, J. C. de Sa, A simple and unified implementation of phase field and gradient damage models, Advanced Modeling and Simulation in Engineering Sciences 5 (1) (2018) 1–24.

[65] E. Martínez-Pañed, A. Golahmar, C. F. Niordson, A phase field formulation for hydrogen assisted cracking, Computer Methods in Applied Mechanics and Engineering 342 (2018) 742–761.

[66] J. W. Dally, W. L. Fourney, G. R. Irwin, On the uniqueness of the stress intensity factor — crack velocity relationship, International Journal of Fracture 27 (3) (1985) 159–168.

[67] D. Broek, Elementary engineering fracture mechanics, Springer Science & Business Media, 1982.

[68] J. F. Kalthoff, On some current problems in experimental fracture dynamics, Tech. rep., Fraunhofer Institute for Mechanics of Materials, Freiburg, Germany (1983).
[69] D. Grégoire, H. Maigre, J. Réthoré, A. Combescure, Dynamic crack propagation under mixed-mode loading-comparison between experiments and X-FEM simulations, International journal of solids and structures 44 (20) (2007) 6517-6534.

[70] G. Molnár, A. Gravouil, 2D and 3D abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture, Finite Elements in Analysis and Design 130 (2017) 27 – 38.

[71] R. von Mises, Mechanik der festen körper im plastisch-deformablen zustand, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1913 (1913) 582-592.

[72] F. Freddi, G. Royer-Carfagni, Regularized variational theories of fracture: a unified approach, Journal of the Mechanics and Physics of Solids 58 (8) (2010) 1154–1174.

[73] R. Alessi, F. Freddi, L. Mingazzi, Phase-field numerical strategies for deviatoric driven fractures, Computer Methods in Applied Mechanics and Engineering 359 (2020) 112651.

[74] A. D. Al, J. Bruchon, S. Drapier, S. Fayolle, Integrating a logarithmic-strain based hyperelastic formulation into a three-field mixed finite element formulation to deal with incompressibility in finite-strain elastoplasticity, Finite Elements in Analysis and Design 86 (2014) 61–70.

[75] F. Aldakheel, C. Miehe, Coupled thermomechanical response of gradient plasticity, International Journal of Plasticity 91 (2017) 1–24.

[76] E. Oñate, R. Owen, Computational Plasticity, Springer Netherlands, 2007.

[77] K. Krabbenhøft, Basic computational plasticity, Technical University of Denmark, 2002.

[78] R. Hill, A theory of the yielding and plastic flow of anisotropic metals, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 193 (1033) (1948) 281-297.

[79] M. Frémond, B. Nedjar, Damage, gradient of damage and principle of virtual power, International Journal of Solids and Structures 33 (8) (1996) 1083 – 1103.

[80] H. M. Hilber, T. J. Hughes, R. L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engineering & Structural Dynamics 5 (3) (1977) 283-292.

[81] C. Miehe, M. Hofacker, L.-M. Schänzel, F. Aldakheel, Phase field modeling of fracture in multi-physics problems. part II. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids, Computer Methods in Applied Mechanics and Engineering 294 (2015) 486 – 522.

[82] B. Dodd, Y. Bai, Adiabatic shear localization: frontiers and advances, Elsevier, 2012.

[83] S. Ham, K.-J. Bathe, A finite element method enriched for wave propagation problems, Computers & Structures 94 (2012) 1–12.

[84] H. Gao, A theory of local limiting speed in dynamic fracture, Journal of the Mechanics and Physics of Solids 44 (9) (1996) 1453–1474.
[85] J. Bleyer, J.-F. Molinari, Microbranching instability in phase-field modelling of dynamic brittle fracture, Applied Physics Letters 110 (15) (2017) 151903.

[86] J. Bleyer, C. Roux-Langlois, J. F. Molinari, Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms, International Journal of Fracture 204 (1) (2017) 79–100.

[87] K. Kishimoto, S. Aoki, M. Sakata, On the path independent integral–J, Engineering Fracture Mechanics 13 (4) (1980a) 841–850.

[88] K. Kishimoto, S. Aoki, M. Sakata, Dynamic stress intensity factors using J–integral and finite element method, Engineering Fracture Mechanics 13 (2) (1980b) 387–394.

[89] D. R. J. Owen, A. J. Fawkes, Engineering fracture mechanics: numerical methods and applications, Pineridge Press Ltd, 91 West Cross Lane, West Cross, Swansea, UK, 1983. 305.

[90] A. J. Rosakis, A. T. Zehnder, On the dynamic fracture of structural metals, International Journal of Fracture 27 (3) (1985) 169–186.

[91] K. Ravi-Chandar, W. G. Knauss, Processes controlling the dynamic fracture of brittle solids, Tech. rep., California Institute of Technology, Pasadena, USA (1983).

[92] E. Sharon, S. P. Gross, J. Fineberg, Energy dissipation in dynamic fracture, Physical Review Letters 76 (1996) 2117–2120.

[93] K. Takahashi, M. Kido, K. Arakawa, Fracture roughness evolution during mode I dynamic crack propagation in brittle materials, International journal of fracture 90 (1-2) (1998) 119–131.

[94] L. V. Zhao, D. Bardel, A. Maynadier, D. Nelias, Velocity correlated crack front and surface marks in single crystalline silicon, Nature communications 9 (1) (2018) 1298.

[95] C. Spyropoulos, Stress intensity factor determination error by the method of caustics, Theoretical and applied fracture mechanics 35 (2) (2001) 179–186.

[96] R. Shen, H. Waisman, L. Guo, Fracture of viscoelastic solids modeled with a modified phase field method, Computer Methods in Applied Mechanics and Engineering 346 (2019) 862–890.

[97] P. Chakraborty, Y. Zhang, M. R. Tonks, Multi-scale modeling of microstructure dependent intergranular brittle fracture using a quantitative phase-field based method, Computational Materials Science 113 (2016) 38 – 52.

[98] G. Molnár, P. Ganster, A. Tanguy, E. Barthel, G. Kermouche, Densification dependent yield criteria for sodium silicate glasses – an atomistic simulation approach, Acta Materialia 111 (2016) 129 – 137.

[99] G. Molnár, P. Ganster, A. Tanguy, Effect of composition and pressure on the shear strength of sodium silicate glasses: An atomic scale simulation study, Physical Review E 95 (2017) 043001.

[100] G. Kermouche, G. Guillonneau, J. Michler, J. Teisseire, E. Barthel, Perfectly plastic flow in silica glass, Acta Materialia 114 (2016) 146 – 153.
[101] N. L. Okamoto, M. Inomoto, H. Adachi, H. Takebayashi, H. Inui, Micropillar compression deformation of single crystals of the intermetallic compound ζ-FeZn13, Acta Materialia 65 (2014) 229-239.

[102] H. Andersson, Analysis of a model for void growth and coalescence ahead of a moving crack tip, Journal of the Mechanics and Physics of Solids 25 (3) (1977) 217-233.

[103] M. Dittmann, F. Aldakheel, J. Schulte, F. Schmidt, M. Krüger, P. Wriggers, C. Hesch, Phase-field modeling of porous-ductile fracture in non-linear thermo-elasto-plastic solids, Computer Methods in Applied Mechanics and Engineering 361 (2020) 112730.

[104] O. Allix, J.-F. Deü, Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading, Engineering transactions 45 (1) (1997) 29-46.

[105] J. Li, D. Leguillon, E. Martin, X.-B. Zhang, Numerical implementation of the coupled criterion for damaged materials, International Journal of Solids and Structures 165 (2019) 93-103.

[106] C. Hesch, A. J. Gil, R. Ortigosa, M. Dittmann, C. Bilgen, P. Betsch, M. Franke, A. Janz, K. Weinberg, A framework for polyconvex large strain phase-field methods to fracture, Computer Methods in Applied Mechanics and Engineering 317 (2017) 649-683.

[107] C. Kuhn, T. Noll, R. Müller, On phase field modeling of ductile fracture, GAMM-Mitteilungen 39 (1) (2016) 35-54.

[108] P. Rodriguez, J. Ulloa, C. Samaniego, E. Samaniego, A variational approach to the phase field modeling of brittle and ductile fracture, International Journal of Mechanical Sciences 144 (2018) 502-517.

[109] H. Ulmer, M. Hofacker, C. Miehe, Phase field modeling of brittle and ductile fracture, PAMM 13 (1) (2013) 533-536.

[110] M. Ambati, R. Kruse, L. D. Lorenzis, A phase-field model for ductile fracture at finite strains and its experimental verification, Computational Mechanics 57 (1) (2016) 149-167.
# Appendix A. Asymmetric energy decomposition

Tensile and compression energies were decomposed by using the method developed in Ref. [24]. For practical reasons the calculation of the stress tensor and the tangent matrix is shown using Einstein notation.

The potential energy density of an elastic body is written in eq. (7), which equation is practically expressed as:

$$ \psi^{el}(\mathbf{u}, d) = \sum_{i=1}^{3} g(d \cdot \alpha_i) \mu \varepsilon_i^2 + \frac{\lambda}{2} g(d \cdot \alpha) \mathrm{tr}(\boldsymbol{\varepsilon})^2, \quad (A.1) $$

where constants $\alpha$ and $\alpha_i$ control whether the potential energy and thus the stiffness is degraded:

$$ \begin{aligned} \alpha_i = 0 & \quad \text{if} \quad \varepsilon_i < 0, \\ 1 & \quad \text{if} \quad \varepsilon_i \ge 0, \\ \alpha = 0 & \quad \text{if} \quad \mathrm{tr}(\boldsymbol{\varepsilon}) < 0, \\ 1 & \quad \text{if} \quad \mathrm{tr}(\boldsymbol{\varepsilon}) \ge 0. \end{aligned} \quad (A.2) $$

Thus, if the elastic principal strain value ($\varepsilon_i$) is positive (in tension) the related potential energy is degraded by function $g$ (see eq. (8)), otherwise the damage has no effect.

From the potential energy the elastic stress is obtained by:

$$ \sigma_{ij} = \frac{\partial \psi^{el}}{\partial \varepsilon_{ij}} = \frac{\partial \psi^{el}}{\partial \varepsilon_k} \frac{\partial \varepsilon_k}{\partial \varepsilon_{ij}} = L_{kn} \varepsilon_n \frac{\partial \varepsilon_k}{\partial \varepsilon_{ij}}, \quad (A.3) $$

where $\varepsilon_n$ are the eigenvalues and $\varepsilon_{ij}$ are the original components of the elastic strain tensor. The elementary stiffness of the material is calculated as follows:

$$ L_{kn} = \frac{\partial^2 \psi^{el}}{\partial \varepsilon_k \partial \varepsilon_n}. \quad (A.4) $$

The eigenvalues are the roots of the characteristic polynomial expressed with the matrix invariants:

$$ \hat{\varepsilon}^3 - i_1 \hat{\varepsilon}^2 + i_2 \hat{\varepsilon} - i_3 = 0. \quad (A.5) $$

This equation was solved with Cardano's method. The derivatives in eq. (A.3) can be expressed as follows:

$$ \frac{\partial \varepsilon_k}{\partial \varepsilon_{ij}} = \frac{\partial \varepsilon_k}{\partial i_n} \frac{\partial i_n}{\partial \varepsilon_{ij}}. \quad (A.6) $$

The derivatives of the eigenvalues with respect to the tensor invariants are obtained through implicit differentiation, e.g.:

$$ \frac{\partial \varepsilon_k}{\partial i_1} = \frac{\varepsilon_k^2}{3\varepsilon_k^2 - 2i_1\varepsilon_k + i_2}. \quad (A.7) $$

The elastic tangent matrix can be calculated by applying further derivation:

$$ H_{ijkl} = \frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}} = \frac{\partial}{\partial \varepsilon_{kl}} \left[ L_{mn} \varepsilon_m \frac{\partial \varepsilon_n}{\partial \varepsilon_{ij}} \right] = \frac{\partial \varepsilon_n}{\partial \varepsilon_{kl}} L_{mn} \frac{\partial \varepsilon_m}{\partial \varepsilon_{ij}} + L_{mn} \varepsilon_m \frac{\partial^2 \varepsilon_n}{\partial \varepsilon_{ij} \partial \varepsilon_{kl}}. \quad (A.8) $$

The second order derivatives of the elastic principal strains are calculated similarly to the first order components:
<figure><img src="image_16.png" /><figcaption>Figure A.15: Relative difference in the decomposed and theoretical stiffness matrix as a function of the discriminant of polynomial equation (A.5).</figcaption></figure>

$$ \frac{\partial^2 \varepsilon_n}{\partial \varepsilon_{ij} \partial \varepsilon_{kl}} = \frac{\partial}{\partial \varepsilon_{kl}} \left[ \frac{\partial \varepsilon_n}{\partial i_m} \frac{\partial i_m}{\partial \varepsilon_{ij}} \right] = \frac{\partial \varepsilon_n}{\partial i_m} \frac{\partial^2 i_m}{\partial \varepsilon_{ij} \partial \varepsilon_{kl}} + \frac{\partial^2 \varepsilon_n}{\partial i_m \partial i_p} \frac{\partial i_p}{\partial \varepsilon_{kl}} \frac{\partial i_m}{\partial \varepsilon_{ij}}. \quad (A.9) $$

Following the manipulations in eq (A.7), one example of a second order derivative is:

$$ \frac{\partial^2 \varepsilon_k}{\partial i_1^2} = \frac{2\varepsilon_k \partial \varepsilon_k / \partial i_1}{3\varepsilon_k^2 - 2i_1\varepsilon_k + i_2} - \frac{\varepsilon_k^2 (6\varepsilon_k \partial \varepsilon_k / \partial i_1 - 2i_1 \partial \varepsilon_k / \partial i_1 - 2\varepsilon_k)}{(3\varepsilon_k^2 - 2i_1\varepsilon_k + i_2)^2}. \quad (A.10) $$

We note that due to the symmetry in the deformation tensor, Voigt notation was used. Thus, 3 by 3 matrices were substituted with a 6 element vector. This way, for example, eq. (A.6) resulted in a 3 by 6, or eq. (A.9) in a (3 × 6 × 6) array.

**Numerical issues.** As it was reported previously by [24] the asymmetric energy degradation introduces a non-smooth transition in the stiffness matrix during iteration. Therefore, in present implementation the update of the stiffness matrix **H** is limited to the first four Newton-Raphson cycles. To further simplify the calculation, if all eigenvalues had the same sign or *d* = 0, no decomposition was applied.

The above algorithm becomes unreliable in case of double roots when the discriminant (Δ) of eq. (A.5) approaches zero. In this case the denominator of eq. (A.7) becomes also zero and the results become incorrect. To demonstrate this, a simple experiment is carried out. Let us consider the following strain state: ε<sub>x</sub> = 1 + δ, ε<sub>y</sub> = -ε<sub>z</sub> = 1 and γ<sub>xy</sub> = γ<sub>xz</sub> = γ<sub>yz</sub> = 0. δ is a small perturbation, which can be used to increase the discriminant, if δ = 0 → Δ = 0. Assuming that *d* = 0 the stiffness matrix was calculated and compared to the theoretical one (non-decomposed, **C**): η = max (|**H** - **C**|) / max (**C**).

Fig. A.15 shows the relative error as a function of the calculated discriminant. A power law relationship can be observed which decreases significantly if the discriminant increases.

As it was suggested by [24], when Δ = 0 a small perturbation should be added to ε. To analyze the effect of this modification a matrix with three well separated eigenvalues will be studied. During the calculation the second eigenvalue is going to be modified by δ, and the resulting stiffness matrices will be compared to the undisturbed one similarly as it was done previously. Fig. A.16 shows the difference introduced by the perturbation as a function of damage. It can be seen, that when the material is undamaged the perturbation causes no error, while at *d* = 1 the difference saturates almost
<figure><img src="image_17.png" /><figcaption>Figure A.16: Relative difference in the stiffness matrices introduced by the perturbation as a function of the damage variable. Without correction (solid symbols) and with the recalculation of the strain tensor shown in eq. (A.11).</figcaption></figure>

to $\delta$. While, this value can be reduced to $\delta/10$ by recalculating the strain tensor based on the new eigenvalues and the original eigenvectors:

$$ \varepsilon' = V\hat{\varepsilon}'V^T \qquad (A.11) $$

where $\varepsilon'$ is the modified strain tensor, $\hat{\varepsilon}$ is a diagonal matrix containing the eigenvalues (with $\hat{\varepsilon}'_2 = (1 + \delta)\hat{\varepsilon}_2$) and $V$ is a matrix containing the eigenvectors of $\varepsilon$.

As a result of these two tests, in the published implementation a small $\delta = 10\%$ perturbation is applied on the second eigenvalue if $\Delta \le 10^{-7}$.

## Appendix B. Details of the implementation in Abaqus/UEL

To implement the solution in Abaqus two element types are used in a layered manner. Each layer connects to the same nodes, but contributes to the stiffness of different degrees of freedom (DOF). A schematic illustration is depicted in Fig. B.17. The first element type solves the displacement problems, hence it has 2 or 3 translational DOFs. While the second element determines the fracture topology and has only one DOF, the phase-field.

In 2D both four-node isoparametric and three-node triangular elements are implemented. Similarly in 3D we considered eight-node brick and the four-node tetrahedral formulations.

In order to visualize the calculated quantities in Abaqus a third layer is added with infinitesimally small stiffness made from a UMAT (user defined material model). It is used to transfer information from the common block and interpolate between integration points. The internal variables are summarized in Table B.1 for both 2D plane strain and 3D elements.

Additionally, as global history output the kinetic energy (ALLKE), the elastic strain energy (ALLIE), the plastic strain energy (ALLPD) and the fracture energy (ALLEE) is available.

The UEL option in Abaqus requires the definition of the elementary stiffness matrix and the residue vector. In eq. (27) we gave the elementary linear equation system which is used to solve the problem iteratively. The d'Alambert force term can be calculated as:

$$ \mathbf{f}^{\text{ine}} = f_i^{\text{ine}} = \int_{\Omega} \rho N_{ij}^u \ddot{u}_i d\Omega, \qquad (B.1) $$
<figure><img src="image_18.png" /><figcaption>Figure B.17: 2D representation of three layered finite element structure in Abaqus. All nodes have three (or in 3D four) degrees of freedom (DOF). The first element contributes to the stiffness of the displacement DOFs, the second element to the damage DOF. For post processing purposes a third layer is included made as a UMAT model, which allows to display state dependent variables (SDVs).</figcaption></figure>

<table><caption>Table B.1: Solution dependent variables (SDV) used to plot the results in two and three dimensions.</caption><thead><tr><th rowspan="2">Variable</th><th colspan="2">Number of SDV in Abaqus</th></tr><tr><th>2D (x, y)</th><th>3D (x, y, z)</th></tr></thead><tbody><tr><td>displacement - u<sub>x</sub>, u<sub>y</sub>, u<sub>z</sub></td><td>SDV1-SDV2</td><td>SDV1-SDV3</td></tr><tr><td>axial strains - ε<sub>x</sub>, ε<sub>y</sub>, ε<sub>z</sub></td><td>SDV3-SDV4</td><td>SDV4-SDV6</td></tr><tr><td>engineering shear strains - γ<sub>xy</sub>, γ<sub>xz</sub>, γ<sub>yz</sub></td><td>SDV5</td><td>SDV7-SDV9</td></tr><tr><td>elastic strains - ε<sup>el</sup></td><td>SDV6-SDV9</td><td>SDV10-SDV15</td></tr><tr><td>plastic strains - ε<sup>pl</sup></td><td>SDV10-SDV13</td><td>SDV16-SDV21</td></tr><tr><td>equivalent plastic strain - ε<sup>pl</sup><sub>eq</sub></td><td>SDV14</td><td>SDV22</td></tr><tr><td>stresses - σ</td><td>SDV15-SDV18</td><td>SDV23-SDV28</td></tr><tr><td>hydrostatic stress - tr(σ)/3</td><td>SDV19</td><td>SDV29</td></tr><tr><td>von Mises stress - σ<sub>eq</sub></td><td>SDV20</td><td>SDV30</td></tr><tr><td>plastic energy - ψ<sup>pl</sup><sub>0</sub></td><td>SDV21</td><td>SDV31</td></tr><tr><td>elastic tensile energy - ψ<sup>+</sup><sub>0</sub></td><td>SDV22</td><td>SDV32</td></tr><tr><td>potential energy - Π<sub>ε</sub></td><td>SDV23</td><td>SDV33</td></tr><tr><td>phase-field - d</td><td>SDV24</td><td>SDV34</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>Figure C.18: Equivalent stress-strain curve for a perfectly plastic case with different $I_\epsilon$ values. With different colors, the different deformation states are highlighted: elastic (green), plastic (blue), critical (orange) and fracturing (red).</figcaption></figure>

where $\rho$ is the density, $\ddot{u}_i$ are the nodal accelerations and $N_{ij}^u$ is the matrix containing the shape functions for the displacement problem. The "lumped" (diagonal) mass matrix is:

$$ \mathbf{M} = M_{kk} = \int_{\Omega} \rho N_{ik}^{u} N_{il}^{u} d\Omega. \qquad (B.2) $$

While, the definitions of variables: $\mathbf{r}^d$, $\mathbf{K}^d$, $\mathbf{f}^{\text{ext}}$, $\mathbf{f}_n^{\text{int}}$ and $\mathbf{K}_n^u$ can be found in Ref. [70].

## Appendix C. Critical plastic strain

In the literature two major groups of ductile phase-field implementations exist. The first group [107, 62, 108, 59] uses the energy of the plastic deformation to sustain damage evolution, while the other [109, 54, 110] defines a critical strain and a penalty function to achieve fracture. The advantage of the first method is that the strain energy function satisfies the laws of thermodynamics, however in plasticity if the material healing it can maintain part of its resistance independently of the plastic strain [98]. On the other hand it is much easier to measure a macroscopic strain at failure, than the energy necessary to open a unit crack surface in ductile cases.

Therefore, in the present implementation both methods are included. As it is shown in eq. (15) the plastic energy is degraded by the damage function, thus in the energy history both elastic and plastic energies are included. Nevertheless, a penalty function is added which induces fracture if a critical plastic strain is achieved. This section is dedicated to the definition of $I_\epsilon$.

The stress-strain curve of an elastic, perfectly plastic solid is shown in Fig. C.18. With different colors we highlighted the elastic ($\psi_{el}^{max}$), plastic ($\psi_{pl}^{max}$), penalty ($\psi_{cr}^{max}$) and fracturing energy components:

$$ \psi_{el}^{\max} = \frac{1}{6G} \left[ \sigma^{\lim} + \frac{1}{2} H (\varepsilon_{eq}^{cr} + \Delta\varepsilon_{eq}^{cr}) \right]^2, \qquad (C.1) $$

$$ \psi_{pl}^{\max} = \sigma^{\lim} (\varepsilon_{eq}^{cr} + \Delta\varepsilon_{eq}^{cr}) + \frac{1}{2} H (\varepsilon_{eq}^{cr} + \Delta\varepsilon_{eq}^{cr})^2, \qquad (C.2) $$

$$ \psi_{cr}^{\max} = \frac{1}{2} I_{\epsilon} \langle \varepsilon_{eq}^{pl} - \varepsilon_{eq}^{cr} \rangle^2 = \frac{1}{2} I_{\epsilon} (\Delta \varepsilon_{eq}^{cr})^2. \qquad (C.3) $$
It can be seen, that even with a large $I_ε$ there is an additional deformation component to $ε_{eq}^{cr}$. For simplicity $Δε_{eq}^{cr}$ is set to 10% of $ε_{eq}^{cr}$. Thus, $I_ε$ is calculated as follows:

$$I_{\epsilon} = \frac{2 \langle \psi_c - \psi_{el}^{\max} - \psi_{pl}^{\max} \rangle}{(\Delta \epsilon_{eq}^{cr})^2} \qquad (C.4)$$

This can, of course, be changed in the UEL script when needed.

Note, that if the user desires not to use the penalty function $ε_{eq}^{cr}$ must be set to 0.

## Appendix D. Recommendation for using the UEL

Abaqus proposes an alternative way, where the user can define its own element. This is called the UEL option. This section is dedicated to explain the built in options of the UEL subroutine included in supplementary materials. With the paper a 2D plane strain and a 3D example is published.

The basic construction of the model should look like a layered finite element; where the first layer is made of elements (U1 (triangular) or U2 (square)) with translational degrees of freedom. This should be followed by layer with phase-field elements (U3 (triangular) and U4 (square)) defined on the same nodes. Finally a third layer of UMAT is recommended, which simplifies post-processing.

For the elements solving the displacement problem there are several options to consider. The element properties should be added in the following order: Young's modulus (E), Poisson's ratio (ν), yield strength ($σ^{lim}$), hardening modulus (H), critical strain ($ε_{eq}^{cr}$), thickness (only for 2D), density ($ρ$), asymmetric energy switch, plasticity switch, fracture surface energy ($g_c$), internal length scale ($l_c$). For every switch variable, the value 1 should be set if the option is used, while 0 is the it should be omitted.

* **Asymmetric energy switch:** If the option is used, the energy is decomposed into tensile and compression parts. However, if the plastic option is applied parallel, only the hydrostatic components are decomposed.

* **Plastic switch:** If it is used, a von Mises yield criterion is applied.

* **Critical strain:** If the penalty function needs to be turned off, the critical strain should be set to zero.

The fracture properties ($l_c$ and $g_c$) are added to this element as well to control the history energy increment and the time step.

For the phase-field element, the properties are: fracture surface energy ($g_c$), internal length scale ($l_c$), thickness (only for 2D), elastic threshold switch. If the elastic threshold is used $ψ_c = g_c/2l_c$, else it is zero.

The element chooses automatically, if static or dynamic step integration is requested.

The implemented UEL Fortran code is governed by the CeCILL-B license. The code is chosen to be published open-source to facilitate the distribution and the ease of use of the phase-field technique among Abaqus users. While on the personal website of the corresponding author (address shown in the Fortran file), regular updates and conversion scripts are available through elaborate tutorials.

## Appendix E. Supplementary data

Supplementary data associated with this article can be found in the online version at URL.