A detailed investigation of the model influencing parameters of the phase‐field fracture approach

Phase‐field approaches to fracture are gaining popularity to compute a priori unknown crack paths. In this work the sensitivity of such phase‐field approaches with respect to its model specific parameters, that is, the critical length of regularization, the degradation function and the mobility, is investigated. The susceptibility of the computed cracks to the setting of these parameters is studied for problems of linear and finite elasticity. Furthermore, the convergence properties of different solution strategies are analyzed. Monolithic and staggered solution schemes for the solution of the arising nonlinear discrete systems are studied in detail. To conclude, we demonstrate the versatility of the phase‐field fracture approach in a real‐world problem by comparing different simulations of conchoidal fracture using structured and unstructured meshes.

model depends on the choice of a length-scale parameter, a mobility parameter and a degradation function. The length-scale parameter is a measure of the width of the crack and acts as a 2-fold parameter affecting the mesh and the material properties. The mobility is introduced as a regularization of the evolution in time, especially in a quasi-static setting. Besides the length-scale parameter whose influence has already been examined in the last years, [8,9] and the mobility parameter, we will focus on the choice of the degradation function. By means of the degradation function the loss of stiffness of the material in the cracked zone is modeled. While several conditions have to be fulfilled, the degradation function is more or less arbitrary. Commonly, in the standard phase-field model a simple quadratic degradation function is applied. However, there exist also different and more complex approaches which apply cubic degradation functions or degradation functions of exponential type. [25,26] The choice of the phase-field parameters influences also the numerical properties, such as nonlinearity and/or ill-conditioning. In the context of phase-field fracture models, the nonlinear coupled systems are usually solved by using either a monolithic [10] or a staggered solution scheme. [8] In this work, we study the convergence properties of nonlinear solution strategies employed in both solution schemes. A convergence study is performed with respect to the series of the phase-field parameters. A sensitivity analysis with respect to length-scale parameter can be found. [27] We expand this analysis by considering different choices of the degradation function and the mobility parameter.
This paper is structured as follows: after introducing the general equations of the phase-field method for both, linear and finite elastic problems, in Section 2, we will consider the influence of various parameters in Section 3. Exemplarily, a simple mode-I tension test is examined with focus on the length-scale parameter , the mobility parameter and the choice of degradation function . In Section 4, we investigate the robustness of solution strategies with respect to the influencing parameters. Finally, we present numerical simulations of a conchoidal fracture on two different geometries.

PHASE-FIELD FRACTURE MODELING
The phase-field model has been proven as a reliable and convenient approach to predict crack initiation and crack propagation.
In this section we present the main equations of the phase-field approach for brittle fracture. Starting with a domain ℬ ⊂ R 3 we consider its deformation within a time interval ∈ [0, ] ⊂ R + with the total time . The crack growth corresponds to creation of new surfaces Γ( ) of a priori unknown position. Thus, the total potential energy is composed of two main parts-one is the bulk energy of the solid with an elastic free Helmholtz energy density Ψ and the other part summarizes the fracture energy contributions which gives in total: In the case of brittle fracture the specific fracture energy can be interpreted as the Griffith's critical energy release rate. The challenging point for solving the surface integral directly is that the new surfaces Γ( ) are unknown a priori. While also classical discretization techniques do not lead to a satisfying solution, an additional field ( , ) ∈ [0, 1] is introduced. The limits of the phase-field parameter characterize the state of the material. In particular, the fractured state is given by = 1 and the solid material is indicated by = 0. The states ∈ (0, 1) do not have any physical meaning. The phase-field is a continuous field and the crack boundaries are 'smeared' over a finite zone of width . Because of these characteristics the phase-field model can be interpreted as a diffuse interface approach. By introducing a surface density function ( ) which depends on the phase-field parameter , the surface integral can be approximated then by Consequently, the total potential energy can be formulated by inserting (2) into (1) such that the optimization problem with respect to the primal variables { , } results locally in There exists no unique choice of the crack density function and so a series of different possibilities have been suggested. [8,9,22,23,28] However, in general a second order approach [8] is applied and is given by F I G U R E 1 Phase-field function in an analytical solution of the second order model = exp The phase-field within an analytical solution of the second-order approach is depicted in Figure 1. Inserting the second-order approach (4) in Equation (3) the total potential energy results in The temporal evolution of the phase-field can be formulated according to Allen-Cahn as a nonconserved quantitẏ where the mobility parameter is the inverse of a viscosity and ( , ) summarizes all possible driving forces. Commonly a dimensionless formulation is meaningful, such that the evolution equation is reformulated aṡ with a retardation time (seconds) and a normalized driving forcē(-). In variational form it readṡ where the normalized energy is given byΨ = ∕ Ψ and is based on a coupling of the phase-field and the mechanical field, and requires an appropriate tension-compression decomposition. The dimensionless driving forcēcan be decomposed into the crack driving forcēitself and a crack resistancē, such that̄=̄−̄holds and the evolution equation is given bẏ In the standard formulation the crack driving force is based on an energy minimization which corresponds for brittle fracture to Griffith's classical model. The crack driving force can also be motivated by already established failure criteria in fracture mechanics. [24] In this case the crack driving forces do not minimize the energy potential. In this work, we focus on the maximum principal stress criterion, namely a crack starts propagating as soon as a critical stress is exceeded. This criterion can be traced back to the Rankine failure criterion, therefore, we denote this approach in the following as Rankine model or stress based formulation. To state an expression of the crack driving force the principal stresses have to be calculated. As opposed to the classical phase-field model the principal stresses are not degraded for determining the crack driving force. Consequently, it is given bȳ= ⟨ with the maximum stress = max( ) with the eigenvalues , ∈ {1, 2, 3} of the stress tensor . For the sake of completeness we mention that in this work the irreversibility of crack propagation is considered by fixing the phase-field at = 1 as soon as the material cracks. [9]

