An iterative staggered scheme for phase field brittle fracture propagation with stabilizing parameters

This paper concerns the analysis and implementation of a novel iterative staggered scheme for quasi-static brittle fracture propagation models, where the fracture evolution is tracked by a phase field variable. The model we consider is a two-field variational inequality system, with the phase field function and the elastic displacements of the solid material as independent variables. Using a penalization strategy, this variational inequality system is transformed into a variational equality system, which is the formulation we take as the starting point for our algorithmic developments. The proposed scheme involves a partitioning of this model into two subproblems; phase field and mechanics, with added stabilization terms to both subproblems for improved efficiency and robustness. We analyze the convergence of the proposed scheme using a fixed point argument, and find that under a natural condition, the elastic mechanical energy remains bounded, and, if the diffusive zone around crack surfaces is sufficiently thick, monotonic convergence is achieved. Finally, the proposed scheme is validated numerically with several bench-mark problems.


Introduction
Fracture propagation is currently an important topic with many applications in various engineering fields. Specifically, phase-field descriptions are intensively investigated. The theory of brittle fracture mechanics goes back to the works of A. Griffith [23], wherein a criterion for crack propagation is formulated. Despite a foundational treatment on the subject of brittle fracture, Griffith's theory fails to predict crack initiation. This deficiency can however be overcome by a variational approach, which was first proposed in [10,20]. Using such a variational approach, discontinuities in the displacement field u across the lower-dimensional crack surface are approximated by an auxiliary phase-field function ϕ. The latter can be viewed as an indicator function, which introduces a diffusive transition zone between the broken and the unbroken material. The enforcement of irreversibility of crack growth finally yields a variational inequality system, of which we seek the solution {u, ϕ}.
In this work, we concentrate on improvements of the nonlinear solution algorithm, which is still a large bottleneck of phase-field fracture evolution problems. Specifically, high iteration numbers when the crack initiates or is further growing are reported in many works [21,29,44,45]. However, in most studies iteration numbers are omitted. Both staggered (splitting) schemes and monolithic schemes are frequently employed. Important developments include alternating minimization/staggered schemes [9,11,12,29,30], quasi-monolithic scheme with a partial linearization [25], and fully monolithic schemes [21,44,45].
The goal of this work is to propose a linearized staggered scheme with stabilizing parameters. In particular, the proposed scheme is based on recent developments on iterative splitting schemes coming from poroelasticity [13,26,32,33]. Iterative splitting schemes are widely applied to problems of coupled flow and mechanics, where at each iteration step either of the subproblems (i.e., flow or mechanics) is solved first, keeping some physical quantity constant (e.g., fixed stress or fixed strain), followed by solving the next subproblem with updated solution information. This procedure is then repeated until an accepted tolerance is reached. Further extensions of this technique involves tuning some artificial stabilization terms according to a derived contraction estimate in energy

Preliminaries
In this section we explain the notation used throughout this article, see e.g. [18,47]  For k ∈ N, k ≥ 0, we denote by W k,p (B) the space of functions in L p (B) admitting weak derivatives up to k'th order. In particular, H 1 (B) := W 1,2 (B) and we denote by H 1 0 (B) its zero trace subspace. Note that we reserve the use of bold fonts for second order tensors. Hence, if u, v ∈ L 2 (B), their inner product is (u, v) := B u(x)v(x)dx, and similarly, if u, v ∈ (L 2 (B)) d then we take their inner product to be (u, v) := B u(x) · v(x)dx. Finally, if u, v ∈ (L 2 (B)) d×d then their inner product is (u, v) := B u(x) : v(x)dx.
We will also frequently apply several classical inequalities, in particular: Cauchy-Schwarz, Young, Poincaré, and Korn. See e.g. [15,24] for a detailed description of these.

