A quasi-monolithic phase-field description for mixed-mode fracture using predictor–corrector mesh adaptivity

In this work, we develop a mixed-mode phase-field fracture model employing a parallel-adaptive quasi-monolithic framework. In nature, failure of rocks and rock-like materials is usually accompanied by the propagation of mixed-mode fractures. To address this aspect, some recent studies have incorporated mixed-mode fracture propagation criteria to classical phase-field fracture models, and new energy splitting methods were proposed to split the total crack driving energy into mode-I and mode-II parts. As extension in this work, a splitting method for masonry-like materials is modified and incorporated into the mixed-mode phase-field fracture model. A robust, accurate and efficient parallel-adaptive quasi-monolithic framework serves as basis for the implementation of our new model. Three numerical tests are carried out, and the results of the new model are compared to those of existing models, demonstrating the numerical robustness and physical soundness of the new model. In total, six models are computationally analyzed and compared.


Introduction
In classical fracture mechanics, three fundamental modes of fractures are defined. They are termed as tensile (mode-I), in-plane shear (mode-II) and anti-plane shear (mode-III), respectively [7]. But the fractures we deal with in reality are mostly mixed-mode. For example, mixed mode-I and mode-II fracture propagation is very common in rocks under compression [13,31]. Because of the complex formation mechanisms and propagation behaviors of mixed-mode fractures, their simulations require mixed-mode fracture propagation criteria and robust numerical methods.
In terms of numerical methods, phase-field methods have seen a tremendous popularity in fracture propagation simulations in recent years. Such variational methods were first proposed by Francfort and Marigo [16]. A subsequent numerical realization was derived by Bourdin et al. [10]. Therein, two equations are coupled: one to solve displacement and the other to solve the phase-field unknown. The latter is a regularized equation with a model parameter that characterizes the thickness of a smoothed transition zone from full fracture to the unbroken material.
In terms of mixed-mode fracture propagation criteria, Shen and Stephansson modified Griffith's criterion (G-criterion) [19], which can only address single-mode fracture propagation, into an F-criterion. Zhang et al. [40] and Bryant and Sun [11] then incorporated the F-criterion into the phase-field fracture model. Accordingly, the crack driving energy in the original phase-field formulation is split into mode-I and mode-II parts. This is achieved by splitting the crack driving stress + into + I and + II , or in some cases, splitting the crack driving strain e + into e + I and e + II . Zhang et al. modified the energy splitting method for single-mode fractures [29] so that the first and second terms of + became 1 3 + I and + II respectively [40]. This method can distinguish primary cracks from secondary cracks in uniaxially compressed limestones, but it overestimates the maximum force response when mode-I loading is dominant [40]. Strobl and Seelig [34] and Steinke and Kaliske [33] proposed more physically sound splitting methods, where the asymmetry effects of tension and compression normal stresses on fracture propagation are fulfilled explicitly 1 . However, determining fracture surface directions remains unsolved in these two studies. Bryant and Sun derived a similar splitting method considering the asymmetry effects of tension and compression normal strains on fracture propagation [11]. The direction of fracture surface is determined by solving a local fracture dissipation maximization problem [11]. In the numerical experiments, we have carried out, the latter three splitting methods [11,33,34] showed some numerical instabilities. Newton iteration steps increase sharply when the crack starts to propagate and in some cases divergence takes place. Freddi and Royer-Carfagni derived a splitting method for masonry-like materials based on structured deformation theory [17]. This method has a similar formulation to those proposed by Strobl and Seelig [34] and Steinke and Kaliske [33], but seem numerically more robust. Recently, another splitting model was proposed by [9]. However, it can only deal with mode-I fracture as it is derived under the principal strain coordinate and there is no mode-II fracture driving energy term in the formulation.
The objective of the current work is two-fold. First, a novel phase-field model for mixed-mode fractures is introduced. Second, this new model is implemented in wellproven robust and efficient numerical framework, allowing for parallel computations and adaptive mesh refinement. We explain both novelties in more detail in the following.
First, we develop a splitting method, which splits the crack-driving energy into mode-I and mode-II parts. Specifically, we change the reference frame of the energy splitting method proposed by Freddi and Royer-Carfagni [17] from principal strain coordinates to local fracture surface coordinates. We then incorporate this method into the mixed-mode phase-field fracture formulation proposed by Zhang et al. [40] and develop a physically sound and numerically robust phase-field model for mixed-mode fracture propagation. The fracture surface direction is determined based on the maximum-dissipation criterion [11]. We further develop a simpler computational algorithm to implement this criterion, as derived in Sect. 3.6. Second, the resulting model is then implemented in a framework developed in [20,21], named pfm-cracks [22]. The resulting approach is robust, accurate and efficient using parallel computations and adaptive mesh refinement employing a predictor-corrector scheme to detect the priori unknown fracture path. Three numerical tests are carried out under different mesh sizes h and transitional zone sizes , and the results of our new model are compared with those of the existing models. These comparisons demonstrate the numerical robustness and physical soundness of our new model. Preliminary results of this work are published in conference proceedings [15].
The outline of this paper is as follows: in Sect. 2, the notation and governing equations for mixed-mode phase-field fracture model are presented. Next in Sect. 3, several energy splitting methods for mixed-mode fractures are reviewed and our new splitting method is derived. In Sect. 4, the detailed algorithm and numerical scheme used to solve our model are presented. Then in Sect. 5, three numerical tests based on five existing splitting methods and our new splitting method are presented and analyzed. Finally, in Sect. 6, the contents of this paper are summarized and conclusions are drawn.