Linear elasticity
Let us start with the main equations of the linear theory, where the strain energy function is formulated as by using the strain ( ) and the elasticity tensor C. In general, it is necessary to distinguish between tension and compression states in a way that only tension induces crack propagation. Thus, to take only the tensile parts into account the energy density is split up into positive and negative components. A degradation function ( ) ∶ [0, 1] → [ , ∞) reduces the cracked material's stiffness by modifying the tensile part of the energy density according to Here the positive and negative parts of the energy are given by Ψ ± = 1 2 ⟨tr ⟩ 2 ± + ± ∶ ± , whereby the positive and negative parts of the strain tensor are based on the spectral decomposition ± = ∑ 3 =1 ⟨ ⟩ ± ⊗ . This decomposition accounts for a degradation of tensile parts and allows contact states of the compressive parts. The exact choice of the degradation function will be analyzed later in Section 3.2. For reasons of simplicity, we denote this decomposition in linear elasticity as -split.

Finite elasticity
By now the phase-field model has been introduced for linear kinematics, but it is also valid for finite deformations. Therefore, the required notations will be introduced in this subsection. The deformation mapping is given as The fields in capital letters refer to the initial configuration, whereas the fields in small letters are used for the current configuration.
The deformation gradient is defined as = ∇ and its determinant as = det . The right Cauchy-Green stress tensor and its isochoric part can be determined as = and = 2∕3̄. Commonly, the elastic strain energy function is formulated with the well-known principal invariants of the right Cauchy-Green deformation tensor for an isotropic material. Namely, these invariants are given by 1 = tr( ), 2 = tr(cof ) and 3 = det . For numerical computations of (nearly) incompressible materials a multiplicative decomposition of the deformation gradient into a volume-changing and a volume-preserving part is used, that is, = 1∕3̄w herēcorresponds to the volume-preserving part of the deformation gradient. Exemplarily, we formulate the strain energy function of the Mooney-Rivlin material model with the material constants 1 , 2 , the bulk modulus and the invariants̄1 and̄2 of̄. This strain energy density is based on the isochoric split of the deformation gradient, where the second term of the strain energy function has to be adapted slightly such that the strain energy function fulfils the requirements of a polyconvex function. [29][30][31] More details about polyconvex terms and functions can be found in the literature. [29,32] Analogous to the linear case only the tensile parts induce the crack growth and the stiffness in the crack has to be degraded artificially. Thus, the strain energy function is decomposed such that it leads to the following form with a degradation function ( ). The positive and negative parts of the invariants used in Equation (14) are given bȳ The formulation of the crack-driving force is of variational type, = Ψ . Please note that the different decompositions of the strain energy function in positive and negative parts and the coupling with the phase-field are quite arbitrary. In finite deformations we can also apply ad hoc crack driving forces proposed by Bilgen and Weinberg [24] for Equation (9) which are not stated again here.

2.3
The formulation of the problem In this subsection we introduce the main equations of the coupled phase-field model. The numerical realization is based on the finite element method and the variational principle. Besides the evolution equation of the phase-field in Equation (6) the simulation solves the balance of linear momentum which is stated as div + = 0 (18) for finite deformations. The equations are formulated in the following for finite elasticity. In linear elasticity, the balance equation is similar to Equation (18), where the first Piola Kirchhoff stress tensor is substituted by the Cauchy stress . The weak formulation of the coupled system reads as follows: The symbol 1 in the above formulation denotes the Sobolev functional space of square integrable functions with square integrable derivatives. The functional ( , ; ) describing the balance of linear momentum is defined as where is a body force and represents a traction force. Analogous to Equation (20) the evolution equation of the phase-field can be also written in a weak form as Equation (19) can be solved by using either a monolithic [10] or a staggered solution scheme. [8] The monolithic scheme requires the solution of the nonlinear coupled problem on each time step, see Algorithm 1. This is challenging, as the coupling between the displacement and the phase-field results in a nonconvex minimization problem. In the staggered solution scheme, the functional (19) is split into two functionals related to the displacement and the phase-field. These two functionals are then minimized alternately until the system converges, see Algorithm 2. We terminate the staggered scheme once the change in the phase-field is smaller than some threshold . Alternative stopping criteria can be found. [33,34] Compared to the monolithic scheme, it is easier to design a reliable solution strategy for the decoupled subproblems arising in the staggered scheme. Hence, the staggered scheme became popular and preferred in practical applications. One disadvantage of it might be that it tends to underestimate the speed of crack evolution. [34] However, the quantification of the speed of crack growth is not meaningful in a quasi-static setting. In order to solve Equation (19) numerically, we discretize the weak form (20) and (21) Find ∈ 1 (ℬ 0 ) ∶ ( , −1 ; ) = 0, ∀ ∈ 0 . 6: Find ∈ 1 (ℬ 0 ) ∶ ( , ; ) = 0, ∀ ∈ 0 . { +1 0 , +1 0 } = { , } 10: end for nonoverlapping elements such that Then both fields and are approximated with the same ansatz functions as wherê( ) and̂( ) are the nodal values of the displacement and the phase-field, respectively. Inserting the approximations (23) and (24) into the weak forms in Equations (20) and (21) leads to the finite element system of equations. By dividing the time interval [0, ] into pairwise disjoint equidistant subintervals such that the time step is given by Δ ∶= +1 − , the performed simulations are based on an Euler-backward scheme. More details about the solution process will be presented in Section 4.

