Assessment of four strain energy decomposition methods for phase field fracture models using quasi-static and dynamic benchmark cases

Strain energy decomposition methods in phase field fracture models separate strain energy that contributes to fracture from that which does not. However, various decomposition methods have been proposed in the literature, and it can be difficult to determine an appropriate method for a given problem. The goal of this work is to facilitate the choice of strain decomposition method by assessing the performance of three existing methods (spectral decomposition of the stress or the strain and deviatoric decomposition of the strain) and one new method (deviatoric decomposition of the stress) with several benchmark problems. In each benchmark problem, we compare the performance of the four methods using both qualitative and quantitative metrics. In the first benchmark, we compare the predicted mechanical behavior of cracked material. We then use four quasi-static benchmark cases: a single edge notched tension test, a single edge notched shear test, a three-point bending test, and a L-shaped panel test. Finally, we use two dynamic benchmark cases: a dynamic tensile fracture test and a dynamic shear fracture test. All four methods perform well in tension, the two spectral methods perform better in compression and with mixed mode (though the stress spectral method performs the best), and all the methods show minor issues in at least one of the shear cases. In general, whether the strain or the stress is decomposed does not have a significant impact on the predicted behavior.


Introduction
Phase field fracture models were first developed in the late 1990s, and have long been established as an efficient method for modeling brittle fracture. Their popularity is due to their remarkable ability to accurately represent various complex fracture phenomena in both 2D and 3D, including crack initiation, propagation, branching, and merging. Phase field models use continuous variable fields to represent material structure, and define the evolution of these fields, and thus the evolution of the structure, by minimizing the total free energy of the system (Chen 2002;Moelans et al. 2008;Steinbach 2013;Tonks and Aagesen 2019). In phase field fracture models, a continuous variable field is used to represent cracks and it is evolved to minimize the total strain energy of the system. Due to the continuous nature of the crack field, cracks have a finite width and the need to track the crack surface is eliminated.
The first phase field fracture models were developed between the late 1990s (Francfort and Marigo 1998; Bourdin et al. 2000;Aranson et al. 2000) and early 2000s (Karma et al. 2001;Henry and Levine 2004;Kuhn and Müller 2008;Hakim and Karma 2009). Later, Borden et al. (2014) built a higher-order element model based on the isogeometric analysis framework proposed by Hughes et al. (2005). Ambati et al. (2015a) proposed a hybrid phase field fracture model that enables a significant reduction in computational cost but is not thermodynamically consistent. Wu (2017) built a unified phase field fracture model for static and quasi-static fracture problems by including multiple degradation functions. Phase field fracture models have been applied in many areas, including cohesive zone modeling (Verhoosel and de Borst 2013;May et al. 2015;Nguyen and Wu 2018), ductile fracture (Ambati et al. 2015b;Ambati et al. 2016;Dittmann et al. 2018;Azinpour et al. 2021), hydraulic fracture (Miehe and Mauthe 2016;Wilson and Landis 2016;Heider and Markert 2017;Zhuang et al. 2020;Heider et al. 2018), and thermomechanical fracture (Kuhn and Müller 2009;Nguyen et al. 2019;Svolos et al. 2020;Rezwan et al. 2021). Phase field fracture methods have also been expanded to model dynamic fracture (Karma et al. 2001;Bourdin et al. 2011;Larsen 2010), with the most widely used approach developed by Borden et al. (2012). Phase field fracture models have been implemented in various finite element software tools, including ABAQUS (Msekh et al. 2015), COMSOL (Zhou et al. 2018), and the open source MOOSE framework (Chakraborty et al. 2016;. In all phase field fracture models, a degradation function is used to separate strain energy that can be stored by cracked material from energy that cannot be. The energy that cannot be stored directly contributes to the formation and propagation of cracks and the energy that can be stored does not. Early models assumed that cracked material could not store any strain energy. Later, methods were developed to decompose the strain energy in a manner representative of the behavior of actual cracked material. The method by Amor et al. (2009) decomposed the deviatoric and volumetric parts of the strain tensor, and only the deviatoric and positive volumetric parts contributed to fracture. The method by 2010) used a spectral decomposition method in which only the part associated with positive eigenvalues of the strain tensor contributed to fracture. Despite the differences in their approaches, both methods successfully capture observed fracture behavior. A newer decomposition method that is similar to that from Miehe et al. applied a spectral decomposition on the stress rather than the strain tensor and was applicable to anisotropic elastic responses . Another new method decomposed the stress tensor by separating stress components normal to the crack surface from those that are tangential .
The strain energy decomposition method is a critical part of a phase field fracture model. However, since there are multiple methods available in the literature, it is difficult to know which method is most applicable for a given problem. The goal of this work is to assess four decomposition methods, three from the literature and one that is proposed here for the first time. The first is the approach from Amor et al. (2009) in which the strain tensor is separated into volumetric and deviatoric parts. The second is the approach from 2010), which applies spectral decomposition to the strain tensor. The third is the approach from , which applies a spectral decomposition to the stress tensor. The fourth model, which is new to this work, separates the stress tensor into deviatoric and volumetric parts. We focus on the spectral and deviatoric decomposition methods because they are the most widely used; the new method that separates the normal and tangential stress components with respect to the crack surface ) warrants future study. Basic quasi-static and dynamic benchmark examples are solved with these four models and their results are compared. In "General equations of phase field fracture models" section, the general governing equations for phase field fracture models are summarized, along with the four decomposition methods. The mechanical behaviors of cracked materials are compared in "Mechanical response of cracked material" section. The performance of the four decomposition methods is compared for quasi-static fracture examples in "Quasi-static fracture examples" section and for dynamic fracture examples in "Dynamic fracture examples and results" section. "Discussion" section discusses our results and this work is concluded in "Conclusions" section.

