title: An accelerated staggered scheme for phase-field modeling of brittle fracture
authors: Erlend Storvik Jakub Wiktor Both Juan Michael Sargado Jan Martin Nordbotten Florin Adrian Radu
text: 
  
  Reputation-driven Decision-making in
  Networks of Stochastic Agents
      
  David Maoujoud1
  Gavin Rens2
  
      December 30, 2023
  ==================================================================
  
  
  
  
  There is currently an increasing interest in developing efficient solvers for phase-field modeling of brittle fracture. The
  governing equations for this problem originate from a constrained minimization of a non-convex energy functional, and the most commonly used solver is a staggered solution scheme. This is known to be robust compared to the monolithic Newton method, however, the staggered scheme often requires many iterations to converge when cracks are evolving. The focus of our work is to accelerate the solver through a scheme that sequentially applies Anderson acceleration and over-relaxation, switching back and forth depending on the residual evolution, and thereby ensuring a decreasing tendency. The resulting scheme takes advantage of the complementary strengths of Anderson acceleration and over-relaxation to make a robust and accelerating method for this problem. The new method is applied as a post-processing technique to the increments of the solver, hence, the implementation can be done with minor modifications to already available software. Moreover, the cost of combining the two acceleration schemes is negligible. The robustness and efficiency of the method are demonstrated through numerical examples.
  
  
  
  
  § INTRODUCTION
  
  Mathematical modeling of brittle fracture propagation is an important and challenging topic in engineering sciences. The main difficulty arises in the transition between the distinct material properties in the fracture and the bulk domain. In this paper, we consider a variational phase-field model, as introduced by Bourdin, Francfort, and Marigo <cit.>. A smooth indicator that marks the broken and unbroken parts of the material regularizes the sharp crack topology. This enables modeling of fractures without conforming meshes or path-tracking algorithms (as in XFEM <cit.>). However, fine meshes are needed to resolve the regularized region between the fracture and the bulk domain. 
  
  
  
  
  
  The system is modeled by minimizing its energy as a function of material displacement and the indicator function. This leads to a system of coupled, nonlinear equations which is challenging to solve. The most common technique, due to its robust nature, is the staggered scheme. This method decouples the system and sequentially updates the displacement and indicator variable by solving their respective subproblems. However, the convergence properties are at times very bad, and iterating to satisfactory precision can result in large numbers of iterations <cit.>. The monolithic Newton method, on the other hand, does not show the same numerical robustness. Therefore, several attempts have been made to find a method that is both fast and robust. A monolithic, modified Newton method was proposed in <cit.>, a monolithic quasi-Newton method of BFGS type was applied in <cit.> and <cit.>, a monolithic line-search Newton method (dependent on the system energy) was applied in <cit.>, and the truncated nonsmooth Newton multigrid method was proposed in <cit.>. In <cit.>, the L-scheme <cit.> was applied in the context of an augmented Lagrangian solver, and a combination of an over-relaxed staggered scheme and the monolithic Newton method was applied in <cit.>. 
  
   In this paper, we propose a novel strategy to accelerate the classical staggered solution scheme solely utilizing two techniques for post-processing increments: Anderson acceleration and over-relaxation. In addition to accelerating the staggered scheme without sacrificing robustness, the new method allows the use of already available staggered scheme solvers with minor modifications to the implementation.
  
  
  
  
  
  Anderson acceleration was first developed in <cit.> for integral equations. Since then, it has seen many applications, including electronic structure computations <cit.> and flow in deformable porous media <cit.>. It is a multi-secant, quasi-Newton method that has been related to a preconditioned GMRES <cit.>. Moreover, the method post-processes the increments of the solver by approximating the inverse of the Jacobian of the system by reusing previous iterations. It can, therefore, easily be applied in combination with splitting techniques such as the staggered scheme while maintaining the decoupled nature of the scheme. 
  
  In <cit.>, the authors show theoretically that the Anderson acceleration improves the convergence rate of linearly convergent schemes, which is the case of the staggered scheme. However, as proved in <cit.>, the convergence is only local. In the case of phase-field modeling of brittle fracture, this is a challenge; when fractures are initiating or propagating, the system state usually jumps drastically between consecutive loading steps. Therefore, a “naive” application of Anderson acceleration is not suitable for this application, as will be demonstrated in the numerical examples of this work. Recently, there has been an increasing interest in modified Anderson acceleration methods to overcome issues of local convergence. In <cit.> a safeguard, based on the residual norm of the problem, is applied to restart Anderson acceleration, and in <cit.> a periodically restarted Anderson acceleration is applied within a Richardson fixed-point iteration to accelerate the convergence of iterative solvers for large sparse linear systems. 
  
  Relaxation was applied to the staggered scheme on a phase-field model of brittle fracture in <cit.>. It is a post-processing method that updates each iterate by relaxing (scaling) its increment. For the purpose of this work, over-relaxation (a scaling larger than one) is of particular interest. This is because the staggered scheme steadily moves towards the final configuration of each loading step, and over-relaxation might move the iterates further during each iteration, potentially accelerating the convergence. For the particular loading steps in which fractures are propagating, the gain can be quite substantial. There is, however, a drawback with over-relaxation: Near the solution of each loading step one might end up over- and undershooting the solution sequentially leading to poor performance.
  
  The most important observation of this paper is the complementary strengths of these two acceleration techniques; Anderson acceleration accelerates close to the solution, while over-relaxation accelerates during loading steps with large jumps in the solution (e.g., during crack propagation). We propose an acceleration algorithm that switches between Anderson acceleration and over-relaxation during each loading-step ensuring convergence at an accelerated rate. This scheme is related to the one in <cit.> where the authors switch between over-relaxation and monolithic Newton. However, for the new acceleration scheme, proposed in this paper, both of the combined acceleration methods function as post-processes to the increments of the standard staggered scheme. In other words, the new acceleration method can be implemented with minor modifications to already available software. Moreover, switching between the two acceleration techniques does not change the sparsity of the underlying linear systems. The switch criterion is based on the history of the residual norms of the staggered solution steps.  
  
  To summarize, the main contributions in this paper are:
  
   
    * Presentation of the difficulties encountered with the application of plain Anderson acceleration and over-relaxation applied to the staggered solution scheme for variational phase-field modeling of brittle fracture. 
   
    * A new acceleration algorithm that exploits the complementary strengths of Anderson acceleration and over-relaxation, utilizing residual norm evolution as a rule for switching between the methods.
   
    * The performance of the proposed acceleration scheme is demonstrated through thorough numerical examples including classical benchmark problems. 
  
  
  The paper is structured as follows: The mathematical model and numerical discretization are presented in Section <ref>. Here, we introduce the energy functional which is subject to minimization together with the discretization. In Section <ref>, the staggered scheme and the acceleration techniques are presented. Both Anderson acceleration and relaxation are described before the combined acceleration scheme is presented together with the inexact Newton modification. Section <ref> contains the numerical study of the accelerations applied to the staggered scheme. We test the staggered scheme both with and without the combinations of Anderson acceleration and relaxation. Moreover, the optimal depth of Anderson acceleration and the choice of relaxation parameter is discussed. Finally, some concluding remarks are made in Section <ref>.
  
  
  
  
  
  
  § MATHEMATICAL PROBLEM
  
  In this section, the mathematical problem that is considered throughout the paper is presented. An elastic medium, represented by the domain Ω⊂ℝ^d with d= 2 or 3, is subject to loading through traction forces, t, along Γ_N and displacement, u_D, along Γ_D to the extent that it might break. Here, Γ_N∪Γ_D=∂Ω are subsets of the boundary of the domain, Ω, and Γ_D has nonzero measure. The state of the material is modeled by Griffith's criterion <cit.>, with constant G_c, and a smooth indicator function (the phase-field variable) φ:Ω→ [0,1] describes the state of the damage to the material. The phase-field is defined to take the value 0 whenever the material is unbroken, and 1 when the material is broken, and a model parameter ℓ determines the width of the regularized zone where the phase-field transitions from 0 to 1.
  
  
  
   §.§ The energy of the system
  
  Following the work of <cit.>, we can express the total energy of the system as a sum of the medium's elastic energy, the surface energy dissipation associated with the broken parts of the material and external work related to traction. 
  Now, let u denote the material displacement and define the total energy functional as 
  
      ℰ( u, φ):=∫_Ωℰ_c(φ)+ ℰ_m( u,φ)    d x-∫_Γ_N t · u  ds
  
  where 
  
      ℰ_c(φ) := G_c/2(φ^2/ℓ+ℓ∇φ·∇φ),
  
  and 
  
      ℰ_m( u, φ) := g(φ)Ψ^+(ε)+ Ψ^-(ε) - b· u.
  
  Here, we have applied the degradation function 
      g(φ):=(1-κ)(1-φ)^2+κ,
   where κ is a “small” constant. Other choices have been proposed in <cit.>. Moreover, the material is assumed to be homogeneous and isotropic, and the elastic strain energy functional 
      Ψ(ε) :=1/2ε:ℂ:ε =μ(ε:ε)+λtr(ε)^2/2,
  
  where ε = ∇ u+∇ u^⊤/2 is the linearized elastic strain tensor and μ and λ are the Lamé parameters, has been decomposed into “tensile”, Ψ^+, and “compressive”, Ψ^-, parts. The additive spectral decomposition
  
      Ψ^±(ε):=μ(ε_±:ε_±)+λ⟨tr(ε)⟩_±^2/2,
  
  proposed in <cit.> has been employed.
  Here, ⟨ a ⟩_± :=12(a±|a|) and ε_± :=∑_i⟨ε_i⟩_± n_i ⊗ n_i where {ε_i} and { n_i} are the principal strain and principal strain directions, respectively. Additionally, the material is unable to heal, and the constraint ∂_t φ≥ 0 is applied accordingly.
  
  
  
   §.§ Time discretized, contiuous-in-space equations
  
  The loading procedure is discretized by the implicit Euler scheme, giving the non-healing constraint at loading step n≥1: 
  
      φ^n(x)-φ^n-1(x) ≥ 0   ∀ x∈Ω .
  
  
  
  
  
  Now, we define the displacement solution space
  V^n = { v ∈(H^1(Ω))^d |   v|_Γ_D= u_D^n}, the displacement test space V^0={ v ∈(H^1(Ω))^d |   v|_Γ_D=0}, and the phase-field solution and test space 
  Q = H^1(Ω).
  Then, the solution ( u^n, φ^n)∈ V^n × Q at loading step n≥1 is given by 
  
      ( u^n,φ^n):=_ u, φ{ℰ( u, φ,  t^n)  |  u ∈ V^n,φ∈ Q}.
  
  Letting ⟨·,·⟩_X denote the usual L^2 inner product over the domain X and denoting 
      σ^±( u ) := ∂Ψ^±(ε( u))∂ε( u),
   we find the variation of the energy (<ref>) with respect to u and φ respectively: 
  
      ℰ_δ u( u,φ, v)    =   ⟨(g(φ)σ^+( u)+σ^-( u)), ε( v)⟩_Ω -⟨ b, v⟩_Ω-⟨ t, v⟩_Γ_N
      ℰ_δφ( u, φ, q)    =   ⟨ g'(φ)Ψ^+(ε),q⟩_Ω + G_c/ℓ(⟨φ,q⟩_Ω +ℓ^2⟨∇φ,∇ q⟩_Ω).
  
  It is now easy to see that the solution to (<ref>), ( u^n, φ^n), satisfies the system of equations
  
      ℰ_δ u( u^n,φ^n, v)    =    0
      ℰ_δφ( u^n, φ^n, q)    =    0
  
  for all v ∈ V^0 and q∈ Q.
  
  The inequality (<ref>) still requires some special treatment, and in this paper, we follow the approach of <cit.> and introduce a history variable; 
  
      ℋ^n := max_k≤ nΨ^+(ε( u^k) ).
  
  A modified version of the variation with respect to φ is defined by
  
      ℰ_δφ( u, φ, q) =  ⟨ g'(φ)ℋ^n,q⟩_Ω + G_c/ℓ(⟨φ,q⟩_Ω +ℓ^2⟨∇φ,∇ q⟩_Ω),
  
  and the solution at loading step n will be defined as the pair ( u^n,φ^n)∈ V^n× Q that satisfies 
  
      ℰ_δ u( u^n,φ^n, v)    =    0
      ℰ_δφ( u^n, φ^n, q)    =    0
  
  for all v∈ V^0 and q ∈ Q.
  
  
   §.§ Spatial discretization
  
  To solve (<ref>)–(<ref>) we apply conforming linear finite elements <cit.>, both for the phase-field variable and the displacement. Let 𝒯^h_Ω={T_k}_k be a decomposition of the domain Ω into simplices, T_k, and define, at loading step n≥1, the spaces 
  
      V_h^n    =   { v_h ∈ (H^1(Ω))^d  |  v_h|_T_k∈ (𝒫^1(T_k))^d ∀  T_k ∈𝒯_Ω,   v_h|_Γ_D= u_D^n },
      
       Q_h    =   {q_h ∈ H^1(Ω)  |   q_h|_T_k∈𝒫^1(T_k) ∀  T_k ∈𝒯_Ω},
  
  and V_h^0 accordingly with zero trace.
  The system of equations to be solved is then: Find ( u^n_h,φ^n_h)∈ V^n_h × Q_h such that 
  
      ℰ_δ u( u^n_h,φ^n_h, v_h)    =    0
      ℰ_δφ( u^n_h, φ^n_h, q_h)    =    0
  
  for all v_h ∈ V^0_h and  q_h ∈ Q_h. Following standard procedures, this naturally translates to the algebraic residual equations 
      Res_ u( u_h^n,φ_h^n)   =   0
      Res_φ( u_h^n, φ_h^n)   =   0,
  
  where Res_ u and Res_φ denote the algebraic residuals corresponding to (<ref>) and (<ref>), respectively.
  
  
  
  
  § STAGGERED SCHEME AND ACCELERATION
  
  The discrete governing equations (<ref>)–(<ref>) are strongly nonlinear and coupled. In this paper, we apply the staggered scheme <cit.> to solve them, decoupling the equations. We let i≥ 1 be the iteration index and define the staggered scheme as: Given φ^n,i-1_h∈ Q_h, find ( u^n,i_h,φ^n,i_h)∈ V^n_h × Q_h such that 
  
      ℰ_δ u( u^n,i_h,φ^n,i-1_h, v_h)    =    0
      ℰ_δφ( u^n,i_h, φ^n,i_h, q_h)    =    0
  
  for all ( v_h,q_h)∈ V^0_h × Q_h and φ^n,0_h := φ^n-1_h. The iterations are terminated when the following stopping criterions are reached:
  
      Res_ u( u_h^n,i,φ_h^n,i)_2   ≤   Tol_Res,Abs,
      Res_ u( u_h^n,i,φ_h^n,i)_2/Res_ u( u_h^n,1,φ_h^n,0)_2   ≤   Tol_Res,Rel, 
       u_h^n,i- u_h^n,i-1_L^2(Ω)+φ_h^n,i-φ_h^n,i-1_L^2(Ω)   ≤   Tol_Inc,Abs,
       u_h^n,i- u_h^n,i-1_L^2(Ω)/ u_h^n,1_L^2(Ω)+φ_h^n,i-φ_h^n,i-1_L^2(Ω)/φ_h^n,0_L^2(Ω)   ≤   Tol_Inc,Res,
  
  for given tolerances Tol_Res,Abs, Tol_Res,Rel, Tol_Inc,Abs and Tol_Inc,Rel.
  Notice that controlling the residuals corresponding to the phase-field equation (<ref>) is redundant due to it being solved second in the staggered scheme by an exact linear solver. 
  
  To solve the nonlinear equation (<ref>) we apply the Newton method with the relative stopping criterion 
  
      Res_ u( u_h^n,i,j,φ_h^n,i-1)_2/Res_ u( u_h^n,1,φ_h^n,0)_2≤Tol_inner.
  
  Here, j≥ 1 is the iteration index for the Newton method and the initial guess is chosen as the previous staggered iteration u_h^n,i,0:= u_h^n,i-1.                                                                                                                                                                                                                                                                                                         
  
  
  
  
  The staggered scheme (<ref>)–(<ref>) is closely related to the alternate minimization method (it differs in the application of the history variable (<ref>)) and is known to be a robust solution method <cit.>. However, it might require a large number of iterations to reach satisfactory tolerances <cit.>. We aim to accelerate this slow convergence and propose a combination of Anderson acceleration and over-relaxation. We note that the staggered solution scheme can be written as the fixed-point iteration
  
      x^n,i_h := 𝒮( x^n,i-1_h) =  x^n,i-1_h+Δ𝒮( x_h^n,i-1)
  
  where 𝒮 is the staggered solution scheme operator, Δ𝒮 is the increment of the staggered scheme and x_h^n,i is the vector [ u^n,i_h; φ^n,i_h ].
  
  Now, we present both the Anderson acceleration and the relaxed staggered scheme and describe their strengths and weaknesses. Then, taking advantage of the strengths of both schemes, a combined scheme is presented. 
  
  
  
   §.§ Anderson acceleration
  
  Anderson acceleration is a multi-secant method that mimics the monolithic Newton method. The acceleration acts as a post-processing procedure that updates the current iterate by a linear combination of the m previous iterates, according to their respective increments. The value of m is free to be chosen and is known as the depth of the acceleration. Moreover, it can be chosen adaptively. At loading step n, the Anderson accelerated staggered scheme of depth m reads:
  
  
  
  Algorithm <ref> is independent of the underlying fixed-point iteration, but is presented for the application to the staggered scheme here. An important feature of Anderson acceleration is that it preserves the decoupled nature of the staggered scheme, hence, the subproblem solvers are unaffected by it.
  
  It is demonstrated in the numerical section that Anderson acceleration improves the convergence when close to the solution. However, the acceleration might deteriorate otherwise. This is especially important to notice for brutal crack propagation, where it sometimes fails to converge at all.
  
  
  
  
   §.§ Over-relaxation
  
  Relaxation applied to each subproblem of the staggered solution scheme was described and applied in <cit.>. The method first calculates the increment Δ u^n,i-1_h obtained by solving equation (<ref>),  before defining the updated iterate as 
      u_h^n,i:= u^n,i-1_h +ωΔ u^n,i-1_h,
  
  where ω∈(0,2) is a parameter. This new iterate u^n,i_h is now passed on to equation (<ref>) and the same procedure is executed for the phase-field resulting in the updated iterate φ^n,i_h. Following standard literature on iterative methods, we refer to the choice ω∈(1,2) as over-relaxation and ω∈(0,1) as under-relaxation. At the n-th loading step the relaxed staggered scheme reads:
  
  
  
  
  Under-relaxation is robust when applied to the staggered scheme, however, it usually slows down the scheme. Over-relaxation, on the other hand, tends to accelerate the loading steps of the staggered solution scheme where cracks occur, while it might slow down the process for quasi-static loading steps.
  
  
  
   §.§ Combining Anderson acceleration and over-relaxation
  
  As neither Anderson acceleration nor over-relaxation should be applied naively to the staggered scheme (<ref>)–(<ref>), due to their mentioned weaknesses, we propose a combined robust acceleration scheme. The key observations that motivate such a method are:
  
   
    * Anderson acceleration is locally accelerating, while over-relaxation might struggle close to the solution.
   
    * Anderson acceleration is applied as a post-processing algorithm to the increments of the staggered scheme, hence, only trivial implementation is required to switch between relaxation and Anderson acceleration. 
   
    * During crack propagation the residuals for the staggered scheme show a stagnating, oscillatory behavior, and during quasi-static steps they are strictly decreasing, see <cit.> and Figure <ref>. Therefore, it is possible to use residual evolution as a rule for switching between the acceleration techniques. 
  
  
  
  
   A new parameter N_ω→ AA∈ℕ, related to the switch from relaxation to Anderson acceleration is defined, and at loading step n the combined accelerated staggered scheme reads: 
  
   
    * Apply with Anderson acceleration of given depth m.
   
    * While the norms of the residuals are strictly decreasing, continue with Anderson acceleration until convergence.
   
    * If the norms of the residuals are not strictly decreasing, switch to relaxation with given parameter ω.
   
    * When the norms of the N_ω→ AA previous residuals are strictly decreasing go back to 1, and restart[Restart means that Anderson acceleration should be applied as it is in the first iteration, i.e., using no information of previous increments and iterates.] Anderson acceleration.
  
  
  
  
  Below, we give a pseudo-code for the new combined acceleration method. Define the residual norm Res_i = Res_ u( u_h^n,i,φ_h^n,i)_2 as in (<ref>), and notice that the application of Anderson acceleration and relaxation in the pseudo-code denotes the i-th step of the accelerations (see Algorithm <ref> and Algorithm <ref>).
  
  
  
  
  
  
  § NUMERICAL EXAMPLES
  
  This section explores the effects of the proposed acceleration methods from Section <ref> applied to the staggered scheme (<ref>)–(<ref>). Both Anderson acceleration and over-relaxation alone are shown to be infeasible acceleration methods when plainly applied to the staggered scheme, while the combined scheme is superior to the unaccelerated scheme for all tests.
  We consider four different test cases which are widely used for numerical studies in the literature:
  
   
    * A domain with a single notch subject to
   
    
    * tensile load;
    
    * shear load.
   
   
    * An L-shaped domain subject to loading.
   
    * Bending of an asymmetrically notched beam with holes.
  
  
  All the numerical examples have been implemented using modules from the DUNE project <cit.>, specifically dune-functions <cit.>.
  
  
  
  
  
   §.§ Single notch test
  
  Two of the most commonly found test cases in the literature are both based on the same single notch geometry <cit.>. They consist of a square domain with a pre-existing crack that penetrates half the domain, see Figure <ref>. The domain is held still at the bottom, and a displacement driven load is applied at the top boundary.
  
  
  
  
  
  
    §.§.§ Single notch tensile test case
  
  A tensile load is applied on the top boundary, and at loading step n we have 
  
      u^n_|_Γ^Top=[   0; u̅n ],
   where the load size u̅ is given in Table <ref> and Γ^Top is the top part of the boundary in Figure <ref>. Due to the load being strictly tensile, there is no need to split the elastic strain energy functional (<ref>) into tensile and compressive parts, which would effectively add nonlinearities to the system. Therefore, the first term in (<ref>) is replaced by ⟨ g(φ)σ( u):ε( u))⟩, 
  for 
      σ( u ) := ∂ψ(ε( u))∂ε( u).
   Material parameter values are chosen as in e.g., <cit.>, and can be found in Table <ref>. We employ a triangular mesh, which has been locally refined in the region where the crack is expected to propagate, see Figure <ref>.
  
  
   
  In this test case, the crack fully propagates in one single critical loading step, see Figure <ref>, in which the crack gradually expands through the domain with increasing staggered iteration count. Figure <ref> shows that, as expected, the staggered scheme under Anderson acceleration alone struggles as a consequence of its local convergence. The combined scheme, however, takes advantage of over-relaxation and its ability to move further each iteration and accelerates this particular loading step significantly. For the rest of the loading steps, the combined scheme accelerates by Anderson acceleration, as its local convergence is sufficient. The total number of iterations for the combined scheme is, therefore, smaller than those of the unaccelerated staggered scheme and the Anderson accelerated staggered scheme.  
  
  The traction vector is defined by 
  
      τ = (τ_x,τ_y) = ∫_Γ^Topσ·ν dS,
   
  where ν is the outward pointing normal vector, and σ is defined in (<ref>). For this problem, the load in the direction of interest is τ_y, and we observe in Figure <ref> that the load-displacement curves remain unchanged after the combined acceleration. This is an important observation that demonstrates that the acceleration method only affects the convergence properties of the solver, not the quality of the solution. The Anderson accelerated staggered scheme, however, does not converge in the maximal prescribed iterations for each loading step and we observe that its load-displacement curve is affected. 
  
  The total number of iterations is displayed in Figure <ref>, and there are several key observations. First of all, the combined scheme accelerates by more than 50 % for large relaxation parameters. We also observe that the depth of Anderson acceleration is not influential as long as it is larger than one. Moreover, there is a trend that more aggressive over-relaxation (higher ω) results in faster computations. Additionally, the combined acceleration scheme is robust with respect to the tuning parameters, m and ω, and exhibits convergence for all tested combinations.
  
  
  
  
  
  
  
  
  
  
  
  Notice that the plots of the number of iterations for depth m= 0 in Figures <ref>, <ref>, <ref> and <ref> are not corresponding to plain relaxation. Here, relaxation is switched on and off depending on residual evolution, turning it into a safeguarded relaxation. The same goes for the plots of plain Anderson acceleration with over-relaxation parameter ω=1. These correspond to safeguarded Anderson accelerations, similar to those that are proposed in <cit.>.
  
  
  
  
  
    §.§.§ Single notch shear case
  
  In this test case, a shear load is applied on the top boundary of a unit square domain with a prescribed crack that halfway penetrates the domain. The displacement boundary condition
  
      u^n_|_Γ^Top=[ u̅n;   0 ]
   is applied at loading step n. The load size u̅ is presented in Table <ref>, and the top part of the boundary Γ^Top is displayed together with more details on the domain and boundary conditions in Figure <ref>.
  The material properties are taken from <cit.> and displayed in Table <ref>. A triangular mesh, which is refined according to where the crack is expected to propagate, has been employed, see Figure <ref>. 
  
  Contrary to the tensile test case, the crack propagation happens gradually over the course of many loading steps, see Figure <ref>. Therefore, solutions at subsequent loading steps do not differ as significantly as for the brutal crack growth in the tensile test case. We expect that the Anderson acceleration is a more suitable choice for accelerating the staggered scheme. Indeed, Figure <ref> shows that even with the naive Anderson acceleration the staggered scheme is quite significantly accelerated. Moreover, the combined scheme is even better, and we see that it reaches convergence in every single loading step. 
  
  In the load-displacement curves, Figure <ref>, the load τ_x from (<ref>) is displayed for each loading step. The plot shows minor differences towards the end of the displacement. This is due to the scheme not converging in its given maximal iterations per loading step (see Table <ref>) for both the unaccelerated staggered scheme and the plain Anderson accelerated scheme in all the loading steps. This is similar to the tensile case where the load-displacement curves (Figure <ref>) also are affected.
  
  In Figure <ref>, we see the total iteration count for several acceleration depths in combination with over-relaxation. The figure shows that the staggered scheme is accelerated significantly for all combinations of Anderson acceleration and over-relaxation as long as the depth is greater than one. In fact, we have more than 80 % reduction in the total number of iterations when choosing a high relaxation parameter. Moreover, for this test case the plain Anderson acceleration is in itself a suitable alternative to the unaccelerated staggered scheme. Notice the difference between plain Anderson acceleration and Anderson acceleration combined with over-relaxation of depth one described in Remark <ref>.
  
  
  
  
  
   
  
  
  
  
  
   §.§ L-shaped domain subject to loading
  
  An L-shaped domain with a displacement boundary condition applied on the right part of the boundary is considered, see Figure <ref> for details. The displacement is uniformly increased on the boundary segment over 800 loading steps. As a result, a crack occurs in the inner corner, propagating into the domain, see Figure <ref>. A uniform quadrilateral mesh with a mesh diameter of 125/32 mm is employed. See Table <ref> for material and computational parameters.
  
  Here, the crack propagation has a character somewhere between the single notch tensile test and the single notch shear test. Crack initiation shows similar behavior as brutal crack propagation, but not as extreme as for the single notch tensile test. A large peak in the number of iterations is experienced when the crack initiates, see Figure <ref>. We observe that both the combined scheme and plain Anderson acceleration with depth m=1 accelerates for all loading steps. Moreover, the only difference between the combined acceleration and Anderson acceleration is in the large peak where the combined scheme outperforms Anderson acceleration. For the rest of the simulation, the staggered scheme converges in relatively few iterations (less than 30 per loading step), but the accelerated method converges faster in almost every loading step. 
  
  Figure <ref> displays the total number of iterations for plain over-relaxation with several relaxation parameters. A parabolic dependence on the relaxation parameter is observed, and choosing it to be too large results in more than three times the number of iterations that are required by the unaccelerated staggered scheme. Therefore, a plain application of over-relaxation is not recommended. The total number of iterations required by the combined acceleration, however, is significantly smaller than those of the unaccelerated staggered scheme, as observed in Figure <ref>. Although the reduction in the number of iterations is not as extreme as for the single notch shear test case the combined scheme accelerates robustly with respect to the tuning parameters. It is clear that any combination of Anderson acceleration and over-relaxation is superior to the unaccelerated staggered scheme, accelerating by approximately 40 %.
  
  The load-displacement curve for this test case is displayed in Figure <ref>. Here, the traction vector (see equation (<ref>)) is calculated on the bottom boundary and the vertical component τ_y is considered. We observe that, as all acceleration schemes converge within each loading step, the curves are completely overlapping. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
   §.§ Asymmetrical bending test 
  
  This test case considers a rectangular domain with three holes, slightly to the left, and a notch in the lower left part of the domain. It is subject to symmetrical displacement loading on the top boundary,
  
      u_h^n|_Γ^Top = [    0; u̅ n ].
  
  The beam is simply supported as shown in Figure <ref>. See Figure <ref> for details on boundary conditions and domain. Experimental results from <cit.> have shown that the crack path should hit the second hole, and we see from the numerical solution, Figure <ref>, that this also happens here. The mesh has been refined in the region where the crack is expected to propagate, see Figure <ref>. The problem parameters are chosen similarily to <cit.>, and are presented in Table <ref>.
  
  Here, we have two “critical” loading steps, in which the crack evolves and a large number of iterations is required, see Figure <ref>. For these loading steps, we see that the plain Anderson acceleration does not accelerate, while the combined acceleration performs very well. 
  
  In Figure <ref> the total number of iterations for the plain Anderson acceleration is displayed for several depths. We clearly observe that the staggered scheme is significantly decelerated for depths larger than one. In other words, Anderson acceleration is not a robust method in itself for this problem. The combined scheme, on the other hand, reduces the total number of iterations for all combinations of over-relaxation and Anderson acceleration, see Figure <ref>. There is, however, a tendency that larger relaxation parameters accelerate more, which is expected due to the brutal nature of the crack propagation in loading step 12, see Figure <ref>.  
  
  The traction vector (<ref>) is here calculated on the top boundary and the component of interest is τ_y. In Figure <ref>, the load-displacement curve is displayed, and the displacement is calculated at the left corner of the top boundary. They are, as expected, overlapping as there are no loading steps for any configurations in which the convergence is not achieved in the given maximal amount of iterations. 
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSION
  
  The staggered solution scheme is, due to its robustness, a popular method for solving variational phase-field models of brittle fracture. As it often requires a large number of iterations to converge we have proposed a method to accelerate it that exploits the complementary advantages of Anderson acceleration and over-relaxation. The acceleration method alternates between Anderson acceleration and over-relaxation according to a switch that depends on the norms of the previous residuals of the scheme. For problems without brutal crack growth, Anderson acceleration is quite efficient. It is, however, unstable for problems with brutal crack growth, and therefore, not a technique that can be applied without modifications. Over-relaxation, on the other hand, works well within regimes of brutal crack propagation, but might struggle when the iterates get close to the solution. The scheme shows robustness with respect to the tuning parameters, Anderson acceleration depth and relaxation parameter, and converges for all combinations. Moreover, there is a tendency that Anderson acceleration depths larger than one are insignificant, and that over-relaxation with parameters of at least 1.6 are the best choices. Therefore, we propose to apply the method with depth one and over-relaxation 1.6, although one might gain some speed in tuning these parameters to specific problems. 
  
  unsrt
  
  
  