Numerical Solution Strategies for a Dynamic Phase Field Fracture Model

This work discusses the efficiency of six strategies for the numerical solution of the coupled system of partial differential equations that arise from a phase field description of dynamic fracture. Efficient numerical treatment of the dynamic phase field fracture model is a prerequisite for the simulation of failure due to brittle fracture in realistic scenarios such as manufacturing.Firstly, the phase field description of fracturing of brittle solids is introduced. Afterwards, three monolithic as well as three staggered finite element solution strategies are outlined and their performance is studied in two benchmark problems.


Introduction
Fracturing of solids has been an important aspect of continuum mechanics for almost a century. Starting with the classical theoretical works [12] and [9], engineers and scientists have developed the theoretical framework of fracture mechanics in order to prevent structural failure or to gain insight into the physical process of fracturing itself. Results have also been obtained for the dynamic case in which inertial effects are accounted for, e.g. see [8].
With the rise of computer technologies, the numerical simulation of fracture has gained importance in the more recent past. These simulations enable the application of concepts of classical fracture mechanics to more complex engineering problems for which no analytic solution exists.
Phase field models have first appeared in the context of phase transition processes where they allow to describe various processes in heterogeneous continua, that range from solidification [3] and phase transformation in a solid [18] to the modelling of ferroelectric materials [19]. The phase field approach to fracture has gained popularity since the late 1990s, starting with Bourdin's regularized description [6] of the variational formulation of brittle fracture [7]. Phase field models of fracture usually only have one phase field, that indicates the state of the material, i.e. whether the material is fractured or undamaged and thus approximates the crack. The fracture process, i.e. the evolution of the phase field over time, is fully governed by a set of coupled partial differential equations, the equation of motion and a phase field evolution equation. This renders tracking of the crack propagation, which is required in other methods, unnecessary. Consequently, the numerical solution of this set of equations is straightforward and can be handled with standard finite element methods. Phase field models have shown some promising results for the simulation of dynamic fracture, see [11], [5] and [17].
However, the advantages of the phase field approach to fracture also come with costs. Since the smooth transition zone in which the phase field varies from broken to unbroken material has to be small relative to the dimensions of the considered body, a very fine spatial discretization is required, at least locally, to resolve this transition properly. Furthermore, in dynamic simulations wave propagation needs to be treated adequately and hence reasonable time step sizes are small. This means that the computational effort for dynamic fracture phase field simulations can become relatively high.
It is therefore of interest to find solution strategies that allow for an efficient numerical treatment of dynamic phase field fracture models. In the literature, there are mainly two approaches. The first way is to solve the two coupled equations simultaneously for the two primary unknown fields, i.e. the displacements and the phase field. This so-called monolithic strategy promises to find the, up to a predefined tolerance, exact solution at each discrete time step. Here, implicit schemes have to be used and the computational effort per time step is relatively large, because of the nonlinearity of the coupled system. Alternatively, staggered solution strategies can be used, see e.g. [11] and [1]. In these strategies, the problem is decoupled and solved in two steps. First, the equation of motion is solved keeping the phase field fixed and subsequently the evolution equation is solved for the phase field variable fixing the displacements. In this case nonlinearities that come from the coupling of the two fields are removed and consequently the subproblems can be solved more efficiently. However, the two steps of the staggered algorithm may have to be repeated several times in order to obtain convergence, i.e. a solution that is insensitive to further so called 'staggered iterations'. A criterion to judge convergence of the staggered algorithm in this sense is proposed in [1]. Nevertheless, it is rather common to perform only one staggered iteration. For such a one-step staggered approach it is necessary to keep the time step rather small which might be required anyway in highly dynamic simulations. Furthermore, the decoupling allows for efficient explicit time integration of the equation of motion which is not possible in the monolithic case.
This work studies the performance of six numerical solution strategies for dynamic phase field fracture problems. The first section introduces the phase field description of dynamic fracture. Subsequently, the chosen numerical solution strategies, three monolithic and three staggered strategies, are presented. Finally, the performance of the schemes is analyzed by means of two benchmark problems.