Governing equations
In this section, we summarize the general governing equations for phase field fracture models. More details can be found in the literature (Amor et al. 2009;. In these models, a continuous variable field d ranging from 0 to 1 is used to describe the crack behavior, where d = 0 means the material is intact and d = 1 means it is fully cracked. Two governing equations are solved to define the crack evolution; the first controls the crack variable evolution and the second the displacement field.
The crack variable field d evolves to minimize the total free energy in the system. The total free energy in a damaged material consists of strain energy and fracture energy, according to where g c is the crack energy release rate, l is a length parameter that defines the width of the crack, and g(d) is the degradation function with k a small number to improve numerical convergence. The strain energy density is divided into the portion that contributes to fracture ψ + and the portion that does not ψ − . The d evolution equation must include an approach to prevent crack healing. In this work we use the approach from  that employs a history variable H that is equal to the largest ψ + value experienced at a given location up to the current time t c , i.e.
Gerasimov et al. (Gerasimov and De Lorenzis 2019) discuss the necessity of an irreversibility constraint to avoid crack healing and summarize alternative approaches to the one used in this work. The damage variable d evolves according to where L is the mobility parameter. We assume linear elastic material response, such that the strain tensor is determined from the displacement vector field u according to The stress in intact material σ 0 is obtained according to where C is the rank four elasticity tensor. In all the examples in this work we assume a homogeneous isotropic elastic constitutive response, such that C ijkl = λδ ij δ kl + μ(δ ik δ jl + δ il δ jk ) throughout the material, where λ and μ are the Lamé elastic constants and δ ij is the Kronecker-delta. In a cracked material, the local stress will depend on the local value of the crack field d according to where σ + is the portion of the stress tensor not present in cracked material and σ − is the portion that is present. Neglecting the body force, the displacement vector field u is determined using for quasi-static fracture and for dynamic fracture, where ρ is the density andü is the acceleration vector (second time derivative of the displacement vector u). Both come with the boundary conditions, where n is the out normal vector on the surface boundary, t is the surface traction, and u is the displacement boundary condition. For dynamic fracture, we assume no initial displacements.

Strain energy decomposition methods
For the phase field fracture model, a decomposition method is needed to determine the positive and negative parts of the strain energy density ψ + and ψ − used in Eq.
(3) to define the evolution of d. The decomposition method must also determine the positive and negative stress tensors σ + and σ − used in Eq. (6) to determine the displacement field u. The strain energy terms and stress tensors are related according to In this section, we summarize the four strain energy decomposition methods that are compared in this work. Throughout the rest of this work, the four methods are called strain deviatoric decomposition (StrainDe), strain spectral decomposition (StrainSp), stress spectral decomposition (StressSp), and stress deviatoric decomposition (StressDe).

Strain deviatoric decomposition (StrainDe)
In the decomposition method from Amor et al. (2009), ψ + and ψ − are separated based on the volumetric and deviatoric parts of the strain tensor. Deviatoric and positive volumetric parts contribute to fracture, i.e. strains that change the shape but not the volume and strains that cause the volume to increase. Negative volumetric parts do not contribute to fracture, i.e. strains that cause the volume to decrease. This is expressed as where K n = λ + 2μ n is the n-dimensional bulk modulus, tr(ε) = ε 11 + ε 22 + ε 33 , the positive and negative operators are defined as x ± = x±|x| 2 , and ε dev = ε − 1 3 tr(ε)I with the second-order identity tensor I. The positive and negative stresses are defined as

Strain spectral decomposition (StrainSp)
The decomposition method from 2010) uses a spectral decomposition of the strain tensor, ε = a ε a n a ⊗ n a , where ε a and n a are the ath eigenvalue and eigenvector, respectively. Positive eigenvalues and a positive sum of all eigenvalues contribute to fracture evolution, while negative eigenvalues and a negative sum of all eigenvalues do not. Positive eigenvalues are associated with tensile strain and negative eigenvalues with compressive strain. However, shear strain results in both positive and negative eigenvalues and therefore only partially contributes to fracture using this method. The positive and negative parts of the strain energy density are defined as and the positive and negatives stresses are defined as