INFLUENCE OF VARIOUS PARAMETERS
Applying the phase-field model several parameters influence the results, in particular, the crack initiation and patterns. In the following we focus on the length-scale parameter , the choice of degradation function and the mobility parameter , respectively, all within the framework proposed in Section 2. The simulations of this study are based on a staggered solution scheme in linear elasticity if not mentioned otherwise.

Length-scale parameter
Let us start examining the regularization length which represents the width of the transition from = 0 to = 1. The length-scale parameter acts as a 2-fold parameter. On the one hand, the length-scale parameter has to be adapted to the given mesh size to be able to resolve the transition. Thus, it can be calculated via the mesh size ℎ whereby the relation > ℎ has to hold true. It is recommended to choose the length-scale parameter as = 1.5, … , 3ℎ depending on the approximation order of the ansatz functions. [8] On the other hand, it serves also as material parameter, because it scales the critical energy release rate which is a critical value for crack propagation. As soon as the elastic energy of the material reaches the fracture energy the crack starts growing.
A simple mode-I tension test is performed in two dimensions with the mesh and the related boundary conditions presented in Figure 2. We consider a (100 × 100) mm 2 squared plate with a prescribed notch in the meshing and apply a refined mesh which consists of 2656 quadratic B-spline elements after three local refinement levels. A detailed derivation of the hierarchical refinement strategy can be found. [21] The lower boundary is clamped in -and -direction, whereas the upper boundary is pulled upwards by an incrementally prescribed displacement of Δ = 10 −8 m. Also both boundaries on the left and the right side are restricted in -direction. The material parameters are given by the Young's modulus = 50 000 N/mm 2 , the Poisson's ratio F I G U R E 2 Schematic boundary conditions used for all calculation of Section 3 (left), the mesh in two dimensions (middle), and the phase-field indicating the crack path (right)

F I G U R E 3 Load-displacement curves with varied length-scale parameters in the variational formulation
F I G U R E 4 Phase-field snapshots for different length-scale parameters to illustrate the dependence of the crack width on the length-scale parameter = 0.2 and the critical release rate = 75 × 10 −3 N/mm. If it is not mentioned otherwise a quasi-static setting is conducted, that is, the evolution equation simplifies to 0 =̄such that the influence of the mobility is omitted. The related crack path is presented in the right-hand side of Figure 2. In order to investigate the influence of the length-scale parameter the tension test is simulated by varying the length scale parameter . For the standard phase-field method there exist a series of investigations. [8,9] In Figure 3 the load-deflection curves are depicted for different length-scale parameters = {1.25, 2.5, 5, 10} mm. It is observed that the displacement at crack initiation does not change within the variation; however, the required force for cracking increases with smaller length-scale parameters. The curves resulting from small length-scale parameters show an abrupt cracking. This trend coincides with the theory that for → 0 the results tend to a sharp crack interface. The load-deflection curves for larger length-scale parameters might indicate that the solution has not been converged in time yet. Please note that the other material parameters are unaltered and different parameter combinations or using another solution scheme might lead to different results. In Figure 4 the crack path itself is shown and the influence of the length-scale parameter with respect to the width of the diffuse interface zone can be observed. Thus, the results are improved for smaller length-scales which corresponds to the fact, that for → 0 the diffuse interface converges to a sharp crack.
In the following, the influence of the length-scale parameter is regarded in the context of the Rankine model by means of the proposed mode-I tension test, see Figure 2. Also in this case the length-scale parameter acts as a material and a regularization which already demonstrates the dependence on the length-scale parameter. The degraded critical stress results in = 9 16 √ 6 . [9,28,35] It has to be remarked that the derivation of this formula is based on the simple assumptions of linear elasticity and a quadratic degradation function and might be expanded for more complex problems. The critical stress influences the crack driving force and is also relevant after crack nucleation. In the left plot of Figure 5 the load-deflection curves are depicted for ∈ {1.25, 2.5, 5, 10} mm. By analyzing these curves it can be concluded that for larger length-scale parameters the load and the prescribed displacement at crack initiation increases. Furthermore, the shape of the curve becomes spikier for smaller length-scale parameters which indicates the property of brittle fracture.
If the critical stress is adapted while the length-scale parameter is varied, the load-deflection curves differ only slightly, see right plot in Figure 5, that is, the force at crack initiation decreases for bigger length-scale parameters. Also in this case the progression of the curve characterizes the brittle behavior better for smaller length-scale parameters. Considering Equation (25) it is obvious that the critical stress increases with smaller length-scale parameters. If the critical stress is not adapted the critical stress does not decrease for bigger length-scale parameters, in particular the driving force remains constant such that the crack starts propagating later. If we vary the critical stress with respect to the different length-scale parameters, this effect is compensated. Applying the adapted critical stress it is observed that the required force for cracking increases for smaller length-scale parameters which is similar to the variational approach. The only difference is that is influenced by √ 3 in the stress-based formulation as opposed to the variational ansatz by . Thus, the critical stress depends to a great deal on the length-scale parameter and by association on the choice of mesh such that the critical stress has to be adjusted with the mesh and the length-scale parameter.
The different effects of varying the length-scale parameter within the variational approach denoted by and within the stress-based approach with constant ( = 31.6228 N/mm 2 ) and adapted critical stresses ( ( )) are depicted in Figure 6. The maximum load decreases with rising length-scale parameters in the Rankine model with the adapted critical stress and the variational approach. As opposed to that by varying the length-scale parameter and constant critical stress the maximum load rises with increasing values of . Analogous the cracked zone depending on the length-scale parameter demonstrates the same behavior. Finally, it can be concluded that the length-scale parameter should be chosen quite small, however, this choice is restricted by numerical possibilities and the relation between the material parameters respectively and the length-scale parameter has to be accounted for.