A Phase Field Model for Dynamic Fracture
We consider a body Ω ⊂ R n , n ∈ {1, 2} with external boundary ∂Ω that is regarded as homogeneous, linear elastic material with Lamé parameters λ and µ and mass density ρ, see Fig. 1a). The deformation and the internal forces of this body are described in terms of the displacement field u (x, t) and the Cauchy stress tensor σ (x, t). The fields have to satisfy Dirichlet boundary conditions on ∂Ω u and traction boundary conditions on ∂Ω t , where n is the outward directed normal vector on the boundary ∂Ω = ∂Ω u ∪ ∂Ω t , as well as initial conditions andu where the dot( * ) = ∂( * ) ∂t indicates the material time derivative. The linearized strain tensor serves as a strain measure. The body may also contain cracks with the associated crack surface Γ.
In order to approximate Γ, a phase field order parameter s(x, t) ∈ [0, 1] is introduced which varies continuously from s = 1 in undamaged material to s = 0 in fully broken material. This smooth representation of a crack is illustrated in Fig. 1b). Existing cracks in the material can be modelled by specifying initial conditions s(x, t 0 ) = 0 (6)  Fig. 1: a) Fractured body with internal discontinuity (sharp crack) Γ b) and smooth representation of the crack by means of a phase field s(x, t).
for the phase field. Following Griffith's idea, the fracture surface Γ is associated with a fracture energy that is dissipated during the creation of the new surface. The phase field is used to approximate this fracture energy in a regularized way, i.e.
with the fracture resistance G c being a constant material parameter. The length-scale parameter ϵ controls the width of the phase field approximated crack. An increase of ϵ causes the width of the transition zone between broken and unbroken material to increase as well. It can be shown that the volume integral on the right-hand side of approximation (7) converges to the surface integral on the left-hand side if ϵ → 0. In order to model the degradation of stiffness in broken material, the phase field s is linked to the elastic energy of the body, i.e.
The constitutive behaviour is determined by the definition of the strain energy density The strain energy density is decomposed in a crack driving, strain part ψ e + that is affected by the degradation function g(s) and a part that is associated with compressive strain states ψ e − . In the literature, there are several propositions to implement this decomposition, which all aim to prevent unphysical fracture behaviour in compressive load states. The model of [2] proposes a split based on a volumetricdeviatoric decomposition of ε whereas [20] actually take the orientation of the crack into account. The approach used in this work follows the publication of [16]. It is based on an eigenvalue decomposition of the elastic strains, such that where ε i is the i-th eigenvalue of ε and n i is the associated normalized eigenvector. With and and Eventually an expression for the negative strain energy and the positive strain energy is formulated as and where the operator tr(·) indicates the trace of a second order tensor and the symbol : is the double contraction of two second order tensors. Utilizing the fracture field s, the degradation function g(s) allows to model the loss of stiffness in broken material by reducing the strain energy accordingly, i.e. it has to satisfy g(1) = 1 and g(0) = 0. The compressive strain energy is not affected by g(s) which models the impenetrability of cracks during crack closure, i.e. no degradation of the 'compressive' stress see (17). An additional consequence of leaving ψ e − unaffected by s is that crack growth is not driven by compressive load states, which will become more apparent in (24). The stress is The wave speed of dilatational waves in the considered elastic medium is The total kinetic energy of the body is assumed not to be affected by the phase field, i.e The work of external forces acting on the boundary ∂Ω reads Finally, the dynamic fracture problem can be stated using Hamilton's principle for arbitrary times t 1 < t 2 , where the Lagrangian is given by For the problem at hand, the Euler-Lagrange equations following from (21) are the equation of motion Applied Mechanics and Materials Vol. 869 103 32 Physical Modeling for Virtual Manufacturing Systems and Processes and the phase field equation as well as the Neumann boundary conditions for the displacement and for the phase field ∇s · n = 0 on ∂Ω.
In (24) it becomes apparent, that the property g ′ (0) = 0 is important to eliminate the driving source term for crack growth when s = 0. In this work, the degradation function is chosen to be A degradation function of this type was first proposed by [4]. An additional constraint on the phase field s is necessary to impose the irreversibility of fracture which is not modelled by the the above equations. We define the irreversibility constraint on the crack field. Herein, t * x is the time when the crack field becomes zero at the location x for the first time. The model extension (28) prevents unphysical crack healing, see [13] for details.