Stress spectral decomposition (StressSp)
The decomposition method from Zhang et al. ) uses a spectral decomposition, as done in strain spectral decomposition ), but it is applied to the undamaged stress σ 0 from Eq. (5) rather than the strain. Thus, σ 0 = a σ a n a ⊗ n a , or in its tensor form, where = diag(σ 1 , σ 2 , σ 3 , ) is a diagonal matrix containing the three eigenvalues and Q =[ n 1 , n 2 , n 3 ] is a matrix containing the three corresponding eigenvectors. In this method, positive eigenvalues of the stress contribute to fracture and negative eigenvalues do not. Similar to the strain spectral method, tensile stresses have positive eigenvalues, compressive stresses have negative eigenvalues, and shear stresses have both positive and negative eigenvalues. The positive and negative stresses are defined as where + contains the positive eigenvalues and − the negative, i.e. ± = diag( σ 1 ± , σ 2 ± , σ 3 ± ). The positive and negative parts of strain energy density are defined as A benefit of this method is that it does not assume that the elastic response is isotropic. The elasticity tensor C can have any crystallographic symmetry.

Stress deviatoric decomposition (StressDe)
For this study, we introduce a new decomposition method that is similar to the strain deviatoric decomposition (Amor et al. 2009). In this method, the undamaged stress σ 0 from Eq. (5) is decomposed into deviatoric and hydrostatic parts. Deviatoric and positive volumetric parts of the stress contribute to fracture, and negative volumetric parts do not. The positive and negative stresses are determined according to where The positive and negative parts of the strain energy density are then defined as As with the stress spectral decomposition, this decomposition can use an elasticity tensor C with any crystallographic symmetry.

Numerical implementations for fracture models
In order to focus our study on the impact of the strain energy decomposition method on the fracture behavior, we have implemented the four decomposition methods in the same phase field fracture code. We use the finite element method to solve the coupled governing equations. For quasi-static simulations we solve Eqs. (3) and (7) and for dynamic simulations we solve Eqs. (3) and (8). In the dynamic fracture simulations, the mass and stiffness matrices are constructed for the simulations: where M is the mass matrix, K is the stiffness matrix, f (t) is the external load, and damping is neglected during the dynamic fracture process. The HHT-α integration method (Hilber et al. 1977;Miranda et al. 1989) is adopted in this simulation for time integration, which is second-order accurate. The specific process for solving the finite element dynamic problem can be found in the original papers (Hilber et al. 1977;Miranda et al. 1989).
The phase field fracture model was implemented using the Multiphysics Object-Oriented Simulation Environment (MOOSE) (Gaston et al. 2009;Tonks et al. 2012;Chakraborty et al. 2016;Permann et al. 2020). MOOSE solves the coupled equations using the finite element method with implicit time integration for both quasi-static and dynamic fracture problems. Newton's method is used to solve the resultant nonlinear problem. For all four decomposition methods evaluated in this work, the positive and negative parts of the strain energy density are used to calculate the history variable in Eq. (2) at each time step, which is substituted into Eq. (3) for phase field variable  7) for quasistatic fracture and Eq. (8) for dynamics fracture. For all the example cases presented in this work, all simulations use the same parameters and solution approaches and the only difference is the strain decomposition method. This ensures that any differences between the methods are due to the methods themselves and not the code used for the simulation.
The values for the crack width parameter l for each simulation was chosen to provide a reasonable balance between accuracy and efficiency. The element size used to resolve the crack was always ≤ l/2, which was determined to be a good resolution by Chakraborty et al. (2016).