Degradation function
In the phase-field model it has to be taken into account that the stiffness of the material vanishes in the cracked zone and that the material cannot sustain tension. To map this effect the material has to be degraded. The decompositions (12) and (14) account for the degradation of the material during cracking by means of a degradation function ( ). For = 1, that is, in the crack, a physically motivated degradation has the form (26) and is responsible for the loss of stiffness in tension but contact conditions in compression. A general contact energy is a function of normal and tangential tractions, Ψ contact ( , ). Here, the energy contribution Ψ contact equals the compressive part of the strain energy Ψ − which leads to the introduced formulation in Section 2.2. However, these tensile/compressive decompositions are quite arbitrary. Please note that the crack driving force in (10) is not influenced by the degradation function anymore.
In this subsection the influence and the choice of degradation function are considered in more detail. The crack initialization is affected as soon as the phase-field parameter increases the first time. [36] The formulations and statements in the previous sections were based on a quadratic degradation function first of all. Now we analyze a polynomial ansatz of degree three [37] : There exist a few conditions which have to be fulfilled by the degradation function such that the formulation is meaningful: Furthermore, the following condition can be also considered by introducing an additional parameter such that In 2014 Borden and Hughes proposed this condition and stated that for very small parameters, that is, → 0, the solution of the governing equations is quite uniform for any state of strain. [37] By making use of these four conditions (28) and (29) the constants , , and can be calculated. Thus, the cubic degradation function results as follows: Because of the physical meaning of the phase-field parameter the degradation function has to be a positive and (monotonically) decreasing function for ∈ [0, 1]. Thus, it can be summarized that for the additional parameter ∈ [0, 2] holds true. In Figure 7 the cubic degradation function (30) is presented for the different choices of ∈ {0.01, 0.1, 0.5, 1, 2} whereby for the choice of = 2 the quadratic degradation function is given. It becomes clear that for a smaller parameter the curve is flatter for = 0.
Conducting the mode-I tension test with the related boundary conditions and material parameters as given in the previous section, see Figure 2, we investigate the influence of varying the degradation function by changing the additional parameter . Both approaches, the standard variational ansatz and the stress based formulation, are considered. In Figure 8 the load-deflection curves are illustrated for different choices of parameter . It can be observed that for smaller numerical values of the crack initialization starts within a slightly larger displacement and a bigger force. This coincides for both models. In general, using the ad hoc crack driving forces the property of brittle fracture can be simulated very well. After crack initiation the force decreases abruptly. For the variational approach this behavior can be studied just for smaller numerical values of . In addition to that it has to be mentioned that the contribution of relatively small strains far away from the crack are lower for smaller numerical values of , see Figure 9. This confirms the desired modeling effect that prior to fracture the degradation function is nearly 1 for small values of , that is, hardly any degradation is applied. It corresponds with observations of Borden et al. [37,38] and Vignollet et al. [36] Regarding the change of the force during modifications of the degradation function it can be said that the parameter and the force depend on each other linearly for the Rankine model, see Figure 10. In the standard variational approach the influence of smaller becomes larger. The related values in percentage terms are summarized in Table 1 for the energy-based and the stress-based approach.

Mobility
In the following we focus on the mobility parameter which influences the phase-field evolution Equation (6). This parameter [m 2 /Ns] is introduced to regularize the quasi-static case and can be interpreted as an inverse kinematic viscosity. In order to evaluate the influence of the mobility we simulate the simple mode-I tension test in the same way as in the previous sections with the boundary conditions and material parameters given in Figure 2. Within this tension test the mobility parameter is varied in the range of ∈ {1, 10, 50, 100, 1000, 10 000} mm 2 /Ns for both crack driving forces-on the one hand the energy-based approach and on the other hand the stress-based model. Because the crack path itself does not change only the load-deflection curves are presented in Figure 11. It becomes clear that for increasing mobility the curves converge concerning the applied force and prescribed displacement, that is, for high numerical values of the mobility parameter the crack starts propagating more rapidly. In the opposite case, for low mobility values the crack growth slows down and the phase-field evolution is locally restricted. Thus, the fracture strength of the material is artificially increased. In order to consider a realistic crack growth the phase-field has to significantly evolve within one step of time discretization. This results in a rule of thumb for the size of mobility which depends on the order of magnitude of the time step assuming it to be increasing linearly with the applied load. We recommend to choose  This is true for the ad hoc crack driving forces and the standard variational approach in a similar way. Another way to circumvent the choice of the mobility parameter might be decreasing the time step such that the influence of the mobility parameter becomes negligibly small. [39]