Numerical solution strategies for the dynamic phase field fracture problem
In order to solve fracture problems numerically, the initial boundary value problem described by the field equations (23) and (24) together with the respective initial and boundary conditions is spatially discretized with quadrilateral bilinear, i.e. for n = 2, or linear bar, i.e. for n = 1, isoparametric finite elements. This results in a coupled system of ordinary differential equations of the type where the arrays of unknowns contains all nodal displacement and phase field variables of all finite element nodes I = 1, ..., N . Details on the spatial discretization can be found in [17]. E.g. in two dimensions, the 2N ×2N matrix M uu is the global mass matrix and the 2N ×1 arrays P u , F u and R u denote the contributions of internal and external "forces" to the global residual of the spatially discretized equation of motion. Similarly, the N × 1 quantities P s and R s denote the global residual of the spatially discretized phase field equation (24). An alternative notation of (29) is where the arrays of unknowns are summarized is the solution of (29) and

Staggered strategy
For the solution of (29), two main strategies exist. The first option is a staggered approach, where the solution is decomposed into two stages. Let (u 0 , s 0 ) be the initial guess for a solution of (29). The actual solution can be found by the staggered algorithm summarized in Fig. 2, which splits the coupled problem in two subproblems that are solved consecutively.
One advantage of this staggered algorithm is that it results in two possibly linear subproblems. The equation of motion is nonlinear with respect to the displacements due to the split of the strain energy (17), which causes the stress to be nonlinearly dependent on the strain state. Hence the equation is linear, whenever ψ e − = 0 or ψ e + = 0, i.e. all strain eigenvalues are either positive or negative. The phase field equation is linear in s whenever g ′ (s) is linear, e.g. g(s) is a polynomial of order two at most. This is however not the case for the present cubic degradation function. Even in the nonlinear case the sub problems are more robust because further nonlinear features, e.g. see (17), are removed by freezing one of the fields. Furthermore, efficient explicit time integration schemes can be used to solve the second order spatially discretized equation of motion, which is not possible for the coupled equations because of the lack of a time derivative in the discretized evolution equation. In principle, a viscous approximationṡ could be used instead of the stationary phase field equation. However, explicit time integration of the coupled system would still be impractical because of the strict stability requirements of the first order viscous phase field equation, i.e. ∆t max ∼ h 2 where h is a characteristic length of the smallest element, see textbooks like [23]. The accuracy of one iteration of the staggered algorithm depends on how well s k approximates s k+1 . If the assumed field s k is close to the final field s k+1 , the displacements u k+1 will also be a good approximation of the correct solution. Thus, one staggered iteration might be enough to find Applied Mechanics and Materials Vol. 869 105 34 Physical Modeling for Virtual Manufacturing Systems and Processes an accurate solution whenever the time step is small. In general however, several iterations of the staggered algorithm need to be performed until the algorithm converges. In [1] it is proposed that the algorithm can be assumed to be converged if the solution is insensitive to further iterations of the staggered algorithm. [1] propose the total energy E e +E s as a quantitative measure to judge the change of the solution per staggered iteration. We follow the idea presented in [1] but use the contribution of inertial and internal forces to the global residual I(u k+1 , s k+1 ), see (31), as a global quantitative measure to judge the change of the solution (u, s) per staggered iteration. First, the change of I in the current staggered iteration is normalized by the total change of I in all previous staggered iterations, i.e.
The number of iterations is normalized as well Eventually the quantity can be computed to judge convergence in the form of the criterion where γ tol is chosen to be 1 The proposed criterion requires at least three staggered iterations per time step, because information of the previous iteration is needed, i.e. there is no convergence check in the first iteration, and convergence cannot be achieved in the second iteration, ∆I 2 norm = 1.0. Typical graphs of ∆I k norm and γ k are shown in Fig. 3.

106
Physical Modeling for Virtual Manufacturing Systems and Processes Applied Mechanics and Materials Vol. 869 35 Alternatively to the staggered strategy, (31) may be solved in a single step for (u, s). In this case, the full coupled problem has to be solved and implicit time integration is most practical. This requires the solution of the complete nonlinear system of algebraic equations in every time step and the computational effort per time step is larger than for each of the subproblems in the staggered case. Utilizing an adequate time integration scheme the algorithm can be set up such that it is unconditionally stable, i.e. there are no strict numerical bounds for the time step size. However, the larger time step has still to be justified physically, i.e. it needs to be small enough to resolve wave propagation and possibly fast crack growth.
The monolithic algorithm is also not as robust as the staggered algorithm due to the inherent nonlinearity of the coupled phase field fracture problem. This means the Newton algorithm will not necessarily find a solution if the time step and therewith the changes in the computed quantities are too large. Hence, it is reasonable to introduce an adaptive time step size control to enhance the robustness of the scheme in situations where a small time step is necessary, e.g. fast crack propagation, and to exploit the numerical efficiency of larger time steps whenever possible.