Mechanical response of cracked material
As discussed above, the strain energy decomposition method also decomposes the stress into the portions that can be supported by cracked material and the portions that cannot. Thus, simply calculating the stress field in a cracked material under various applied loads provides a simple means to begin comparing the four methods.
We apply displacements to fully cracked material in quasi-static conditions to induce tension, compression, and simple shear stresses. The cracked sample is shown in Fig. 1, where the crack is represented with phase field parameter d = 1, with a crack thickness of 0.01 mm.
In these simulations, we only determine the response of the cracked material by solving Eq. 7, and do not evolve the crack variable field d. The material properties used in this section are taken from  and are shown in Table 1.

Tension tests of cracked materials
First, we determine the strain distribution and max tensile stress under tension, where the y-displacements are fixed on the bottom boundary, the x-displacements are fixed on the left and right boundaries, and the top boundary has a fixed y-displacement of 0.01 mm.
Since our material has a horizontal crack across its entire width, an applied tensile displacement should just separate the two regions divided by the crack; thus, all of the strain should be localized in the cracked region and no tensile stress should be generated. The strain ε yy distributions across the domain for the four decomposition methods are shown in Fig. 2a-d. From the figure, we can see the max strains are close to 1 and are localized in the cracked region. To quantify the stress within the material, we plot the maximum stress throughout the domain from each method on a bar plot, shown in Fig. 2e. All methods predict a small stress within the material (around 0.28 MPa, due to the value k = 1.0 × 10 −6 used in the degradation function), since the strain is fully localized in the crack and the cracks cannot support substantial tensile stress. These results indicate that the four methods predict the correct mechanical behavior for a cracked material in tension.

Compression tests of cracked materials
Second, we determine the mechanical behavior in the cracked material under a compressive load. In order to capture the pure compressive load, the y-displacements are fixed on the bottom boundary, the x-displacements are fixed on the left and right boundaries, and the top boundary has a fixed y-displacement of -0.01 mm, where the sign means a compressive displacement load.
Even though the material has a horizontal crack across its entire width, the crack can still support a compressive load such that the material should distribute the strain across Fig. 2 Mechanical response of the cracked material under tension predicted by the four methods: (a) -(d) strain distribution for the strain deviatoric, strain spectral, stress deviatoric, and stress spectral methods, respectively; (e) maximum tensile stress throughout the domain the material as if it were not cracked. The strain ε yy distribution and the minimum compressive stress for these four methods are shown in Fig. 3e. The compressive strain distributions show that the two deviatoric methods have localized strain in the cracked region, while the two spectral methods have strain that is uniformly distributed. In addition, the compressive stresses are predicted to be much larger by the spectral methods than by the deviatoric methods (Fig. 3e), since the deviatoric methods cannot transmit the majority of the compressive stress through the cracked material. Thus, the spectral decomposition methods predict the mechanical behavior under compression more accurately than the deviatoric methods.

Shear tests in cracked materials
Third, we determine the mechanical behavior in the cracked material under a shear load. To capture simple shear, the y-displacements are fixed on all four boundaries, the x-displacements are fixed on the bottom boundary, and the top boundary has a fixed x-displacement of 0.01 mm.
Since our material has a horizontal crack across its width, an applied shear displacement should slide the top region along the bottom region. Thus, the shear strain should be localized in the cracked region and no shear stress should be generated. The shear strain ε xy distributions for the four methods are shown in Fig. 4a-d. The shear strain for the two deviatoric methods is localized in the cracked region as expected, while the two spectral methods have the largest shear strain in the cracked region but some shear strain throughout the material. The maximum shear stress σ xy for the four methods is shown in Fig. 4e; the deviatoric methods have no stress, while the spectral methods have significant stress. In addition, the two spectral methods have different values for the maximum stress; this is because the manner in which they use the positive and negative eigenvalues to determine σ + and σ − are different (see Eqs. (17) and (19)). Thus, the deviatoric