SOLUTION PROCESS
Both solution schemes, monolithic and staggered, give rise to large scale systems on each time step. This is due to the fact, that in order to resolve smoothed fracture zones accurately, the mesh has to be sufficiently fine. As a reminder, the requirement ℎ ≤ has to be fulfilled with the mesh size ℎ and the length-scale parameter . In order to solve such a large scale system efficiently, multilevel solution strategies can be used. Multilevel methods provide optimal complexity per step by employing a hierarchy of usually nested finite element spaces, also called levels. The solution process then combines smoothing and the F I G U R E 11 Stress-displacement curves with varied mobility parameters for the variational approach (left) and the stress-based formulation (right) coarse grid correction step. Smoothing is used to reduce the high-frequency error related to each level, while the coarse grid correction step eliminates the low-frequency error remaining on the coarsest level. Linear geometric and algebraic multilevel methods have been developed for more than 50 years and they are used extensively in many areas. An introduction to linear multilevel methods can be found. [40,41] In the context of phase-field fracture problems, linear multilevel methods were applied. [27,42,43] Nonlinear multilevel methods, such as full-approximation scheme (FAS), [44] nonlinear multigrid, [45] optimization multigrid (MG-OPT), [46] or recursive multilevel trust region (RMTR) [47] provide the generalization to the linear variants. They enhance the convergence by computing coarse level corrections based on the nonlinear coarse level representation of the original/fine-level problem. Regarding the phase-field fracture problems, the nonlinear multilevel method was recently applied, [48] in which a significant speed-up strategy is demonstrated and compared with a single-level solution strategy.
In the forthcoming section, we discuss how to employ both linear and nonlinear variants of the multilevel method in order to solve phase-field fracture problems. Our discussion refers to both solution schemes, the monolithic and the staggered scheme.

Monolithic solution scheme
The monolithic scheme, Algorithm 1, requires the solution of a highly nonlinear, ill-conditioned, large-scale problem in each time step. A standard approach to deal with nonlinearities is Newton's method as it is simple and locally quadratically convergent. [49] The main drawback of the method is its global convergence which relies on a good initial guess. In order to ensure convergence to local minimizers, globalization methods such as trust region, [50] or line-search [49] can be applied. In this work, we employ a line-search globalization strategy specifically tailored for the phase-field fracture problems. [34] A single iteration of Newton's method is defined as where the iteration number is denoted by subscript and ∈ R + is the step-size. All quantities in (32) have the following block structure The solution vector ∈ R contains the coefficients { } and { } related to the displacement and phase-field, respectively. The blocks of the residual vector ∈ R consist of the terms corresponding to (20) and (21). The entries of the Jacobian matrix ∈ R × are obtained by computing the directional derivatives of (19).

Field-split preconditioner with linear multigrid
The main computational cost of Newton's iteration (32) is associated with computing the search direction , for example, with solving the linear system −1 . Here, we employ the MINRES Krylov solver, [51] as the Jacobian matrix is symmetric, but not necessarily positive definite, especially far from the solution. Since the Jacobian is also poorly conditioned due to strongly varying stiffness induced by spatially varying fracture localizations, we employ the preconditioner −1 to speed up convergence. We take the advantage of the block structure of and employ a preconditioner based on a Schur complement method. [52] Let the matrix = − ( ) −1 be the Schur complement of with respect to , see Benzi et al. [52] for details. We assume that the matrix is nonsingular, then the inverse of Jacobian ( ) −1 can be expressed explicitly as Computing the explicit inverse ( ) −1 is computationally expensive. The most tedious part is computation of the Schur complement , as it requires solving ( ) −1 . To simplify the solution process, we approximate with , thus ≈ . The preconditioner −1 obtained in this way has the same block structure as ( ) −1 described by (34), but we replace −1 with ( ) −1 . An application of this preconditioner −1 requires one application of ( ) −1 and two applications of ( ) −1 .
In our implementation, both actions of the inverse are approximated by two V-cycles of algebraic multigrid method. [53] In order to make the algebraic multigrid as efficient as possible, we supply near-nullspace vectors, for example, eigenvectors associated with the eigenvalues of small magnitude. [41] Since the block is similar to the standard elasticity operator, we provide the near-nullspace vectors related to the rigid body modes. Additionally, the operator has near-nullspace associated with the independent rigid body motions of the broken regions. [27] We do not provide these vectors, as they are computationally expensive to obtain. For the block , we supply near-nullspace specified by a constant vector.