Time integration
For both the monolithic as well as the staggered strategy we consider time discrete approaches. Thus discrete time steps are introduced recursively as where t n , t n+1 and ∆t are the previous discrete time, the current discrete time and the current time step, respectively. By this, (29) at time t n+1 becomes [ In order to express the current nodal 'accelerations' in terms of current unknowns, time integration schemes are employed. In this work, we will only consider one-step time integration schemes which only require knowledge of quantities from the previous time step and not of even earlier time steps. As mentioned before, implicit time integration schemes are the adequate choice, if the monolithic solution strategy is used, but can be used in a staggered solution strategy as well. In this work, two different time integration schemes are considered. The Newmark method uses with parameters 0 ≤ 2β ≤ 1, 0 ≤ γ ≤ 1 for time integration. In [21] the amplitude error and stability have been investigated for linear problems. The amplitude error is zero for γ = 0.5 for problems with no damping. In case of damping effects the amplitude error is of order O(∆t 2 ). Stability is ensured for 2β ≥ γ ≥ 0.5 for linear problems. The first order accurate fully implicit Euler scheme uses the following approximations However, it envolves severe numerical dissipation for larger time steps. ) is the solution of (41) Fig. 4: Newton-Rapshon algorithm. The residual norm ∥ * ∥ rel is a relative quantity that allows the tolerance tol to be chosen independently of the problem size.
The respective equations gained from time integration can be substituted in (41) to yield a nonlinear system of algebraic equations for the current unknowns which is solved with the Newton-Raphson algorithm summarized in Fig. (4). The computational effort of the Newton-algorithm per iteration is dictated by the computation of the residual, i.e. step 1, the computation of the tangent matrix, i.e. step 2 and the solution of a set of linear equations. The amount of necessary Newton iterations is determined by the nonlinear features of the considered problem. The previously mentioned automatic time step control that is utilized in this work adjusts the time step according to the number of iterations that was needed in the previous time step for the Newton algorithm to converge. If the number of Newton iterations of the previous time step exceeds eight, the time step is decreased to half of its previous value. If the number of Newton iterations falls below four, the time step is doubled for the next time step. As mentioned above, explicit time integration schemes may only be used in a staggered solution strategy as a means to solve the equation of motion. The most popular explicit time integration method in structural mechanics is the central difference schemė which when substituted in (40) 1 , evaluated at time t n , results in the explicit relation for the current unknowns u n+1 , i.e.
In the first time step the quantity u −1 has to be initialized as Since equation (47) is linear in u n+1 it can be solved by inverting the matrix M uu . This is trivial if this matrix has a diagonal structure, which can be achieved by so called lumping strategies. The global mass matrix is assembled from all element mass matrices, i.e.
where n e is the total number of elements. For four node quadrilateral finite elements the nodal mass matrices have the structure as described in [10,22]. In this work, the lumped mass matrix is used whenever the central difference method is employed. The central difference scheme is not unconditionally stable. As mentioned in [23] and [22], the critical time step for the second order problem (40) 1 is governed by where c d is the dilatational wave speed and h is the characteristic length of the smallest element. The factor δ further reduces the allowable time step and accounts for nonlinear effects. The quantity ∆t cf l = h c d which is determined by the physics of the problem and the spatial discretization may be referred to as the Courant-Friedrichs-Lewy (CFL) time step in analogy to the convergence condition described in [15]. Thus, the central difference requires small time steps for stability reasons. However the staggered approach requires small time steps as well in order to achieve sufficient accuracy. The amplitude error of the central difference scheme is also of second order, see [22].
In a staggered approach, the decoupled phase field equation reduces to a nonlinear set of algebraic equations i.e. time integration is not necessary and the Newton-Raphson algorithm, see Fig. 4, is used to solve P s (u k+1 n+1 , s k+1 n+1 ) = 0 (57) for s k+1 n+1 , i.e. to perform step 2 of the staggered algorithm displayed in Fig. 2. The first benchmark problem considers a fractured 1D-bar with an initial crack as depicted in Fig. 5. The bar is constrained at its left end and subjected to a pressure pulse at the right end. For this problem we assume 1D-conditions, i.e. the displacements have the form u = u(x, t)e x , where e x is the unit vector in x-direction. With λ = 0 the equation of motion reduces to which resembles the 1D-wave equation, with wave speed c = √ E ρ . The 1D phase field equation (24) becomes The specimen is subjected to a sinusoidal pressure pulse t * (t) = σ 0 Y (t)e x at its right end. The amplitude of the pulse is σ 0 = −1.0