Quasi-static fracture examples
In this section, we conduct quasi-static fracture simulations with cracked samples by solving the phase field evolution equation, Eq. (3), and the quasi-static equilibrium equation, Eq. (7).

Single edge notched tension test
The first example is the single edge notched tension test (SENT), where we apply all four methods to study the fracture behavior under tension; the sample geometry and boundary conditions are shown in Fig. 5. The domain size is 1 mm × 1 mm, with an increasing positive y-displacement applied across the top boundary and the x-and y-displacements fixed on the bottom boundary. The material properties used for the simulations are shown in Table 1.
The final crack paths for the four methods are shown in Fig. 6a -d. As shown in the figures, the crack paths are all very similar, indicating that all decomposition methods work well in mode I fracture. To provide a quantitative comparison, we plot the crack length as a function of the applied strain, as shown in Fig. 6e. The two deviatoric methods predict crack nucleation at a slightly smaller strain than the two spectral methods. The stress spectral method predicts faster crack propagation than the strain spectral method, such that it approaches and eventually reaches the crack length predicted by the two deviatoric methods. Thus, the strain spectral method predicts slower crack propagation than the other three methods. It is not surprising that the four methods predict similar behavior in a tensile test, since they all predict the correct response for cracked material under a tension load, as shown in Fig. 2.

Single edge notched shear test
The second example is a single edge notched sample under shear loading, where the geometry and boundary conditions are shown in Fig. 7. The domain size is 1 mm × 1 mm, with increasing x-displacement applied across the top boundary, as well as no y-displacement, and the x-and y-displacements fixed on the bottom bound-  ary. The material properties are the same as we used for the previous example (see Table 1).
The crack evolution under a shear load is shown in Fig. 8a -e; the crack profile and propagation varies significantly between the various methods, so we include the crack paths at various times. We also plot the crack length versus the x-displacement of the top surface for the four methods in Fig. 8f. At a displacement of 0.007 mm (Fig. 8a), the cracks from the two deviatoric decomposition methods and the stress spectral method are identical, but the crack from the strain spectral method is much shorter. The cracks from the deviatoric methods and the stress spectral methods remain similar up to a displacement of about 0.0085 mm. The stress deviatoric method then predicts acceleration of the crack and formation of a second crack along the bottom boundary, as shown in   The two cracks meet and the crack evolution stops at a displacement of 0.0094 mm. The strain deviatoric and the stress spectral methods continue to predict similar crack propagation up to displacement of 0.0097 mm, at which point the strain deviatoric method then predicts the rapid propagation of the crack until it meets a crack that forms along the bottom boundary and stops propagating (Fig. 8c). The crack from the strain spectral method has nearly caught up to the crack from the stress spectral method at a displacement of 0.015 mm (Fig. 8d). Once the crack from the stress spectral method reaches the bottom boundary, it then extends slightly along the boundary. Though its propagation is slower, the final crack from the strain spectral method is very similar to that from the stress spectral method (Fig. 8e). In the earlier shear example, the deviatoric methods predicted the behavior of the cracked material under a shear load more accurately than the spectral methods, as shown in Fig. 4e, but in this case the spectral and deviatoric methods predicted fairly similar behavior.