Nonlinear multilevel method
Nonlinear multilevel methods provide a possibly computationally cheaper alternative to the standard globalization schemes for Newton's method. However, the majority of nonlinear multilevel methods is suitable only for convex minimization and is not well suited for the minimization of nonconvex functionals, such as our phase-field fracture energy (3). The only globally convergent nonlinear multilevel method for nonconvex minimization problems, called RMTR, was introduced by Gratton et al. [47] and further developed by Groß and Krause [54] and Gratton et al. [55] The RMTR method provides the global convergence by inserting the trust region globalization strategy [50] into the nonlinear multilevel framework. As other multilevel methods, the RMTR method alternates between nonlinear smoothing and coarse grid step. The nonlinear smoothing is used to reduce high-frequency error related to each level, while the coarse grid correction step aims at eliminating the low-frequency components of the error. The RMTR method performs both, smoothing and coarse-grid step, by minimizing nonlinear level dependent objective functions. The minimization on each level is carried out by the trust region strategy. Once the approximate minimization on the given level is finished, the obtained coarse level correction is brought back to the fine level. However, the transferred correction is not added immediately to the current iterate, but only if it provides a decrease in the original/fine level objective function. Rigorous details about how to apply the RMTR method to the phase-field fracture problems can be found. [56]

Staggered solution scheme
The staggered solution scheme requires the solution of two subproblems, see Algorithm 2. The first subproblem is associated with finding the updated displacement field . This requires solving a problem similar to standard elasticity, however, with strongly varying stiffness parameter induced by evolving fracture. Here, we can apply either the RMTR method, see Section 4.1.2 or the Newton's method with the line-search, see Section 4.1. The linear systems arising on each Newton's iteration have a similar structure as the block from (33). To solve these systems efficiently, we employ the conjugate-gradient (CG) method [57] preconditioned with the algebraic multigrid. The second subproblem consists of finding the updated phase-field . Structurally, the second subproblem is akin to a generalized Helmholtz problem, also with spatially varying coefficients. As the phase-field subproblem is linear, we employ CG method preconditioned with an algebraic multigrid. In both cases, the algebraic multigrid was configured as described in Section 4.1.

Effect of the phase-field parameters on the convergence of the solution strategies
We investigate how the model parameters influence the solution process. For benchmarking purposes, we consider a simple 2D tension test setup on a uniformly refined squared plate of size 100 mm × 100 mm with a prescribed notch as depicted in Figure 2, where we prescribe incremental displacements Δ = 5 × 10 −5 m on the upper boundary in -direction. The results presented here were obtained by using the energy-based approach with a linear elastic material model. The discretization was performed on uniform mesh with 153 600 nodes. We employed the Lagrange P1 finite elements. The material parameters were chosen as described in Section 3.1. If not specified otherwise, we employ the quadratic degradation function, the length-scale parameter as = 2ℎ, where ℎ is the mesh size. By choosing the mobility parameter → ∞ we can approximate the quasi-static limit and obtain a rate independent model. [58] This coincides with the simulations in the previous Section 3.
Despite the robustness and the efficiency of the RMTR method, its practical implementation is a technically demanding task, which requires a suitable finite element framework. As a result, in practice the method is not as prevalent as standard linear multigrid which is applied to the linear systems arising during Newton's iteration. For this reason, the parameter study performed in this section considers only Newton's method with line-search. However, we demonstrate the overall performance of the RMTR method in Section 5 for the selected example of the conchoidal fracture. We will investigate the convergence properties of the RMTR method with respect to the phase-field parameter and present a comparison with the standard Newton's method.
The linear and nonlinear solution strategies employed in this section were implemented inside of the UTOPIA library. [59] The implementation highly relies on the PETSc library [60] for both linear algebra and linear solvers. During all tests, we used following stopping criterion: ‖ ‖ < 10 −7 , where ‖ ‖ is Euclidean norm of the residual .

Robustness with respect to the length-scale parameter
As pointed out in Section 3.1, the choice of the length-scale parameter is limited by the mesh size ℎ. In particular it is recommended to choose the length-scale parameter as = 1.5ℎ, … , 3ℎ. The mesh dependence of causes difficulties while designing efficient multilevel strategies, for detailed discussion see Kopaničáková and Krause. [48] This is especially the case for geometric multilevel methods, where a hierarchy of nested finite element spaces is based on a hierarchy of meshes with different resolution. Here, we overcome this difficulty by employing the algebraic multilevel method, which assembles the coarse level operators by using information stemming from the quadrature of the fine level problem and corresponds to the so-called Galerkin assembly. [40] We investigate the effect of the length-scale parameter on the convergence of nonlinear as well as linear solution strategies. Table 2 summarizes the results obtained for the monolithic solution scheme. We clearly observe that the number of cumulative nonlinear iterations over all time steps increases as decreases. In order to evaluate the computational cost required per time step, we measure the average number of nonlinear iterations over all time steps. As it can be observed from Table 2, the average number of nonlinear iterations also decreases as increases. The same behavior is observed for the average number of preconditioned MINRES iterations. Here, the average is computed with respect to all Newton's iterations encountered during all time steps. In the end, we study how the mesh size ℎ influences the convergence properties of the solution strategies. For this reason, we consider the hierarchy of three meshes obtained by uniformly refining the initial mesh with 153 600 nodes and mesh size ℎ 0 = 0.0025 m. As illustrated in Table 2, the number of all iterations increases rapidly with the increased mesh resolution.
The same benchmarks as described above were performed for the staggered solution scheme. Firstly, we study the behavior of the elasticity subproblem, Table 3. We observe that the number of accumulated and average nonlinear iterations decreases with increased . In contrast to the monolithic scheme, the average number of linear iterations decreases with decreasing . Secondly, we focus on the phase-field subproblem, Table 4. Here, we study only the average number of the CG iterations over all time steps, as the underlying problem is linear. As depicted in Table 4, the average number of linear iterations decreases while the length-scale parameter increases. As for the monolithic scheme, the number of linear and nonlinear iterations increases rapidly with decreased mesh size. This pattern can be observed for both subproblems, that is, elasticity and the phase-field.