Governing equations for mixed-mode
phase-field fracture models

Notation
In this subsection, we introduce the basic notation used in the paper. Let B ⊂ ℝ 2 denote the whole two-dimensional domain, inside which a lower-dimensional sharp crack C ⊂ ℝ exists. Following the idea of Bourdin et al. [10], a phase-field variable ∈ [0, 1] is introduced, where = 0 represents the fully damaged state and = 1 represents the undamaged state. Then C is approximated using F ⊂ B with the help of elliptic functional proposed by Ambrosio and Tortorelli [2,3]. The boundary of the regularized fracture is denoted as F , whose position is dependent on the characteristic length . Now the entire domain can be divided into the intact part I ⊂ B and the fractured part F . More specifically, Figure 1 gives an illustration of the geometric setting notation. After applying Dirichlet boundary condition on the domain outer boundary B , we are now faced with a twofield problem depending on the displacement u ∶ B → ℝ 2 and phase-field variable ∶ B → [0, 1] . Finally, we denote the L 2 scalar product: (1) for vectors, (m, n) ∶= ∫ B m ⋅ ndx ; (2) for second-order tensors, (M, N) ∶= ∫ B M ∶ Ndx . (1) 1 Although these two papers are not aimed at solving mixed-mode fracture problem, the proposed splitting methods are capable of splitting + into + I and + II .

Governing equations
In this subsection, we start with the original phase-field fracture model and then introduce its evolvement into the mixedmode phase-field fracture model. Francfort and Marigo generalized Griffith's theory to form a variational model capable of treating crack nucleation, initiation and growth [16]. Bourdin then introduce a phase-field variable to numerically solve the model [10]. Miehe et al. [30] and Amor et al. [4] went one step further and proposed thermodynamically consistent phase-field fracture models in which elastic potential energy is split into fracture-driving and non-fracture-driving parts. Based on the aforementioned literature, the governing equations for thermodynamically consistent phase-field fracture model can be constructed starting from the energy functional E(u, ): where G c is the critical energy release rate, and is a small numerical parameter to prevent singularity in Jacobian matrix. e(u) denotes Cauchy strain. + (u) and − (u) denote the crack-driving and non-crack driving stresses, respectively, and they are fully partitioned from the total stress (u): The material is assumed to be linear-elastic: where and are Lamé coefficients. Since we are dealing with two-dimensional cases in this paper, we want to emphasize that has different formulations in plane stress and plane strain settings: (u) = + (u) + − (u).
(5) (u) = tr(e(u))I + 2 e(u), in plane stress setting in plane strain setting, while has the same formulation in both settings: where E denotes Young's modulus and denotes Poisson's ratio.
To determine the crack propagation path, we need to find the u and which minimize E(u, ) . In the meantime, an inequality constraint must be enforced so that crack irreversibility is fulfilled. For the following, let us introduce three function spaces: on . This minimization problem with inequality constraint can be described as Formulation 1 Let the time step sequence be denoted by the index n = 0, 1, 2, … . Find vector-valued displacements u ∶= u n ∈ {u D + V} and a scalar-valued phase-field variable ∶= n ∈ W in such that The variational form of the constrained minimization problem in Formulation 1 constitutes the governing equations for phase-field fracture model by differentiating E(u, ) with respect to u and : Formulation 2 Let the time step sequence be denoted by the index n = 0, 1, 2, … . With u 0 and 0 as the initial conditions, at each following time step t n , find the displacement and phase-field variables {u, } ∶= {u n , n } ∈ {u D + V} × W in by solving: and Here u D ∶= u D (t) is the prescribed Dirichlet boundary condition at each time step.
Several stress splitting methods have been proposed to compute + (u) and − (u) in Formulation 2, among which the most prevalent ones are those proposed by Miehe et al. [30] and Amor et al. [4]. The splitting method proposed by Miehe et al. [30] is as follows: (1 − ) 2 + + (u), e(w) + + (w), e(u) where e p1 and e p2 are eigenvalues of the strain tensor e, and the matrix P consists of normalized eigenvectors, where < ⋅ > are Macaulay brackets. The splitting method proposed by Amor et al. [4] is as follows: with where d is the dimension number of the studied case, which in this paper is d = 2.
One problem with Formulation 2 and the above splitting methods is that they cannot explicitly account for mixed-mode fracture propagation. The reasons are as follows.
The phase-field evolution equation Eq. (11) in Formulation 2 is based on the maximum strain energy release rate criterion or G-criterion, which is that the where G is the energy release rate describing the crack driving force and G c is the critical energy release rate describing the crack resistance. According to fracture mechanics, there are three basic types of loadings that a crack can experience, namely mode-I, II and III loadings, and each type of loading corresponds to an energy release rate and a critical energy release rate [5]. In this paper, we only consider mode-I and II loadings because we only deal with two-dimensional cases. When dealing with single-mode cracks, i.e. cracks driven by only one basic type of loading, G and G c in Eq. (19) can be replaced by G * and G * c , respectively, where * stands for I or II, corresponding to mode-I or II energy release rates and critical energy release rates. When dealing with mixed-mode cracks, i.e. cracks driven by a combination of two types of loadings, G in Eq. (19) is a sum of two basic types of energy release rates, and G c is a composite critical energy release rate. For most materials, the composite critical energy release rate G c is hard to obtain. Because G c is a weighted average of G Ic and G IIc but the weight coefficients are hard to determine. Moreover, the splitting methods described by Eqs. (12)-(18) cannot explicitly compute each basic type of (13) − = (tr(e)− < tr(e) >)I + 2 (e − e + ), crack propagates when ∶ G ≥ G c , crack driving stress. The + (u) computed is only a composite crack driving stress and does not have a clear fracture mechanics meaning.
To solve the first problem, Shen and Stephansson [32] suggested to use the F-criterion instead of the G-criterion to predict mixed-mode fracture propagation. The F-criterion can be described as [32]: where where denotes the angle between the current crack surface direction and the crack propagation direction. It can be seen from Eq. (21) that when using the F-criterion, instead of a composite critical energy release rate, we have clearly defined mode-I and II critical energy release rates, i.e. G Ic and G IIc are used. And these two critical energy release rates can be measured using standard methods in laboratories [35,36]. Also, Shen and Stephansson found that the F-criterion gave more accurate fracture propagation paths than G-criterion when they are simulating crack propagation paths in rocks subject to compression with an inclined crack [32].
To solve the second problem, a splitting method capable of computing mode-I and mode-II crack driving stresses should be used. A detailed description of such a method will be given in Sect. 3. In the works of Zhang et al. [11] and Bryant and Sun [11], the F-criterion is incorporated into Eqs. (10) and (11), and the phase-field model for mixed-mode fractures is proposed: With u 0 and 0 as the initial conditions, at each incremental time step n = 1, 2, 3 … , find the displacement and phase-field variables {u, } ∶= {u n , n } ∈ {u D + V} × W in by solving: and Remark 1 If the studied material has the property G Ic = G IIc , Formulation 3 is the same as Formulation 2 since (20) crack propagates when: max 1 3 + I (u) + + II (u) = + (u) . However, it has been recognized that most materials have G IIc much larger than G Ic due to the difference of failure mechanisms [28]. Therefore, the phase-field model for mixed-mode fractures is necessary for most materials.
The remaining problem now is how to determine the mode-I crack driving stress + I and mode-II crack driving stress + II .