Governing equations
What follows is a brief description of the phase field approach for quasi-static brittle fracture propagation, see e.g. [20,30] for more details. Consider a (bounded open) polygonal domain B ⊂ R d , wherein C ⊂ R d−1 denotes the fracture, and Ω ⊂ R d is the intact domain, and a time interval (0, T ) is given with final time T > 0. By introducing the phase field variable ϕ : B × (0, T ) → [0, 1], which takes the value 0 in the fracture, 1 in the intact domain, and varies smoothly from 0 to 1 in a transition zone of (half-)thickness ε > 0 around C, the evolution of the fracture can be tracked in space and time. Using the phase field approach, the fracture C is approximated by Introducing the displacement vector u : B×(0, T ) → R d , the model problem we consider arises as a minimization problem: An energy functional E(u, ϕ) is defined according to Griffith's criterion for brittle fracture [23], which is then sought to be minimized over all admissible {u, ϕ}. From this minimization problem, the Euler-Lagrange equations are obtained by differentiation with respect to the arguments, yielding a variational equality system. Finally, a crack irreversibility condition must be enforced (the crack is not allowed to heal), which takes the form ∂ t ϕ ≤ 0. Thus, the variational equality system, which is the previously mentioned Euler-Lagrange equations, is transformed into a variational inequality system, which reads as follows: Find (u(t), ϕ(t)) ∈ V × W := (H 1 0 (B)) d × W 1,∞ (B) such that for t ∈ (0, T ] there holds where G c > 0 is the critical elastic energy restitution rate, 0 < κ << 1 is a regularization parameter, the purpose of which is to avoid degeneracy of the elastic energy (equivalent with replacing the fracture with a softer material), and g(ϕ) := (1 − κ)ϕ 2 + κ is a standard choice for the degradation function (see e.g. [39,45]. Note that g(ϕ) → κ when approaching the fracture zone). The body force acting on the domain B is b : B × (0, T ) → R d , and |Ce(u)| 2 := Ce(u) : e(u) is the elastic mechanical energy, where e(·) := (∇(·) + ∇(·) )/2 is the symmetric gradient, and C = [C ijkl ] ijkl is the fourth order tensor containing the elastic material coefficients, where each C ijkl ∈ L ∞ (B).
We assume that C satisfies the usual symmetry and positive definiteness properties, i.e., (Cu, v) = (u, Cv), and (Cu, u) 1 2 defines an L 2 -equivalent norm, i.e., there exists constants λ m , λ M > 0 such that In order to facilitate the following developments we assume continuity in time for {u, ϕ, b}. Let now 0 = t 0 < t 1 < · · · < t N = T be a partition of the time interval (0, T ), with time step δt := t n − t n−1 , and denote the time discrete solutions by u n := u(·, t n ), (2.3) ϕ n := ϕ(·, t n ). (2.4) The irreversibility condition now becomes ϕ n ≤ ϕ n−1 (using a backward Euler method), and the time-discrete version of the problem (2.1a)-(2.1b) reads as follows: where b n := b(·, t n ). The last term in the phase field equation (2.5b) is a penalization to enforce the irreversibility condition, thus transforming the variational inequality (2.1b) into a variational equality, with penalization parameter γ > 0, and where Ξ ∈ L 2 (B) is given (in practice Ξ will be obtained by iteration, cf. Section 4). Note that we also used the notation [x] + := max(x, 0). From here on, we shall refer to (2.5a) as the mechanics subproblem, and to (2.5b) as the phase field subproblem. Regarding the degradation function g, it is easily seen to satisfy the following Lipschitz condition: The time-discrete system (2.5a)-(2.5b) was analyzed in [36], and there it was shown that at least one global minimizer (u n , ϕ n ) ∈ V × W exists, provided b n ∈ (L 2 (B)) d , for each n. We mention also that the analysis of a pressurized phase field brittle fracture model can be found in [34,35].

Iterative scheme
In this section we introduce the iterative staggered solution procedure for the fully discrete formulation of (2.5a)-(2.5b). To this end, let T h be a simplicial mesh of B, such that for any two distinct elements of T h their intersection is either an empty set or their common vertex or edge. We denote by h the largest diameter of all the elements in T h , i.e., h := max K∈T h diam(K), and let V h × W h ⊂ V × W be appropriate (conforming) discrete spaces. We continue now with the same notation for the variables and test functions as before (omitting the usual h-subscript), since we will from here on mostly deal with the discrete solutions. For each n, the iterative algorithm we propose defines a sequence {u n,i , ϕ n,i }, for i ≥ 0, initialized by {u n−1 , ϕ n−1 }. The iteration is then done in two steps: First, the mechanics subproblem is solved, with the degradation function held constant. Then, the phase field subproblem is solved, with the elastic energy held constant. Note that there are also artificial stabilizing terms which are held constant during solving of the subproblems. Introducing the stabilization parameters L u , L ϕ > 0 (to be determined later), the iterative algorithm reads as follows: • Step 1: Given (u n,i−1 , ϕ n,i−1 , b n ) find u n,i such that Step 2: Given (ϕ n,i−1 , u n,i , ϕ n−1 ) find ϕ n,i such that where, in order to avoid the [·] + -bracket, we also introduced the function η i ∈ L ∞ (B) defined for a.e. x ∈ B by