Robustness with respect to the choice of the degradation function
In this section, we study the numerical behavior of our solution strategies with respect to the degradation function ( ). The degradation function considered here is defined as in Section 3.2, where we are free to choose a parameter ∈ [0, 2]. The results obtained for the monolithic scheme are demonstrated in Table 5. Results for the staggered solution scheme can be found in Tables 6 and 7 for the elasticity and phase-field subproblems, respectively. The obtained results illustrate that the number of accumulated iterations over all time-steps increases as decreases. This is not very surprising, as for the smaller values of , the crack initialization requires larger displacement, see Section 3.1. Consequently, the number of time-steps also increases, given our displacement driven loading scenario.
Further, we are interested in the average number of nonlinear iterations over all time steps. For the monolithic scheme, we can see that the average number of nonlinear iterations increases rapidly for smaller values of . In contrast, the results obtained for the staggered solution scheme do not show a significant increase in the number of nonlinear iterations. These measurements suggest that the choice of the degradation function affects only the coupling relation between the displacement and the phase-field. This observation comes along with the known theory, which identifies the coupling term ( )Ψ + ( ) in (12) as nonconvex in both displacement and the phase-field simultaneously. [61,62] Especially, the derivative of (20) with respect to , block in (33), is known to create the most challenges while designing reliable and efficient solution strategies. [63] As a result, the choice of the degradation function ( ) influences only the nonlinearity of the coupled problem arising from the monolithic scheme. At the end, we investigate the influence of concerning the convergence properties of the linear solution strategies. For the monolithic case, Table 5, the average number of linear iterations slightly increases as decreases. For the elasticity subproblem resulting from the staggered solution scheme, Table 6, the average number of linear iterations seems to be unaffected by the choice of ( ). Table 7 demonstrates results obtained for the phase-field subproblem. We report only the average number of CG iterations over all time steps, as the problem at hand is linear. As we can observe, the average number of the CG iterations seems to be also unaffected by the choice of the degradation function.

Robustness with respect to the mobility parameter
In this section, we study how the mobility parameter influences the convergence properties of the presented solution strategies. We alter the mobility parameter in the range of ∈ {10, 50, 100, 1000, 10 000} m 2 ∕Ns. The results obtained for the monolithic scheme are depicted in Table 8. We observe that the cumulative number of nonlinear iterations increases as the mobility decreases. This effect is not surprising as the lower value of the mobility parameter speeds down the phase-field evolution. Similarly, the average number of nonlinear and linear iterations increases as decreases. This behavior was also expected, as the lower values of the mobility parameter artificially increase the stiffness of the material. For the mobility parameter = 10 000 m 2 ∕Ns, we obtain results more similar to the results gathered in the previous section, for example, where → ∞. In a next step, we analyze the staggered solution scheme. Tables 9 and 10 present the numerical results for the elasticity and the phase-field subproblem, respectively. Here, we observe the same behavior as for the monolithic scheme. The number of all measured iterations increases as the value of decreases.