Physical Modeling for Virtual Manufacturing Systems and Processes Applied Mechanics and Materials Vol. 869 39
where t e = 0.25 L c . The absolute value of the stress pulse is distinctly smaller than the critical stress such that no further fracture is expected even after reflection of the elastic wave at the boundaries, see [14] for details on the critical stress σ c . The pressure pulse is initially travelling through the bar, i.e. with ξ = x − L 2 + ct the analytic solution is thus The bar is discretized by 200 homogeneously distributed elements with linear shape functions, i.e. the element size is h = 0.005L. The length scale parameter is ϵ = 0.01L. The initial crack is modelled by defining respective initial conditions s(x = L 2 , 0) = 0 for both nodes of the center element. In this example, no automatic time step control is used for the monolithic schemes (S1, S2, S3) and only one staggered iteration is performed for the staggered schemes (S4, S5, S6). Time steps are chosen to be 1∆t cf l for all schemes except for S6, the staggered strategy with explicit time integration, for which the time step is 0.5∆t cf l in order to ensure stability.
This numerical experiment is suited to judge the numerical dissipation involved in the solution strategies as well as their capability to accurately model the wave speed, because an analytic solution is available for comparison. Furthermore, because of the subcritical load, the focus is on analyzing the efficiency of the various schemes in handling the equation of motion (58). Since the time step is relatively small and the phase field is not expected to change very much, the restriction to only one staggered iteration (S4, S5) and the lack of the automatic time step control (S1, S2, S3) are regarded as acceptable. Of course, larger time steps can in principle be used for all solution strategies with implicit time integration of the equation of motion (S1-S5), which will be considered in the second numerical benchmark. Fig. 6 shows snap shots of σ xx at different times. Intitially, the numerical solutions all agree reasonably well with the analytic solution except for S3 for which the large amount of numerical dissipation of the Euler scheme becomes obvious, see Fig. 6 a). Furthermore, the staggered solution strategies S4, S5 yield virtually identical stress distributions compared to their monolithic counterparts at all times. The phase field does not evolve at all for all solution strategies. This is to be expected because of the subcritical loading and is the reason why S1, S4 and S2, S5 are not distinguished and only one graph of the phase field s(x, t) is shown in the plots. One can also observe a slight amount of dissipation in the Newmark-based strategies which manifests in a smeared out shape of the stress pulse compared to the analytic solution. This effect is larger in S2/S5 than in S1/S4. Furthermore, first oscillations appear for S2/S5 and S1/S4 where again the effect is larger for S2/S5. The stress distribution computed by S6 matches the analytic solution almost perfectly.
Just before the reflection of the initial pulse at the left boundary further differences become apparent, see Fig. 6 b). While the effects of numerical dissipation increase for S3, the strategies S2/S5 also deviate significantly from the analytic solution. The amplitude of the stress pulse decreases notably, but also the peak position does not match the analytic solution anymore. Thus, the wave speed is also not computed correctly. This is not the case for the standard Newmark scheme S1/S4. Caused by the initial oscillations which result in a partly positive stress, the stress waves are not completely transmitted by the crack but the tensile parts are reflected. Again, S6 almost yields the exact analytical stress distribution. Subsequently, the pulse is reflected at the fixed left end and returns in positive x-direction as a compressive wave. Consequently, there is again only minor interaction of the pulse with the crack in this stage of the simulation. Additionally, the crack field s is displayed as the dash-dotted line.
After reflection at the right end of the bar, see Fig. 6 c), the differences between the various scheme become more obvious, but show the same trends. Now, the main pulse is positive and thus will be reflected at the crack. Up to this point, only S2/S5 show significant oscillations. The strategies S6 and S1/S4 show no or minor oscillations respectively, as well as no notable decrease of the amplitude.
Significant oscillations can be observed for all strategies except S3 after reflection of the tensile pulse at the crack, see Fig. 6 d). Furthermore, the shape of the original pulse is not retained. Fig. 7 shows the runtimes for a simulation of the fractured bar in the time interval 0 ≤ t ≤ 4.0 L c obtained on a laptop computer with a 2.3 GHz Intel Core i7 processor. The monolithic schemes S1/S2 perform distinctly better than their staggered counterparts S4/S5. Thus, the algorithmic decoupling of the equation of motion and the phase field evolution equation does not hold any advantages for the implicit schemes in this example. This may be explained by the fact that, due to the subcritical load, the solution of the phase field equation is stationary and therefore a strong coupling of the phase field equation to the displacements does not exist anyway. However, the staggered algorithm performs significantly better if explicit time integration of the equation of motion is employed (S6). Although twice the number of time steps is required for S6, it outperforms all other strategies except the monolithic Euler strategy S3. In light of the short comings of S3 mentioned above and the almost exact match of S6 with the analytic solution in the first stages of the simulation, the performance of S6, which has to be attributed to the efficiency of the central difference scheme in solving the equation of motion, is even more impressive.
The second benchmark problem considers a tension loaded specimen with an initial crack as displayed in Fig. 9 a). Unlike in the first problem, the load is large enough to cause the crack to grow. The final crack pattern computed by S6 with ∆t = 0.95∆t cf l is displayed in Fig. 9 b). This benchmark allows to evaluate the performance of the different solution strategies in a more complex where v 0 = 1.0 c L and c = √ 2µ ρ , see also Fig. 8. The body is spatially discretized by a regular mesh of 200×100 quadrilateral finite elements with bilinear shape functions, i.e. the element edge length is h = 0.01L. The length-scale parameter is chosen to be ϵ = 0.02L. In this numerical benchmark, we also make use of the automatic time step control mentioned above (S1, S2, S3) and employ the staggered convergence criterion (38) for S4 and S5. The maximum number of staggered iterations is set to k max = 200. The maximum time step for the monolithic solution strategies (S1, S2, S3) is varied from 1∆t cf l to 40∆t cf l but can of course be decreased by the autmatic time step control. For the implicit staggered strategies (S4, S5) the time step is varied from 1∆t cf l to 40∆t cf l whereas for  Fig. 9: a) Setup for the crack branching benchmark: specimen with an initial crack, subjected to a tensile stress and b) plot of the phase field at t = 6.0∆t cf l computed by S6.
S6 we choose time steps in the range 0.2∆t cf l to 1.0∆t cf l . Convergence with respect to the time step size is obtained when the solution is insensitive to a further decrease of the time step. Fig. 10 shows the elastic energy E e , see (8), computed by the different solution strategies. The scheme S1 fails to converge for times larger than t ≈ 2.7 − 3.5[L/c] except for the smallest time step limit 1∆t cf l , see Fig. 10 a). The elastic energy increases significantly for the computations that do not converge. Thus, only a maximum time step of 1∆t cf l produces a meaningful result if S1 is used. The monolithic solution strategy with a modified Newmark scheme (S2) is less sensitive to the choice of the upper bound of the time step, see Fig. 10 b). All time step limits yield useful results. There are minor differences however such that for our purposes we consider simulations with a time step limit of 5∆t cf l or smaller as converged. The remaining monolithic solution strategy S3 only yields a converged solution for the smallest time step limit due to the severe numerical dissipation of the Euler scheme. In contrast to S1 however, the simulations quantitatively and qualitatively are in reasonable agreement with the converged solution, see Fig. 10 c).
The staggered strategies with implicit time integration S4 and S5 show similar behaviour as their monolithic counterparts, see Fig. 10 d) and e). The strategy S4 yields reasonable results only for the smallest time step. Again, the strategy with the modified Newmark scheme for time integration S5 allows to use larger time steps. Here, useful results are obtained for time steps that are smaller than or equal to 5∆t cf l . The staggered scheme with explicit time integration of the equation of motion S6 yields converged results for all chosen time step sizes except for 1.0∆t cf l , see Fig. 10 f). Thus, the maximum admissible time step for S6 is set to be 0.95∆t cf l . Fig. 11 displays the actual time step chosen by the automatic time step control of the three monolithic solution strategies S1-S3. Initially the time step is kept constant at the prescribed maximum admissible time step size. Just as the crack starts to grow the time step size is reduced step-wise and significantly in all schemes. An exception to this are the two smallest time step limits 1∆t cf l and 2∆t cf l for which the time step is kept constant throughout the whole simulation. Eventually, none of the strategies allows for an actual time step larger than 5∆t cf l in the later stages of the numerical experiment. Recall, that we considered the simulation with a time step limit of 10∆t cf l to be converged for S2. Interestingly, a final time step size of ≈ 3∆t cf l seems to be sufficient for S2, see the yellow line in Fig. 11 b). However, if S1 is employed, an even smaller time step 2∆t cf l fails to produce a converged solution, see the red lines in Fig. 11 a) and Fig. 10 a). Fig. 12 displays the number of iterations needed for the different strategies to converge for S4 and S5. Since the maximum number of staggered iterations is set to k max = 200, the staggered algorithm is assumed to have not converged for all computations reaching this limit. For S4, this is the case for the largest two time steps 20∆t cf l and 40∆t cf l , see Fig. 12 a). Interestingly, the computation with  the largest time step does not fail to converge for S5, see Fig. 12 b), but the simulations with 10∆t cf l and 20∆t cf l do. However, for the time step sizes 1∆t cf l (S4) and 5∆t cf l (S5) which are considered to be converged in the sense of an insensitivity of the solution to a further decrease of the time step size, the staggered algorithm also converges. The average, maximum and total number of iterations of the staggered algorithm for these two cases are summarized in Table 1. S5 needs a larger number of iterations in order to converge but allows for a larger time step as well. Thus, the total number of iterations is significantly smaller for S5. Fig. 13 displays the the obtained runtime of the different strategies versus the prescribed time step (S4, S5, S6), or time step limit (S1, S2, S3), respectively. The runtime of the largest time step that yields a converged solution in the sense that the solution is insensitive to a further decrease of the time step is marked by a black circle. As in the previous benchmark, the performance of S6 is remarkable. Although the time step is much smaller than for the other strategies, the total runtime is around half the value of the next best performing scheme, i.e. S2.
Further insight can be gained by the quantitative analysis presented in Table 2 which shows the average and maximum number of residual and tangent computations per time step as well as their total number. These operations together with the solution of the linear systems of equations that appear in every Newton iteration -their number is the same as the number of tangent computations -determine the computational effort and thus the runtime. It has to be kept in mind however that the residual computations, tangent computations and solving operations for the staggered algorithms are performed for the relatively small decoupled subproblems whereas the same operations in the monolithic scheme are performed for the full coupled problem and consequently are more computational expensive. Thus, the number of operations is naturally larger for the staggered algorithm, but this does not necessarily imply that the computational effort is also larger. Considering only the monolithic schemes S1-S3, we find that S2 performs best, because of the larger time step that can be chosen. S2 is only outperformed by the staggered strategy with explicit time integration S6. Interestingly, S6 requires more tangent, solves and residual computations than S2. However, as mentioned above, these operations are performed for the decoupled evolution equation (57), which reduces the computational effort significantly compared to the respective operations for the fully coupled problem (44).