Stress splitting methods for mixed-mode phase-field fracture models
In this section, we first review several stress splitting methods that have been proposed in the literature. Then we introduce in Sect. 3.6, our modifications to obtain a robust and physically sound model.

Stress splitting à la Zhang et al. [40]
Based on the method proposed by Miehe et al. [29], i.e. Eqs. (12)-(15), Zhang et al. [40] directly extracted + I and + II from Eq. (12) : This splitting method has the disadvantage of overestimating the maximum force response under tension-dominant loading [40]. We can give an illustration of this disadvantage via a simple example. Assuming an element at the fracture tip is subjected to uniaxial tension stress xx in the x direction, the strain state of the element can be described as In this case, the crack driving energy should be pure mode-I, i.e. + I ∶ e > 0 and + II ∶ e = 0 . But if the previous splitting method is used, we obtain + II ∶ e = 2 (e xx ) 2 , which is even comparable to + I ∶ e = (1 − 2 ) (e xx ) 2 . If the studied material has G IIc > G Ic , this method will underestimate the crack driving energy and thus overestimate the force response.

Stress splitting à la Amor et al. [4]
Inspired by the work of Zhang et al. [40], we can also compute + I and + II by modifying the stress splitting method à la Amor et al. [4]. Specifically, we extract + I and + II from Eq. (16): We again use the uniaxial tension setting to test this modified splitting method. In this case, + II ∶ e = (1 + ) 2 (e xx ) 2 , which is comparable to + I ∶ e = (1 − ) 2 (e xx ) 2 . Therefore, this method also has a force overestimation problem if the studied material has the property G IIc > G Ic . [34] and Steinke/Kaliske [33] The splitting methods proposed by Strobl and Seelig [34] and Steinke and Kaliske [33] do not have the overestimation problem. These two methods, although the former is presented in the strain form and the latter is presented in the stress form, are basically the same. They are both derived under the local crack coordinates and both have the assumptions that -The fracture can sustain compressive, but not tensile, normal stress perpendicular to its surface. -The fracture cannot sustain shear stress tangential to its surface, i.e. the fracture surface is frictionless.

Stress splitting à la Strobl/Seelig
For the convenience of implementation, we only present the strain-form formulation [34]. In 2D case, let n and s denote the directions normal and tangential to the fracture surface, respectively. Their formulation reads: where We again use the uniaxial tension example to test this splitting method. If the fracture propagation direction is perpendicular to the uniaxial loading, which is a fair assumption in homogeneous isotropic material, then e + II = 0 , + II = 0 , and correspondingly + II (u) ∶ e(u) = 0 . Therefore, this splitting method does not have the overestimation problem in uniaxial tension test. However in our experiment, this method is not numerically robust as the Newton iteration step has a sharp increase when crack starts to propagate and sometimes divergence takes place. The reason for this is not fully clear yet. But the singular strain tensors in Eqs. (30) and (31) might lead to a large condition number in the local Jacobian matrix and thus result in an ill-conditioned global Jacobian matrix. Consequently, the nonlinear and linear solvers show bad convergence histories. The performance comparison of Newton solvers incorporating different stress splitting methods are presented in Sects. 5.1 and 5.2.