Convergence analysis
We now proceed to analyze the convergence of the scheme (3.1a)-(3.1b). Our aim is to show a contraction of successive difference functions in energy norms, which implies convergence by the Banach Fixed Point Theorem (see e.g. [14]). To this end we define the following difference functions 3) where {u n , ϕ n } denotes the (exact) solutions to (2.1a)-(2.1b) at time t n . Using the symmetry properties of C, the following set of difference equations are then obtained by subtracting (3.1a)-(3.1b) solved by {u n , ϕ n } from the same equations solved by the iterate solutions: Furthermore, we introduce the following assumption related to the elastic mechanical strain.
Assumption 1 (Boundedness of elastic strain). We assume there exists a constant M > 0 such that Moreover, we assume that M is large enough such that the above bound holds also for the iterate elastic strain, i.e., Note that M is nothing else than an upper bound for the elastic strain in the system for the converged solution, which is arguably finite for any reasonable problem. Note also that with sufficient regularity of the domain, coefficients, source terms, and initial data, the above assumption is satisfied, i.e., the problem (2.5a)-(2.5b) admits a solution u n ∈ (W 1,∞ (B)) d , thus implying the existence of M . Alternatively to introducing the constant M , we could introduce instead a so-called 'cut-off operator' in the iterate equations (3.1a)-(3.1b), as seen in e.g. [40,41]. Note that in all numerical tests to be done in the next sections, we provide figures validating the second part of this assumption (cf. Section 5.4). With the above definitions, we state our main theoretical result.
if L u , L ϕ > 0, and if the model parameter ε > 0 is sufficiently large such that

9)
where ξ := (M λ max /λ min ) 2 > 0, and where c P , c K > 0 are the Poincaré and Korn constants, respectively, depending only on the domain B and spatial dimension d.
Proof. We begin by taking v = e i u and ψ = e i ϕ in (3.5a) and (3.5b), respectively, add the resulting equations together and obtain where we used the following inner product identity Discarding some non-negative terms from the left hand side of (3.10), using the fact that ess sup x∈B ϕ n (x) ≤ 1, in addition to the Lipschitz property of the degradation function g (2.6), yields where we also invoked the Assumption 1 in the last line, and applied the Cauchy-Schwarz inequality. Using the Young inequality, the properties of elastic tensor (2.2), and rearranging, leads to for some constants δ 1 , δ 2 > 0. Choosing (3.14) Next, by applying the Poincaré inequality on e i ϕ , and by applying successively the Poincaré and Korn inequalities on e i u , we obtain where c P , c K are the (squares of the) Poincaré and Korn constants, respectively (depending only on the domain B and spatial dimension d). Finally, employing these bounds on the left hand side of (3.14) yields Thus, for (3.16) to be a contraction estimate, ε must satisfy the following second order inequality Setting the left hand side of (3.17) equal to zero yields a second order polynomial, the discriminant of which must satisfy one of the following three statements: 1. If then P (ε) = 0 has two distinct positive real roots ε 1 , ε 2 > 0, in which case (3.16) is a contraction for ε ∈ (0, ε 1 ) ∪ (ε 2 , ∞).
Remark 3.1 (Convergence rate). According to the above proof, if the scheme is not converging for a given value of ε, then a larger or a smaller value may be chosen to rectify the situation. However, since crack surfaces become singular as ε → 0 (thus necessitating finer meshing, i.e., h → 0), we choose to state Theorem 3.1 with the condition that ε be large enough. We note also that due to some unknown constants in the convergence rate (3.8), it is not known whether this rate is optimal. Furthermore, working with a large ε is substantiated by the theory of phase field fracture being based on Γ convergence [2,3]. Applying this to phase field fracture was first done in [10]. Specifically, the setting is suitable when h = o(ε); namely when ε is sufficiently large.