Three-point bending fracture
The third quasi-static example is a three-point bending simulation, which is a common test used for brittle materials and has been analyzed theoretically (Bittencourt et al. 1996). It has also been modeled using phase field fracture simulations Ambati et al. 2015a).
The geometry and boundary conditions used for our three-point bending simulations are shown in Fig. 9. The domain size is 8 mm × 2 mm and the x-and y-displacements are pinned at single nodes at the bottom left and bottom right corners, and a displacement load is applied at a single node at the center of the top boundary. The material properties used in the three-point bending fracture simulations are taken from  and are shown in Table 2. To reduce the computational cost, we use a heterogeneous mesh. We use a fine square element of size 0.01 mm in a 0.2 mm wide strip at the center of the domain, marked in red in Fig. 9. The rest of the mesh has a coarser element size of 0.025 mm × 0.01 mm. The time step dt = 5 × 10 −5 s is used for the simulations. The final crack paths using the four methods are shown in Fig. 10a -d. Both spectral decomposition methods predict propagation of a crack from the notch vertically upward, as is expected for this test. Damaged regions also form at the three points where the boundary conditions are applied, though these damaged regions are larger for the strain spectral method than for the stress spectral method. The crack length grows faster with displacement and reaches a higher value with the strain spectral method than with the stress spectral method, as shown in Fig. 10e. However, this is because the crack length is also accounting for the damaged regions that form where the load and boundary conditions are applied; the length of the crack starting at the notch is the same for the two methods. There is no crack propagation at the notch with the two deviatoric decomposition methods, only damage forming around the node where the load is applied. The compressive strain is localized at this damaged region, and does not propagate through to the rest of the material so that a crack does not form at the notch. We repeated the simulations with much smaller time steps (dt = 1 × 10 −5 s and dt = 1 × 10 −6 s) and using a uniform mesh, but the results were unchanged. This localization of the strain at the location where the compression load is applied predicted by the deviatoric methods is likely related to our finding in the cracked material under a compression load in which the deviatoric methods do not transmit the majority of the compressive stress through the cracked material (Fig. 3).
In our three-point bending simulation, the load and boundary conditions were applied at a single node, which is not physically accurate. In reality, the load is applied and the sample is supported across some finite width of material. To more accurately account for Fig. 10 Fracture results from the three-point bending simulations with the load and boundary conditions applied at a single node: (a) -(d) final crack paths for the strain deviatoric, strain spectral, stress deviatoric, and stress spectral methods, respectively. (e) crack length evolution for the four methods the actual conditions of such a test, we repeat the simulations but apply the top load and bottom boundary conditions at three, five, and seven nodes, rather than at a single node. The deviatoric methods still do not predict crack propagation at the notch when three and five nodes are used, but they do with seven nodes. The final crack paths from the four decomposition methods using seven nodes are shown in Fig. 11a -d and the crack length evolution comparison for these four methods are shown in Fig. 11e. The crack paths for the two deviatoric methods are identical, with the crack propagating through the full thickness of the material, and the crack length evolves identically. The deviatoric methods have early crack formation at the bottom boundary conditions and then initiation of the center crack at a displacement slightly larger than 0.04 mm. The spectral methods predict center crack initiation at around 0.05 mm and the center crack never reaches the top surface. The strain spectral method predicts large cracked regions around the bottom boundary conditions and where the load is applied on the top surface, so that it predicts the largest crack length. Thus, with the load and boundary conditions applied at seven nodes, the deviatoric methods show slightly better performance in the three-point bending simulation than the stress spectral method and much better than the strain spectral method.

L-shaped panel test
As the final quasi-static example, we simulate a mixed-mode failure in an L-shaped panel, where the experimental results are taken from Winkler (2001). It has been previously modeled using the phase field fracture method (Ambati et al. 2015a;Wu et al. 2020).

Fig. 11
Fracture results from the three-point bending simulations with the load and boundary conditions applied at seven nodes: (a) -(d) final crack paths for the strain deviatoric, strain spectral, stress deviatoric, and stress spectral methods, respectively. (e) crack length evolution for the four methods The geometry and boundary conditions for the simulation and experimental results are shown in Fig. 12. The x-and y displacements are fixed on the bottom boundary and a point load is applied as indicated in Fig. 12b. The material properties are taken from (Ambati et al. 2015a) and are shown in Table 3. To reduce the computational cost, we use adaptive mesh refinement, where the element size varies from 3.0 mm to 0.625 mm.
The final crack paths for the four methods are shown in 13a -d. From the results, we can see that no cracks form with the deviatoric decomposition methods, but a crack does form from the inner corner with the two spectral methods. The crack propagates nearly straight outward from the corner with the strain spectral method and it curves and eventually propagates horizontally with the stress spectral method. The crack path from the stress spectral method compares better with the experimental results (see Fig. 12(b)). The evolution of the crack length with time is shown in Fig. 13e; the crack propagation rate predicted by the two spectral methods is identical.

Dynamic fracture examples and results
In this section, we continue our comparison of the four strain energy decomposition methods using two dynamic fracture examples.