NUMERICAL SIMULATION OF CONCHOIDAL FRACTURE
Conchoidal fracture is a special type of brittle fracture where the crack initializes inside of the body under the surface. Commonly it occurs in fine-grained or amorphous materials like rocks, minerals, or glasses. Thus, it acts as an appropriate example to show the applicability of the proposed models because the location of crack initialization constitutes a difficulty in numerical simulations, in particular, when there exists no initial notch, kerf, or precrack. The main characteristic of conchoidal fracture is that the shape of the crack surface looks like a mussel shell, that is, the loading leads to a curved breakage with a rippled surface. This rippling can be traced back to the work of Wallner [64] who observed "Wallner lines" in cracked specimen. It has to be remarked that the fractures do not follow any natural planes of separation but rather the solid material cracks arbitrarily by cleavage. This model has been proposed as a benchmark problem within the scope of the DFG Priority Program 1748. The conchoidal fracture example was investigated in detail with respect to the length-scale parameter. [42] In this contribution, we expand the examinations to the influence of the degradation function. In addition, we also study the effect of the degradation function and the length-scale parameter from the perspective of the nonlinear solution strategies. We present the conchoidal fracture examples for two geometries. Firstly, we consider a block of stone material and of size 4 × 4 × 2 where on the upper boundary a square surface is pulled upwards with an edge length of 2 = 1 m. The upper boundary condition is based on a displacement-driven deformation with prescribed incremental displacements of 0.001 mm. The lower boundary is constrained in all directions, that is, = 0. The whole setting is depicted in Figure 12. The finite element discretization is based on 30 × 30 × 30 eight-node brick elements such that the length-scale parameter of = 0.2 m is adapted. We do not use a refinement strategy in order to exclude any influences of the mesh. The material parameters are given by the Lamé constants = = 100 000 N/mm 2 and a critical energy release rate of = 1 N/mm. Second, we demonstrate the applicability of the phase-field method on the unstructured geometry. For this reason, we consider the same boundary value problem setup, but on the geometry of a real rock, see the right plot of Figure 12. The rock example uses the linear material model described in Section 2.1, while the brick example applies the nonlinear Neohooke material model, see Section 2.2. Considering the choice of crack driving forces, in both cases the typical shape of crack surface results and the rippling of the crack surface can be observed, see Figure 13. Furthermore, in both models the crack starts propagating inside of the body.
Next, we investigate the influence of the proposed parameters of the phase-field model within the conchoidal fracture example based on the structured geometry. Here, we refer to detailed statements regarding the length-scale parameter in the publication. [42] It can be summarized that the location of crack initialization depends heavily on the choice of the length-scale parameter which supports the results concerning the meaning of the length-scale parameter in the previous section. Moreover, we simulate the conchoidal fracture example within different degradation functions for both crack driving forces-on the one hand the standard variational approach and on the other hand the ad-hoc crack driving force. The related load-displacement curves are shown in Figure 14. It can be observed that-similar to the tension test in the previous section-by decreasing the parameter the applied force at crack initialization is increased. This effect occurs for both choices of crack driving force. Regarding the crack path itself, the related phase-field snapshots are depicted in Figure 15. In the energy based approach a further characteristic is that stresses far away from the crack tip influence the crack propagation for larger parameter . Thus, the crack path itself also changes within the standard model. The surface of the crack is more flat for the quadratic degradation function and becomes more curved for decreasing parameters , as it can be observed in Figure 15. Therefore, the results applying the cubic degradation function in the variational formulation indicate the characteristic of the rippled curved breakage of conchoidal fracture more pronounced as the results applying the stress-based approach. In contrast to that, this effect is not observed within the phase-field model based on the Rankine crack driving forces. Here, the stresses are responsible for crack growth such that only the maximum stress leads to crack propagation. Both strategies have in common that the crack surface shows a curved breakage for all degradation functions.
Eventually, we investigate the convergence properties of solution strategies with respect to the length-scale parameter and the choice of the degradation function ( ). Our study considers two nonlinear solution strategies, namely Newton's method configured as described in Section 4.1 and the RMTR, described in Section 4.1.2. [47] Both methods are used to solve nonlinear coupled problems arising within the monolithic solution scheme. We terminate the solution process, if the stopping criterium ‖ ‖ < 10 −7 is satisfied. Here, the symbol ‖ ‖ represents the Euclidean norm of the residual on the th iteration. Table 11 demonstrates the comparison of Newton's method with the RMTR method by varying the choice of the degradation function. The comparison is performed in terms of the cumulative number of nonlinear iterations/V-cycles which transfers the fine level residual to the coarse level within the multigrid method over all time steps. [40,42] As we can see, the behavior of the RMTR method is very robust with respect to the choice of parameter . The number of V-cycles increases very slightly as decreases. However, the performance of Newton's method deteriorates rapidly with lower values of . This behavior of Newton's method is in agreement with measurements obtained for the simple tension test, Section 4.3.1.
The comparison of Newton's method with the RMTR method by varying the choice of degradation function ( ) is illustrated in Table 12. We again observe that the performance of the RMTR method with respect to the length-scale parameter is very robust. The number of accumulative V-cycles stays almost constant as the ratio between and mesh size ℎ changes. In contrast, the convergence of Newton's method deteriorates as the approximation of the crack surface becomes more accurate, for example, as decreases.
From the obtained results we can conclude, that the nonlinear multilevel method, that is, RMTR, outperforms the standard single level Newton's method for all presented test cases. Hence, the use of the nonlinear multilevel method, RMTR, is significantly more efficient than applying multilevel based linear preconditioner on each Newton's iteration. In addition, the RMTR method shows the robustness with respect to the phase-field parameters, when applied to the conchoidal fracture example.

CONCLUSION
The phase-field fracture approach is meanwhile an established way to simulate crack propagation. However, the related parameters have to be chosen carefully. In this contribution the influence of the model parameters (length-scale parameter, mobility, degradation function) has been investigated. The recommended choice of the length-scale parameter is very small such that the diffuse transition zone becomes very small and converges to a sharp crack interface. Since the mesh has to be able to resolve the transition zone, the mesh size limits the choice of . Thus, the cost-benefit ratio has to be considered thoroughly. For large values of the mobility the numerical regularization effect vanishes. Conversely, too low values of the mobility lead to a suppressed phase-field evolution. Furthermore, the choice of the degradation function influences the phase-field model concerning the prescribed displacement and the required force at crack initialization. For smaller parameters the required force and the prescribed displacement increase. All these characteristics of the model parameters perform very similar regardless of the model of the crack driving force. The focus has been set on the standard variational crack driving force and the Rankine failure criterion. Moreover, the solution process has been investigated. Basis of the discussions is the multilevel method for both linear and nonlinear theory considering both solution schemes, the monolithic and staggered solution scheme. To sum up the main results it is observed that the number of iterations within Newton's method decreases for a decreasing length-scale parameter, a decreasing parameter and an increasing mobility parameter. As an exception the preconditioned CG method has to be mentioned in which the behavior within the monolithic and the staggered scheme is reversed. Consequently, all solution strategies are influenced by the considered model parameters.
Finally, we point out the applicability of the phase-field model to a three dimensional conchoidal fracture example applying structured and unstructured meshes.