Algorithm
In practice, we apply the stabilizations and penalizations proposed in the previous sections as outlined below. It is well-known (e.g., [37]) that the choice of γ is critical. If γ is too low, crack irreversibility will not be enforced.
On the other hand, if γ is too large, the linear equation system is ill-conditioned and influences the performance of the nonlinear solver. For this reason, γ is updated in at each iteration step. Better, in terms of robustness, is the augmention in such an iteration by an additional L 2 function Ξ, yielding a so-called augmented Lagrangian iteration going back to [19,22]. For phase-field fracture this idea was first applied in [42]. Thus, combining the staggered iteration for the solid and phase-field systems with the update of the penalization parameter Ξ yields the following algorithm: At the loading step t n . Choose initial Ξ 0 . Set γ > 0. repeat Iterate on i (augmented Lagrangian loop) Solve two-field problem, namely Solve elasticity in Problem (3.1a) Solve the nonlinear phase-field in Problem (3.1b) Set: (u n , ϕ n ) := (u n,i , ϕ n,i ). Increment t n → t n+1 .
For the stabilization parameters L u , L ϕ , we have the following requirements (somewhat similar to γ): If the stabilization is too small, the stabilization effects vanish. If the stabilization is too large, we revert to an unacceptably slow convergence, and potentially, may converge to a solution corresponding to an undesirable local minimum of the original problem. In order to deal with these issues, we employ here a simple, yet effective strategy: We draw L := L u = L ϕ from a range of suitable values and compare the results, i.e., L ∈ {1.0e−6, 1.0e−3, 1.0e−2, 1.0e−1}. Moreover, we include also the configurations L u = 0, L ϕ > 0 and L u = L ϕ = 0 in all the numerical tests to be done in the following.

Nonlinear solution, linear subsolvers and programming code
Both subproblems (phase field and mechanics) may be nonlinear. In our theory presented above, we assumed a standard elasticity tensor. However, the model (3.1a)-(3.1b) is too simple for most mechanical applications. More realistic phase-field fracture applications require a splitting of the stress tensor (based on an energy split) in order to account for fracture development only under tension, but not under compressive forces. Consequently, we follow here [31] and split σ into tensile σ + and compressive parts σ − : σ + := 2µ s e + + λ s < tr(e) > I, σ − := 2µ s (e − e + ) + λ s tr(e)− < tr(e) > I, where the elasticity tensor C has been replaced by the Lamé parameters, µ s and λ s . Moreover, I is the d × d identity matrix, and < · > is the positive part of a function. In particular, for d = 2, we have where λ 1 (u) and λ 2 (u) are the eigenvalues of the strain tensor e := e(u), and v 1 (u) and v 2 (u) the corresponding (normalized) eigenvectors. Finally, the matrix P is defined as P := P(u) := [v 1 |v 2 ]; namely, it consists of the column vectors v i , i = 1, 2. We notice that another frequently employed stress-splitting law was proposed in [4].
The modified scheme reads: • Step 1: given (u n,i−1 , ϕ n,i−1 , b n ) find u n,i such that Step 2: given (ϕ n,i−1 , u n,i , ϕ n−1 ) find ϕ n,i such that These modifications render the displacement system (4.2a) nonlinear, for which we use a Newton-type solver. The phase field equation is also nonlinear due to the penalization term and the stress splitting. Our version of Newton's method is based on a residual-based monotonicity criterion (e.g., [17]) outlined in [45][Section 3.2]. Inside Newton's method, the linear subsystems are solved with a direct solver; namely UMFPACK [16]. All numerical tests presented in Section 5 are implemented in the open-source finite element library deal.II [5,6]. Specifically, the code is based on a simple adaptation of the multiphysics template [43] in which specifically the previously mentioned Newton solver is implemented.

Numerical experiments
In this section, we present several numerical tests to substantiate our algorithmic developments.