Tensile dynamic fracture
For the first dynamic fracture example, a single edge notched plate is modeled under a dynamic tensile load. The geometry and boundary conditions are shown in Fig. 14. A traction load is applied with an applied y-displacement on the top and bottom boundaries. This tensile dynamic fracture example has been studied numerically (Song et al. 2008;Ha   and Bobaru 2010; Borden et al. 2012) and therefore has a known solution. The material parameters used in this simulation are taken from (Borden et al. 2012) and are shown in Table 4. Also, plane strain is assumed. A uniform mesh with element size h = 0.25 mm is used. The final crack paths are shown in Fig. 15a -b and the length of the crack (only including the upper branch after branching) versus time is shown in Fig. 15e. The crack paths from the four methods are all very similar. The crack length evolution from the four methods is also very similar. These results are consistent with those shown in Figs. 2 and 6, verifying that the four decomposition methods all perform well with tensile loading.

Dynamic shear fracture example
For the second dynamic fracture example, we use a shear fracture configuration that has been tested experimentally by Kalthoff (Kalthoff and Winkler 1988;Kalthoff 2000). It has also been studied numerically (Song et al. 2008;Borden et al. 2012). The geometry and boundary conditions are shown in Fig. 16, where the y-displacements are fixed on the bottom boundary, and a velocity boundary condition is applied in the x-direction on the left boundary below the initial crack, where v = 16.5t m/s when t < 1 × 10 −6 s and v = 16.5 m/s when t ≥ 1 × 10 −6 s. The material properties are taken from (Borden et al. 2012) and are shown in Table 5, plane strain is assumed, and a uniform element size of h = 0.125 mm is used. The final crack paths for these four models are shown in Fig. 17a -d. The experiments found that the crack propagated from the initial crack tip at an angle of approximately 70° (Kalthoff 2000). All four methods predict the formation of a crack somewhat similar to that seen in the experiments; the crack from the deviatoric methods has an angle of around 64°and the crack from the stress spectral method has an angle of around 67°. Fig. 15 Fracture results from the dynamics tensile fracture simulations: (a) -(d) final crack paths for the strain deviatoric decomposition method, strain spectral method, stress deviatoric decomposition method, and stress spectral method, respectively; (e) the crack length evolution for the four methods The crack from the strain spectral method has an initial angle of about 71°, but the crack curves and when it reaches the top boundary it has an angle of 56°. The deviatoric decomposition methods also predict the formation of a secondary crack below the initial crack tip which was not seen in the experiment. The crack length evolution is shown in Fig. 17e. The two deviatoric methods predict crack nucleation sooner than the spectral methods, though the initial crack is actually the secondary downward crack. Once the crack propagates upward, it propagates faster than predicted by the spectral methods. Also, the final crack lengths predicted by the deviatoric decomposition methods are longer than those predicted by the two spectral methods due to the secondary crack. The final length of the crack from the strain spectral method is longer than that from the stress spectral method, due to its curvature. As with the other shear tests (Figs. 4 and 8), the deviatoric methods predict different behavior than the spectral methods. In this case, the crack path predicted by the stress spectral method is the most similar to the experimental crack path (Kalthoff 2000).