Stress splitting à la Bryant/Sun [11]
Bryant and Sun derived a similar splitting method based on the assumptions that [11]: -The fracture can sustain compressive, but not tensile, normal strain perpendicular to its surface. -The fracture surface cannot sustain shear strain tangential to its surface.
Although this method is originally presented in energyform [11], it is presented here in an equivalent strain-form so that the differences between different splitting methods are more clear: Then + I and + II are calculated using Eq. (29). This method, similar to the one à la Strobl/Seelig and Steinke/Kaliske, can address the uniaxial example well, but is not numerically robust due to the potential ill-conditioned Jacobian matrix.

Stress splitting à la Freddi/Royer-Carfagni [17]
Freddi and Royer-Carfagni derived an energy splitting method based on the structural strain theory under principal strain coordinates [17]. Its assumption is that the structured, or non-elastic, part of the strain can be represented by a symmetric positive-semidefinite tensor. Let p1 denotes the direction of the eigenvector corresponding to the most-positive eigenvalue. And p2 denotes the direction of the other eigenvector. The method's formulation reads: Equation (35) is very similar to Eq. (30) if p1 and n coincide. It can also be inferred that Eq. (35) is numerically more robust than Eq. (30). Because in the most "dangerous" case where both principal strains are positive, the crack driving strain tensor e + is not singular anymore. This inference will be verified by the numerical test results in Sect. 5. This splitting method can capture the tensile fracture geometry under compression loading [18]. But it does not account for mixed-mode fractures, as there is only e + instead of e + I and e + II in the formulation.

Our stress splitting model
Based on the previous explanations, we perform a heuristic modification resulting in a novel stress splitting model. To this end, the reference frame is changed from principal strain coordinates to local crack coordinates, and e + I and e + II are determined to account for mixed-mode fractures. The new formulation reads: Comparing Eq. (36) with Eq. (30), we see that the main difference between our stress splitting method and the stress splitting method á la Strobl/Seelig and Steinke/Kaliske is that we employ an additional if/else condition when e nn + +2 e ss > 0 and e ss > 0 . Adding this condition influences our stress splitting method in two ways: -From the physical point of view, this additional condition brings in one more assumption that, when the fracture is subject to normal tensile stress perpendicular to its surface and normal tensile strain parallel to its surface, the total elastic energy originated from the normal strain should be dissipated to form fracture surfaces. -From the numerical point of view, this additional condition makes our splitting method more robust than the method à la Strobl/Seelig and Steinke/Kaliske. In the case where e nn + +2 e ss > 0 and e ss > 0 , the strain tensor matrix in the first condition of Eq. (30) might have a large condition number, resulting in an ill-conditioned Jacobian matrix. However, in the same case, the first condition of Eq. (36) gives an strain tensor matrix with smaller condition number, thus resulting in a better-conditioned Jacobian matrix.
Again, + I and + II are calculated using Eq. (29). As for the fracture surface direction determination, we solve an additional local dissipation maximization problem [11]. Let denote the angle between n and p1 , we define the normalized crack driving energy as then the normal direction of crack surface can be determined by After a simple derivation, we obtain In all three situations, the local maximum points of F(u, ) are located at = 0 or = 4 or both. Therefore, instead of implementing complicated global maximum search algorithm, we can compare F(u, 0) and F(u, 4 ) and choose the that maximizes F(u, ) as n . The details about the implementation will be presented in Sect. 4.4.

Finite-element implementation and solution algorithms
In this section, we first briefly introduce our spatial discretization strategy. Then we introduce a primal-dual active set method used to treat crack irreversibility. This approach can be characterized as a semi-smooth Newton method [23] and is then combined with a classical Newton � G Ic sin 2 if e nn + +2 e ss > 0 and e ss ≤ 0, technique, resulting in a combined Newton algorithm [20,21]. Then our specific forms of the Newton residual and the Jacobian matrix are presented. Finally, we derive a new algorithm for implementing our suggested splitting method. The resulting implementation is an extension of the open-source phase-field fracture code [21] (pfmcracks [22]) that uses finite elements from deal.II [6].

Spatial discretization
A Galerkin finite element method is used to spatially discretize the 2D domain. Quadrilateral elements with bilinear shape functions ( Q 1 ) for both u and are used. Consequently, we work with the discrete spaces V h ⊂ V , and W h ⊂ W . For the locally adapted meshes, we work with hanging nodes, see, e.g. the work of Carey and Oden [12].

Primal-dual active set method and a combined Newton algorithm
In this subsection, we follow closely the work of Heister et al. [20] and recapitulate the main steps for Newton's method solving an unconstrained minimization problem.
To this end, the primal-dual active set method is used to treat crack irreversibility and the final algorithm which combines the former two methods into one iteration at each time step.