Single edge notched tension test
This test was applied for instance in [31]. The configuration is displayed in Figure 1. We use the system (3.1a)-(3.1b). Specifically, we study our proposed iterative schemes on different mesh levels, denoted as refinement (Ref.    In particular, we see the initial crack build in the geometry in the right part where the domain has a true discontinuity. In the left part, the domain is cracked using the phase-field variable. Here, the displacement variable is still continuous since we are using C 0 finite elements for the spatial discretization.

Single edge notched shear test
The configuration of this second setting is very similar to Example 1 and was first proposed in a phase-field context in [31]. We now use the model with strain-energy split (4.2a)-(4.2b). The parameters and the geometry (see Figure  1) are the same as in the previous test case. The boundary condition is changed from tensile forces to a shear condition (see also again Figure 1): As quanitity of interest we evaluate F x in (5.2). Our findings are shown in the Figures 9, 10, 11, 12, 13, 14, and 15. The major difference to Example 1 is that the scheme is converging even with L u = 0, as computationally justified in Figure 11. As in Example 1, the load-displacement curves are close to the published literature and, again, the proposed L scheme is robust under mesh refinement (see Figures 12 -15).     The results are very comparable to the published literature. In particular, it is nowadays known that the proposed Miehe et al. stress splitting does not release all stresses once the specimen is broken (see [1]) and it is also known that we do not see convergence of the curves when both h and ε are refined (see [25]).

L-shaped panel
For the configuration of this third example we refer to [1,29,44], which are based on an experimental setup [46]. We use again the model with strain-energy split; namely (4.2a)-(4.2b). Moreover, in this test a carefully imposed irreversibility constraint is important since the specimen is pushed, pulled, and again pushed (see Figure 16 for the loading history on the small boundary part Γ u ). In the pulling phase the fracture vanishes if the penalization is not strong enough.
The geometry and boundary conditions are displayed in Figure 1. In contrast to the previous examples, no initial crack prescribed. The initial mesh is 1, 2 and 3 times uniformly refined, leading to 300, 1200, 4800 mesh elements, with h = 29.1548 mm, 14.577 mm, 7.289 mm, respectively.
We increase the displacement u D := u y = u y (t) on Γ u := {(x, y) ∈ B| 470 mm ≤ x ≤ 500 mm, y = 250 mm} over time, where Γ u is a section of 30 mm length on the right corner of the specimen. We apply a loading-dependent, non-homogeneous Dirichlet condition (see also Figure 16  We use µ s = 10.95 kN/mm 2 , λ s = 6.16 kN/mm 2 , and G c = 8.9 × 10 −5 kN/mm. The time (loading) step size is δt =10 −3 s. Furthermore, we set k = 10 −10 h[mm] and ε = 2h. As before, we observe the number of Newton iterations and we evaluate the surface load vector τ = (F x , F y ) := Γu σ(u)ν ds, with normal vector ν, and now we are particularly interested in F y . The crack path at the chosen time step snapshots in Figure 17 corresponds to the published literature [44,29,1]. The load-displacement curves and the number of iterations for different L and corresponding mesh refinement studies are displayed in the Figures 18, 19, 20, 21 and 22. Observe that stabilizing the mechanics subproblem has no effect in this example. At left, the load-displacement curves displaying the evolution of F y versus u y are shown. At right, the number of staggered iterations is displayed.

Verification of Assumption 1
In this last set of computations, we verify whether Assumption 1 holds true in our computations. We choose some prototype settings, namely on the coarest mesh level Ref. 4 and L u = L ϕ = 1e − 6. In Figure 23, we observe that ess sup x∈B |e(u n (x))| varies, but always can be bounded from above with M > 0. The value of ess sup x∈B |e(u n (x))| is the final strain when the L-scheme terminates. The minimum and maximum values shows that there are no significant variations in ess sup x∈B |e(u n (x))| during the L-scheme iterations with respect to the finally obtained value.

Conclusions
We have proposed a novel staggered iterative algorithm for brittle fracture phase field models. This algorithm is employing stabilization and linearization techniques known in the literature as the 'L-scheme', which is a generalization of the Fixed Stress Splitting algorithm coming from poroelasticity. Through theory and numerical examples we have investigated the performance of our proposed variants of the L-scheme for brittle fracture phase field problems. Under natural constraints that the elastic mechanical energy remains bounded, and that the model parameter ε is sufficiently large (i.e., that the diffusive zone around crack surfaces must be sufficiently thick), we have shown that a contraction of successive difference functions in energy norms can be obtained from the proposed scheme. This result implies the algorithm is converging monotonically with a linear convergence rate. However, in the convergence analysis there appears some unknown constants which makes the precise convergence rate, as well as the precise lower bound on ε unknown.
We provide detailed numerical tests where our proposed scheme is employed on several phase field brittle fracture bench-mark problems. For each numerical example we provide findings for different values of stabilization parameters. For most cases we let L u = L ϕ > 0, but for comparison we include also for the stabilization configurations L u = 0 with L ϕ > 0, and L u = L ϕ = 0. For the test cases presented here, there is only Example 1 where L u = 0 does not work. This might be due to the very rapid crack growth, which sets Example 1 apart from Examples 2 and 3. In this regard, we conclude that further work is needed to find an optimal configuration of L u and L ϕ . For all numerical test we also provide computational justification for the assumption of bounded elastic mechanical energy. Furthermore, a slight dependency on h in the iteration counts is observed in the numerical tests, but this is expected since we use ε = 2h, and as our analysis demonstrates, the convergence rate is dependent on ε. The variation in iteration numbers with mesh refinement is in any case sufficiently small enough that we conclude our algorithm is robust with respect to mesh refinement.
Moreover, due to the iteration spikes at the critical loading steps, we have included, for comparison, several results in which the iteration has been truncated (labeled LF I in Examples 1-3). Due to the monotonic convergence of the scheme, this strategy still produces acceptable results, while effectively avoiding the iteration spikes. We therefore conclude, at least for the particular examples presented here, that a truncation of the L-scheme can be employed for greatly improved efficiency with only negligible (depending on the situation at hand, of course) loss of accuracy.