Discussion
It is clear from the results from the various cases we have investigated in this work that the strain energy decomposition method used in a phase field fracture model can have a significant impact on the fracture results. This impact causes different mechanical responses of cracked material to applied loads, different fracture rates, and different crack paths. In this section, we discuss the results from the stress and strain decomposition methods. We also highlight how the different approaches perform for material in tension, compression, and shear.  The original decomposition methods from Amor et al. (2009) and 2010) both acted on the strain, and nearly all other work from the phase field fracture literature decompose the strain. However, decomposing the stress makes it much easier to include elastic anisotropy ) and could facilitate the use of other constitutive models besides linear elasticity. The stress spectral method has been introduced previously ), but the stress deviatoric method is introduced for the first time in this work. Our results indicate that whether the strain or the stress is decomposed does not have a significant impact on the predicted behavior. The two deviatoric methods are unable to propagate compressive stress through a crack, while the spectral methods are able to. The two spectral methods are unable to localize the shear strain in the crack, while the two deviatoric methods are able to. In addition, the crack paths and crack length versus displacement curves from the single edge notched shear tests are similar for the deviatoric methods and for the spectral methods, but are distinct from each other. The deviatoric methods predict no crack propagation in three-point bending simulations until the load was applied at seven nodes, while the spectral methods predict crack propagation no matter how the load is applied. The deviatoric methods did not predict crack propagation in the mixed-mode L-shaped panel case, while the two spectral methods did (though stress and strain spectral methods predicted different crack paths). In the dynamic shear test, the deviatoric and spectral methods again predict similar behavior but distinct from each other. Our results here clearly indicate that the stress can be decomposed without significantly impacting the predicted fracture behavior. When using a deviatoric decomposition method, our new deviatoric decomposition of the stress should be considered as a possible option, especially if elastic anisotropy is used.
Our results demonstrate that when modeling the fracture behavior under a tensile load, all the methods perform well and either spectral or deviatoric methods provide reasonable results. Our results in Fig. 2 show that all the methods correctly localize tensile strain within the crack. They also all predict similar behavior in a quasi-static tensile test (Fig. 6) and in a dynamic tensile test (Fig. 15).
When modeling the fracture behavior under a compressive load, our results indicate that spectral decomposition will provide more accurate predictions. Our results in Fig. 3 show that the spectral methods accurately allow cracked material to fully support a compressive stress, while the deviatoric methods only support a portion of the compressive stress. Also, the spectral methods accurately predict the fracture in a three-point bending test (Fig. 10), while the deviatoric methods do not accurately predict the behavior unless the applied load and boundary conditions are applied at seven nodes (Fig. 11).
With a shear load, our results indicate that both deviatoric and spectral methods provide reasonable predictions, though both approaches do predict inaccurate behavior in one case. Our results in Fig. 4 show that the deviatoric methods accurately localize the shear strain within the crack, while the spectral methods allow the crack to support some shear stress. Both the spectral and deviatoric methods predict similar crack paths in a quasi-static shear test, though the deviatoric methods predict much faster crack propagation for the later half of the crack (Fig. 8). Both spectral and deviatoric methods also predict similar behavior in a dynamic shear test, though the deviatoric methods predict a secondary crack that was not seen in experiments (Fig. 17).
For the L-shaped panel mixed mode case, the spectral methods provide reasonable predictions, while the deviatoric methods do not predict any crack formation. Thus, in cases where the load type is not known, our results indicate that a spectral method has a better chance of giving reasonable predictions.
Though all four decomposition methods performed well in many of the cases investigated in this work, none of them performed well in all cases. Thus, it is clear that new strain energy decomposition methods could be a valuable contribution to phase field fracture models. The benchmark cases we have used in this work should be considered as a means of evaluating any new decomposition methods, to allow the results to be directly compared to the four methods we have assessed.

Conclusions
This paper provides a comparative study of four strain energy decomposition methods for phase field fracture models. Three of the methods come from the literature, and one, the stress deviatoric method, is new in this work. The methods are compared by comparing the predicted mechanical behavior of cracked material, three quasi-static fracture cases, and two dynamic fracture cases.
In general, all four methods perform well in tension, the two spectral methods perform better in compression, and all the methods perform reasonably well in shear. All four methods correctly predict all of the tensile strain localized within the cracked material, the spectral methods correctly predict cracked material can support compressive stress, and the deviatoric methods correctly predict localization of shear strain within cracked material. From the quasi-static and dynamic fracture cases, all the methods predict similar behavior for cases under tension, including both the crack path and the fracture propagation rate. The two spectral methods accurately predict the crack path in a three-point bending test, while the two deviatoric methods require that the load be applied across at least seven nodes for the crack to propagate. The two spectral methods predict crack formation in the mixed mode L-shaped panel fracture case, though the stress spectral method results better match experimental results. For the quasi-static shear case, the two spectral methods predict similar crack paths but the stress spectral method predicts slower crack propagation. The two deviatoric methods predict similar crack paths, though they are somewhat different than those from the spectral methods. For the dynamic shear case, the two spectral methods again predict the correct crack paths and similar crack propagation rates. The two deviatoric methods also predict the correct crack paths and a similar rate, but predict a secondary crack below the initial crack tip. In all cases, whether the strain or stress is decomposed does not have a large impact on the model predictions.