Conclusions
This work discusses the efficiency of six strategies for the numerical solution of the coupled system of partial differential equations that arise from a phase field description of dynamic fracture. In the phase field fracture model, the dynamic fracture process can be described by two coupled partial differential equations. These two equations govern the evolution of the two primary fields which are

116
Physical Modeling for Virtual Manufacturing Systems and Processes Applied Mechanics and Materials Vol. 869 Fig. 13: Runtimes of the different schemes for the second benchmark versus the prescribed time step limit (S1, S2, S3) and time step (S4, S5, S6) respectively. the displacement and the phase field. The equation of motion primarily describes the evolution of the displacements but is coupled to the phase field via a degradation function. A further nonlinearity of the equation of motion is caused by the tension-compression split. The evolution equation on the other hand primarily describes the evolution of the phase field but is coupled to the displacements via the strain energy density that acts as a driving force. Spatial discretization of the problem by standard finite elements with bilinear shape functions yields a system of coupled ordinary differential equations. In order to solve this system two main strategies exist in the literature. Monolithic strategies solve the problem simultaneously for the two unknown fields at each discrete time step. In this case, only an implicit time integration scheme is viable due to the elliptic character and the lack of a time derivative in the evolution equation. In principle, these strategies allow for larger time steps that exceed the Courant-Friedrich-Lewy limit ∆t cf l dictated by the mesh and the wave speed of the problem. So-called staggered solution strategies on the other hand decouple the phase field fracture problem by freezing one of the primary fields and solving one of the two equations for the other field. Here, several advantages are obtained. First, by freezing one field some of the effects that cause nonlinearity of the field equations are removed in the subproblems. Second, the problem size for the decoupled subproblems is decreased. Third, flexibility is gained in terms of the available time integration methods for the equation of motion. An efficient explicit method can be used that also minimizes the computational effort that comes from dealing with the remaining nonlinear effects, i.e. the tension-compression split. However, the accuracy of these staggered strategies is critically dependent on the chosen time step. For explicit time integration this is not a problem at all since stability requires a very small time step anyway and thus only one staggered iteration, i.e. the subsequent solution of the equation of motion and the phase field evolution equation, is sufficient to guarantee accurate results. If implicit schemes are used for the time integration of the equation of motion, it is in principle possible to choose larger time steps. In this case, it is not sufficient Applied Mechanics and Materials Vol. 869 117 46 Physical Modeling for Virtual Manufacturing Systems and Processes to perform only one staggered iteration but the process has to be repeated several times. In order to judge the convergence of the staggered algorithm, a convergence criterion is proposed based on the contribution of internal 'forces' and accelerations to the global residual.
In this contribution we study three monolithic strategies: time integration by the standard Newmark method β = 0.25 (S1), time integration by a modified Newmark method β = 0.5 (S2) and time integration by the implicit Euler method (S3). These strategies adopt an automatic time step control that adjusts the time step with regard to the Newton iterations needed in the previous time step. Furthermore, three staggered strategies are discussed as well: time integration by the standard Newmark method β = 0.25 (S4), time integration by the standard Newmark method β = 0.5 (S5) and time integration by the explicit central difference method (S6).
The first benchmark problem is a 1D-bar with a stationary crack that is subjected to pulse load. As expected, strong numerical dissipation of the Euler scheme (S3) can be observed. Furthermore, the monolithic strategies (S1, S2) yield virtually identical results compared to their staggered counterparts (S4, S5). This is due to the fact that the phase field is stationary in this example due to the subcritical load. The modified Newmark schemes (S2, S5) also reveal some shortcomings. Their numerical dissipation is less severe than that observed for the Euler scheme (S3), but is notable. S1, S4 as well as S6 qualitatively yield the best results, when compared to the exact analytical solution. From a computational effort point of view, S6 shows the second shortest runtime and is only outperformed by S3 which struggles with the numerical dissipation problems mentioned above.
The second benchmark problem considers a tension loaded 2D specimen with an initial crack. This time the load is large enough to actually cause crack propagation. For this problem, a time convergence study is performed, i.e. simulations are run for varying upper time step limits (S1-S3) and time steps (S4-S6). As a measure to judge insensitivity of the computational results to a further decrease in the time step, graphs of the total elastic energy are recorded. The strategies S1, S3 and S4 only yield useful results for a time step size of 1∆t cf l where especially S1 and S4 yield useless results for larger time steps. The strategies with a modified Newmark scheme S2, S5 however allow for larger time steps, i.e. 5∆t cf l . The explicit scheme S6 requires smaller time steps of around 0.95∆t cf l . However, in spite of this small time step size S6 shows the smallest runtime, with S2 and S5 being the second and third fastest strategies, respectively.
It is concluded that the monolithic strategy with a modified Newmark scheme for time integration is a viable strategy to solve the dynamic fracture phase field problem. The other two monolithic strategies are found to be less useful because of the small admissible time step size (S1) or the significant numerical dissipation (S3). If a staggered approach is considered, explicit time integration of the equation of motion should be preferred because of the efficiency of this approach. Since a small time step (∼ 0.95∆t cf l ) has to be used for stability reasons, it is not necessary to perform more than one staggered iteration in order to obtain accurate solutions. Staggered strategies with the studied implicit time integration schemes on the other hand seem to hold no advantages over the explicit time integration strategy. The admissible larger time steps are more than compensated by the necessary iterations of the staggered algorithm.