Newton's method for a minimization problem
First, we rewrite Formulation 1 into a more compact form: where U ∈ {u D + V} × W in is the solution of current time step, and Ū is the solution from last time step.
The unconstrained minimization problem in Eq. (41) can be solved using Newton's method. A Newton iteration sequence U 0 , U 1 , … , U N can be constructed in each time step, with where the update U k is computed as the solution of the linear system: if at each Newton iteration step, the condition is fulfilled, and let the initial guess U 0 =Ū , then the crack irreversibility constraint in Eq. (42) can be fulfilled since

A primal-dual active set method
In the following, the notation is shorten by dropping the index k of Newton iteration, i.e. U ∶= U k , and setting . Thus, the system described by Eqs. (44) and (45) can be rewritten as a minimization problem: Following the work of Hintermüller et al. [23], the minimization problem Eq. (46) can be solved using a primal-dual active set strategy. Using a Lagrange multiplier ∈ 0 × W * (where W * is the dual space of W), the minimization problem can be described using a equation system: where for a given c > 0 . The primal-dual active set strategy replaces the condition C( U, ) = 0 with U = 0 on in the active set A and = 0 in the inactive set I . In other words, the active set is the subdomain where the constraint is applied and no PDE is solved. The inactive set is the subdomain where constraint is already fulfilled and PDE is solved. The algorithm for computing active set reads: Algorithm 1 In each time step, with U 0 and 0 as initial condition, repeat for k = 1, 2, … until the active set A k does no longer change: 1. Compute active set: 2. Compute U k ∈ V × W and k ∈ 0 × W * using: A discretized version of step 2 leads to a linear system with block structure: Using quadrature only on the support points of k h , B block becomes diagonal. The equations B T U k h = 0 is not actually solved, but handled via linear constraints used to eliminate equations in the G block on active set A k . The eliminated equations are those where the ith entry of k h is non-zero. Therefore, the linear equation system simplifies to where Ĝ and F stem from G and F by removing the constrained rows from the system.
Finally, each entry of k h can be computed using: Then k h is used to compute the active set A k+1 . The index i is in the active set A k+1 if and in the inactive set I k+1 otherwise.

A combined semi-smooth Newton solver
In the implementation, two Newton iterations, namely the active set computation iteration and the PDE solving iteration, are combined into one iteration to calculate the variable U k . The combined iteration contains a back-tracking line search to improve the convergence radius. Moreover, Eq. (55) is rewritten using the earlier notation: Finally, the algorithm combining two iterations can be described as

3
Algorithm 2 In each time step, with U 0 , U 0 and 0 as initial conditions, repeat for k = 1, 2, ⋯ until the active set A k does no longer change and R U k h < TOLR U 0 h : . Eliminate rows and columns in A k from G and F to obtain Ĝ and F . 5. Solve linear system Ĝ U k =F. 6. Find a step size 0 < ≤ 1 such that Remark 2 In the stopping criterion in Algorithm 2, the residual on the inactive set R U k h instead of full residual R U k h is used. The former is obtained by eliminating the constrained rows in the latter.

Formulation of the Newton residual and the Jacobian matrix
In this subsection, we introduce the specific forms of the Newton residual and Jacobian matrix. A monolithic scheme in which all equations are solved simultaneously is used to solve the aforementioned PDE system. However, the energy functional in Eq. (3) and its derived mixed-mode form are not convex with respect to both u and , but convex with respect to only u or if the other is fixed. The critical term is the cross term (1 − ) 2 + + (u), e(u) . Heister et al. linearized in the cross term by a time-lagging interpolation scheme ≈̃∶=̃( n−1 , n−2 ) [20]. A simple sub-iteration with a few cycles can be adopted to increase the accuracy comparable to fully monolithic schemes [38].
However, in this work, we choose a more conservative linearization scheme that ≈̃∶= n−1 to further increase the stability of the solver. Finally, the Newton residual (right-hand side) is written as ,

Numerical tests
In this section, we compare five existing energy splitting methods with our new splitting method. Therein, we compare the simulated crack geometry and load-displacement curves. The numerical performance is evaluated by observing Newton iteration numbers.
� if e nn + +2 e ss > 0 and e ss > 0, if e nn + +2 e ss > 0 and e ss ≤ 0, 0 if e nn + +2 e ss ≤ 0,  The uniaxial compression test with an inclined notch is chosen to show that our new method can best capture the mixedmode fracture geometry observed in the experiment. In each test, three uniformly refined meshes are used to guarantee the reliability of the results. Two adaptively refined meshes are also used to demonstrate the compatibility of our new method with adaptive mesh refinement. To guarantee that the fracture profile width becomes smaller as mesh becomes finer-which is the main reason to use fine mesh-we let to change adaptively as the mesh is refined. To be more specific, we set = 2h min , where h min is the minimum mesh size.

Configuration and boundary conditions
The configuration of this test is illustrated in Fig. 2. We emphasize that the fracture generated under this configuration should be pure mode-I and the incorporation of G IIc , whatever value it has, should not change the crack propagation behavior. In the uniform mesh cases, the initial domain is five, six and seven times uniformly refined, leading to 4096, 16384 and 65536 mesh cells, and the minimum mesh size h min = 0.022mm, 0.011 mm and 0.0055 mm, respectively. In the adaptive-refined mesh cases, the domain is first five times uniformly refined and then one and two times adaptively refined. The threshold for adaptive mesh refinement is set as threshold = 0.5.
For the displacements, a non-homogeneous timedependent Dirichlet boundary condition is applied on the top boundary: where t denotes the total loading time. And the incremental time step size is set as t = 1.0 × 10 −4 s. For the phase-field function, homogeneous Neumann boundary conditions are applied on the entire boundary.

Quantities of interest
Apart from crack geometry, we also check the average stress normal to the top boundary: where top ∶= {(x, y) ∈ B | 0 mm ≤ x ≤ 1 mm, y = 1 mm} is the top boundary, n is the outward unit vector normal to the top boundary surface and H is the Hausdorff measure. We also want to emphasize that, due to the influence of boundary conditions, the fracture profile enlarges when the fracture touches the boundary. This phenomenon has been reported by previous studies in single-edge notched tension test [29], single-edge notched shear test [20], and pressurized crack propagation test [37]. Similarly, we observed this phenomenon in all of our three tests, as shown in Figs. 3, 11 and 20.  -The Newton solver incorporating mixed-Freddi's method has better performance than that incorporating Freddi's method. Between the 63th and 114th time steps mixed-Freddi's case has comparable Newton step numbers with Miehe's case. While Freddi's case has Newton step number two times more than that in Miehe's case. Strobl's cases, the Newton solver stagnates before the load decreases to zero. -In all but mixed-Miehe's cases, the load-displacement curves of the same splitting method are consistent under three different meshes. But the maximum load response increases as mesh refinement increases and decreases. This phenomenon has also been observed in the work of Zhang et al. [41] However in mixed-Miehe's case, the finest mesh, seven times global refinement, has the intermediate load response, indicating the instability of the method.

Discussion of the load-displacement curves
In Fig. 6, we observe the load-displacement curves of the six splitting methods when the domain is seven times refined. It shows that -All but Bryant's and mixed-Miehe's methods have almost the same load-displacement curves. Particularly, the load-displacement curves of Strobl's, Freddi's and mixed-Freddi's are almost identical that they overlap with each other. The only difference among them is that the curve terminated earlier in Strobl's case, which can also be observed in Fig. 5. Mixed-Freddi -In Bryant's case, the load-displacement curve terminates early due to the Newton solver non-convergence.
In mixed-Miehe's case, the maximum load response is about 37% higher than that in other cases. This phenomenon verifies the overestimation problem of mixed-Miehe's method, as mentioned in Sect. 3.1.
We remark that although Bryant's method worked well in uniaxial compression tests with single and multiple notches in [11], it appeared to be not numerically robust in all of our tests, namely the single-edge notched tension, single-edge notched shear and uniaxial compression with an inclined notch. But Bryant's model is still a significant reference for this paper as it gives inspiration for how to determine the crack surface direction. Figure 7 shows the crack geometry of mixed-Freddi's splitting method under five times global refinement plus two adaptive refinement cycles (5+2). Comparing Fig. 7 with the last row in Fig. 3, we observe that the crack in 5+2 refinement case has almost the same propagation progress and smeared width as those in 7+0 refinement case. In Fig. 8 the load-displacement curves of mixed-Freddi's method under 5+0, 5+1 and 5+2 refinement are displayed. Comparing Fig. 8 with the last sub-figure in Fig. 5, we can see that the 5+1 and 5+2 refined cases have very similar load responses with the 6+0 and 7+0 refined cases.

Discussion of the adaptive mesh refinement
In conclusion of this first numerical example, mixed-Freddi's method works well with adaptive mesh refinement. Here, almost the same crack geometry and load response of uniform refinement case can be captured, but with much less degrees of freedom (DoFs). For example, the 7+0 case has 198531 DoFs through the entire computation process, while the 5+2 case has only 12771 DoFs at the 1st time step and 18903 DoFs at the 114th time step. Illustration of meshes for the 7+0 case (only uniformly refined) and 5+2 case (five times uniformly refined and two times adaptively refined) for the single-edge notched tension test are displayed in Fig. 9.

Configuration and boundary conditions
The configuration of this test is illustrated in Fig. 10. This test is widely used to check whether an energy splitting method is capable of only inducing fractures in the tensile-loaded region. If the splitting method fails to replicate this behavior or if no splitting method is used, we will get an incorrect branched crack geometry, as shown in for instance the work of Bourdin et al. [10] We emphasize that the fracture generated under this configuration is also pure mode-I although shear load is applied at the top  boundary. The maximum principal stress in the tensileloaded region drives the crack to propagate [11].
The boundary condition is similar to that in the tension case, but the loading direction is now tangential to the top boundary: Correspondingly, we now check the average stress tangential to the top boundary: where s is the unit vector tangential to the top boundary surface.

Material and model parameters
The mesh refinement strategy, mechanical and numerical parameters are the same as those in the single-edge notched tension test presented in the previous subsection. Figure 11 displays the crack geometry of different splitting methods when the domain is seven times uniformly refined. We make the following observations:

Discussion of our findings with respect to crack geometries
(79) u x (t) =ūt withū = 1 mm/s.
-All methods except Bryant's can correctly capture the curved fracture path propagating from the center to the lower-left corner. The Newton's solver incorporating Bryant's method does not converge when the fracture initiates at the center so the computation ends at u = 8.9 × 10 −3 mm. -The fracture propagation in mixed-Miehe's case is noticeably lagged compared with that in Miehe's case. While the fracture propagation in mixed-Freddi's case is almost, if not exactly, the same as that in Freddi's case.

Discussion of Newton's performance
In Fig. 12, the Newton iterations are monitored for different splitting methods when the domain is seven times uniformly refined. The Newton solver performances of different splitting methods are similar to those in the tension case. In detail: -The Newton solvers incorporating Bryant's, Strobl's and Freddi's methods still have much higher iteration numbers at certain time steps than those incorporating other splitting methods. -The Newton solver incorporating mixed-Freddi's method has approximately 10 Newton iteration steps at each time step when the crack is propagating, showing a better performance than the solver incorporating Freddi's method. Figure 13 shows the load-displacement curves of the six splitting methods under three uniformly refined meshes. Again, the findings are similar to those in single-edge notched tension test. In detail, -All but Bryant's and Strobl's methods can capture the full loading-unloading process. -In all cases, the maximum load response increases as mesh refinement increases and decreases. -After the crack reaches the boundary, the tangential load, instead of going to zero, starts to increase. This nonphysical phenomenon is caused by the elastic response of the cracked specimen clamped at the undamaged lower right portion of the boundary [1].

Discussion of the load-displacement curves
We remark that if Amor's method [4] and Bryant's method [11] are used, the stresses can go completely to zero in the single-edge notched shear test when the crack reaches the boundary [1,11]. However, in our shear tests: (1) Amor's method results in an undesirable sub-horizontal crack geometry; (2) Bryant's method is not numerically robust and only works in very limited cases. Figure 14 displays the load-displacement curves of the six splitting methods when the domain is seven times refined. We observe -All but Bryant's and mixed-Miehe's methods have similar load-displacement curves. Specifically, the loaddisplacement curve of mixed-Freddi's method overlaps with that of Strobl's method, and it only diverges from that of Freddi's method after the crack reaches boundary and load starts to re-increase. -In Bryant's case, the load-displacement curve terminates early due to the Newton solver non-convergence.
In mixed-Miehe's case, the maximum load response is about 59% higher than that in other cases.

Discussion of the adaptive mesh refinement
Similar to the single-edge notched tension test, Figs. 15 and 16 show that employing adaptive mesh refinement reproduces the findings obtained with uniform refinement, but with much less DoFs. Illustration of meshes for the 7+0 case (only uniformly refined) and 5+2 case (five times uniformly refined and two times adaptively refined) for the single-edge notched shear test are displayed in Fig. 17.

Uniaxial compression test with an inclined notch
In this test, we only present and compare Miehe's, mixed-Miehe's, Freddi's and mixed-Freddi's methods since Mixed-Freddi  Bryant's and Strobl's methods appear to be less robust in the previous tests.

Preliminary considerations
As it is already reviewed by Wong [39], various experiments have been done to study the crack geometry in specimens containing a single flaw under uniaxial compression. To choose one experiment from them as the benchmark for our numerical test, we make the following considerations: -We want to find an experimental setting which is featured with distinct mode-I and mode-II fractures. -Since we only consider two-dimensional situations in the present study, we want to find an experimental configuration satisfying the plane stress setting. More specifically, the dimension of the specimen in the z direction should be much smaller than the dimensions in the x and y directions. And the crack should penetrate through the specimen in the z direction. -Since we do not consider the friction force on the fracture surface, we want to find an configuration with an open flaw, i.e. a notch, so that this frictionless simplification has minor influence on the comparison between simulation and experiment results.
A first possibility is the limestone experimental setting presented by Ingraffea and Heuze [26], which is chosen as the numerical test case in the work of Zhang et al. [40]. However, in this setting, the primary and secondary cracks are all mode-I cracks [25,26]. Instead, we choose the marble experimental setting used by Huang et al. [24], in which distinct mode-I and mode-II fractures are induced. As shown in Fig. 18, the crack geometry of this setting can be summarized as -Mode-I cracks: 1-primary forward tensile cracks (PFTCs), 2-secondary forward tensile cracks (SFTCs) and 3-backward tensile cracks (BTCs). -Mode-II cracks: 4-forward shear belts (FSBs) and 5backward shear belts (BSBs).
The formation of FSBs and BSBs is caused by the coalescence of micro-cracks [24]. These micro-cracks are essential to increasing the production rates in the oil and gas industry as they can largely increase the permeability of rocks (see for instance the work of Zoback and Byerlee [42]).

Configuration and boundary conditions
The configuration of this test is illustrated in Fig. 19. The uniformly refined base mesh is first generated with 7952 mesh cells. On this basis, (1) in the uniformly refined mesh cases, the domain is zero, one and two times uniformly refined, leading to 7952, 31808 and 127232 mesh cells, and the minimum mesh size h min = 1.17 mm, 0.59 mm and 0.29 mm , respectively. (2) In the adaptively refined mesh cases, the domain is further one and two times adaptively refined. The threshold for adaptive mesh refinement is set as threshold = 0.5.
For the displacements, a compressive Dirichlet boundary condition is applied at the top boundary: And the incremental time step size is set as t = 1.0 × 10 −3 s. For the phase-field function, homogeneous Neumann boundary conditions are applied on the entire boundary.

Materials and model parameters
We consider plane stress setting and use the mechanical parameters as follows: Young's modulus E = 63.5 kN/mm 2 [24] and Poisson's ratio = 0.21 . In the cases using Miehe's method and Freddi's method, G c = 11.0 × 10 −6 kN/mm, otherwise G Ic = 11.0 × 10 −6 kN/mm and G IIc = 10G Ic [40]. The parameters , and the Newton solver tolerance are set to the same values as in the previous two tests.

Quantities of interest
Apart from the crack geometry, we also check the average compressive stress normal to the top boundary: In this test, the zero and one time uniformly refined cases have similar crack propagation behaviors. But the two times uniformly refined case behaves differently. Therefore, we present the crack geometries of both one and two times uniformly refined cases.

Discussion of our findings with respect to crack geometries
In Fig. 20 (37) is non-zero and drives the mode-II micro-cracks to propagate. The only discrepancy between the results of mixed-Freddi's simulation and the results of experiment is that there are no SFTCs. This discrepancy is acceptable because in the experiment SFTCs propagate in a very similar manner as PFTCs, and the initiation positions of SFTCs and PFTCs are very near to each other. Figure 21 displays the crack geometries when the domain is two times uniformly refined. The main differences between Figs. 20 and 21 can be summarized as follows: -In Miehe's and mixed-Miehe's cases, the PFTCs start branching at stage(b). Thus, the crack geometries in these two cases become more similar to each other. -In Freddi's and mixed-Freddi's cases, PFTCs stop propagating after having reached certain lengths. The load response continue increasing until the solver dies. The specimen does not rupture at the end of the simulation.
The reason for crack branching in Miehe's and mixed-Miehe's cases might be that the cracks become thinner as decreases, so the originally merged neighboring cracks now become distinct from each other and thus exhibiting the branching behavior. The reason for the non-rupture phenomenon in Freddi's and mixed-Freddi's cases is more complicated. One physical interpretation is that the crack smeared width can represent the width of micro-crack zones surrounding the propagating macro-cracks. In rocks, there are usually wide micro-crack zones around macro cracks  [14]. But in more isotropic and homogeneous materials like poly(methyl methacrylate) (PMMA), the macro cracks are sharper and there are very thin or even no micro-crack zones [8]. In other words, when other parameters are fixed, large can better describe the cracking behavior of rocks while small is more suitable for PMMA. Therefore, in the one time refined case, the simulated crack geometries are similar to those in the marble experiment, where the specimen ruptures as cracks propagate to the boundary. While in the two times refined case, the simulated crack geometries are similar to those in the PMMA experiment [25][26][27], where PFTCs stop propagating after having reached certain length and specimen does not rupture even when load increases to a very high level. We also want to emphasize that in this extreme case, the mixed-Freddi's method is numerically more robust than Freddi's method, as the PFTCs can propagate further in mixed-Freddi's case. Figure 22 monitors the load-displacement curves of the four splitting methods under three times uniformly refined meshes. We observe -As it is already stated in Eq. (5), we use a linear elastic constitutive law in the phase-field model. However, in Miehe's and mixed-Miehe's cases, the load-displacement curves exhibit softening behavior before peak loading, which might be caused by the wide-smeared crack profile. In Freddi and mixed-Freddi's cases, the load-displacement curves are straight lines from zero loading to peak loading, which is a typical behavior of linear-elastic brittle material. Therefore, we can conclude that Freddi's and mixed-Freddi's method can better capture the constitutive behavior in this test.  -When the domain is zero and one time uniformly refined, the simulation continues until the load starts to decrease (when cracks reach boundaries) in all four cases. When the domain is two times uniformly refined, the simulation still continues until the load starts to decrease in Miehe's and mixed-Miehe's cases. While in Freddi's and mixed-Freddi's cases, the simulation ends before the load decreases. Again we can see that mixed-Freddi's method is numerically more robust than Freddi's method, as the load-displacement curve extends further in mixed-Freddi's case.

Discussion of the adaptive mesh refinement
As before, Figs. 23 and 24 show that tests using adaptive refinement can reproduce the results in tests using uniform refinement, but with much less DoFs, as shown in Fig. 25.

Conclusions
In this work, we developed on the basis of existing studies [11,17,33,34,40], a phase-field fracture model for capturing mixed-mode fracture propagation behavior. We used the phase-field fracture governing equations which incorporated the F-criterion [40]. Additionally, we changed the reference frame of the energy splitting method proposed by Freddi and Royer-Carfagni [17] from principal strain coordinates to local fracture surface coordinates so that the crack-driving energy was further split into mode-I crack driving part and mode-II crack driving part. The fracture surface direction was determined based on the maximum-dissipation criterion [11], and a simpler algorithm is proposed by us. Our computational framework was based on the implementation of pfm-cracks [22], where a quasi-monolithic scheme was used to solve the governing equations and adaptive mesh refinement and parallel computations were employed to reduce computation time. In this framework, a semi-smooth Newton method was used to solve the non-linear PDEs and treat the crack irreversibility. We hence provided the detailed algorithm to compute the right-hand side and Jacobian matrix of Newton's method. Finally, to show the superiority of our new model over existing models, we used two classical numerical tests to demonstrate its numerical robustness and one practical test to demonstrate its physical soundness. We also showed the good compatibility of our new model with adaptive mesh refinement for reducing the computational cost. In the future, we are interested in extending the new model to three-dimensional scenarios so that more practical settings from engineering applications can be treated.  otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.