Dynamic and weighted stabilizations of the L-scheme applied to a phase-field model for fracture propagation

We consider a phase-field fracture propagation model, which consists of two (nonlinear) coupled partial differential equations. The first equation describes the displacement evolution, and the second is a smoothed indicator variable, describing the crack position. We propose an iterative scheme, the so-called L-scheme, with a dynamic update of the stabilization parameters during the iterations. Our algorithmic improvements are substantiated with two numerical tests. The dynamic adjustments of the stabilization parameters lead to a significant reduction of iteration numbers in comparison to constant stabilization values.



There are no comments yet.


page 1

page 2

page 3

page 4


Tusas: A fully implicit parallel approach for coupled nonlinear equations

We develop a fully-coupled, fully-implicit approach for phase-field mode...

A convergent numerical scheme for a model of liquid crystal dynamics subjected to an electric field

We present a convergent and constraint-preserving numerical discretizati...

An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on LR NURBS

We present an adaptive space-time phase field formulation for dynamic fr...

Variable-Order Fracture Mechanics and its Application to Dynamic Fracture

This study presents the formulation, the numerical solution, and the val...

Efficient Solvers for Nonstandard Models for Flow and Transport in Unsaturated Porous Media

We study several iterative methods for fully coupled flow and reactive t...

Time discretization of a nonlocal phase-field system with inertial term

Time discretizations of phase-field systems have been studied. For examp...

Backprop Evolution

The back-propagation algorithm is the cornerstone of deep learning. Desp...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

This work is an extension of BrWiBeNoRa19 in which an L-type iterative scheme (see list2016study ; MR2079503 ) with stabilizing parameters for solving phase-field fracture problems was proposed. In BrWiBeNoRa19 , the stabilization parameters were chosen as constants throughout an entire computation. With these choices, the convergence of the scheme has been proven rigorously. The resulting approach performs well in the sense that an unlimited number of iterations compared to a truncated scheme yields the same numerical solution. The results were validated by investigating the load-displacements curves. Moreover, the robustness of the scheme w.r.t. spatial mesh refinement was shown. Nonetheless, the iteration numbers (for an unlimited number of iterations) remained high.

In this work, we propose and compare two extensions of the aforementioned scheme. First, we update the scheme parameters dynamically. Second, we use an adaptive weight depending on the fracture location inside the domain. For the latter idea, we use the phase-field variable to weight locally.

The outline of this work is as follows: In Section 2 the model is stated whereas Section 3 presents the dynamic choice of the stablization parameters. In Section 4, we present two numerical tests to study the performance of the proposed scheme.

2 The phase-field fracture model

We consider the crack propagation model proposed in BrWiBeNoRa19 . Let be a -dimensional, polygonal and bounded domain and a maximal time. We use common notations, and in particular denotes the space of functions on having essentially bounded weak derivatives in any direction, while contains the functions having square integrable weak derivatives, and vanishing at the boundary of (in the sense of traces). For the ease of writing, we let and . For (almost every) location and time

the vector-valued displacements are denoted by

. The fracture and its propagation within are modeled with the help of a phase field variable

, which approximates the characteristic function of the intact region of

. Written in weak form, the fracture propagation model in resumes to finding such that for : &• Step 1: given find such that
&a_u(u^n,i,v) := L_u(u^n,i - u^n,i-1,v) + ( g(φ^n,i-1) σ^+(u^n,i), e(v) ) + ( σ^-(u^n,i), e(v) ) = 0,
& ∀v ∈V_h,
&• Step 2: given find such that
&a_φ(φ^n,i,ψ) := L_φ(φ^n,i - φ^n,i-1,ψ) + G_c ε(∇φ^n,i, ∇ψ) - Gcε (1 - φ^n,i, ψ)
&    + (1-κ)(φ^n,i σ^+(u^n,i):e(u^n,i),ψ) + ( Ξ+ γ[φ^n,i - φ^n-1]^+, ψ)= 0,
& ∀ψ∈W_h. Here, we note that the ‘time’ appears only the irreversibility constraint , yielding an incremental problem and which is regularized using an augmented Lagrangian penalization as proposed in WheWiWo14 . Here, is an function and a positive parameter.

Furthermore, in the above, is a (small) phase-field regularization parameter, is the critical elastic energy restitution rate, and is a regularization parameter used to avoid the degeneracy of the elastic energy. The latter is similar to replacing the fracture with a softer material. Next, is the degradation function, and

is the strain tensor.

The stress tensor in the above is split into a tensile and compressive part,

where stands for the positive cut of the argument. Further, with

being the matrix containing the unit eigenvectors corresponding to the eigenvalues of the strain tensor

. In particular, for one has and

3 The -scheme with dynamic updates of the stabilization parameters

The iteration iter2_mod-iter1_mod is essentially the scheme proposed in BrWiBeNoRa19 , in which the stabilization parameters and are taken constant. To improve the convergence behaviour of the scheme, we propose a dynamic update of these parameters.

Dynamic update at each iteration / constant in space:  The iteration discussed in BrWiBeNoRa19 uses constant parameters and . With this choice, the convergence has been proved rigorously. However, the number of iterations can remain high. High iteration numbers for phase-field fracture problems were also reported in GeLo16 ; Wi17_SISC . To improve the efficiency, we suggest in this work to update and at each iteration :

Inspired by numerical continuation methods in e.g. AllGe90 , one would naturally choose a large and to obtain a decreasing sequence , updated until a lower bound is reached. However, this seems not to be a good choice in phase-field fracture since the system does not have a unique solution. Consequently, with increasing the iterations would oscillate in approaching one or another solution, and the algorithm convergence deteriorates. For this reason, we propose the other way around: the closer the iteration is to some solution, the larger the stabilization parameters is chosen, so that the iterations remain close to this solution. We choose , yielding up to a maximal .

On the specific choice of the parameters:  A possible choice for is (), while

. This heuristic choice and may be improved by using the solution within the iteration procedure, or a-posteriori error estimates for the iteration error. Moreover,

is motivated as follows. Higher values greater than would emphasize too much the stablization. On the other hand, too low values, do not lead to any significant enhancement of the convergence behaviour. We substantiate these claims by also using and in our computations.

Dynamic update using the iteration:  An extension of the strategy is to adapt the -scheme parameters in space by using the phase-field variable . We still take , but now . Away from the fracture, we have and essentially only the elasticity component iter1_mod is being solved. On the other hand, the stabilization is important in the fracture region, for which we take

Recalling that the fracture is characterised by , it becomes clear that the stabilization parameters are acting mainly in the fracture region. Finally, to improve further the convergence behaviour of the scheme we adapt at each iteration. In this case we take

At the loading step
Choose , , as well as and . Set .
     Let ;
     Solve the two problems, namely
          Solve the (nonlinear) elasticity iter2_mod
          Solve the nonlinear phase-field iter1_mod
Set   .
Increment   .
Algorithm 1 Dynamic variant of the L-scheme for a phase-field fracture

The final algorithm:  The algorithm is based on the iterative procedure for phase-field fracture originally proposed in WheWiWo14 . Therein, the inequality constraint is realized by an augmented Lagrangian iteration. Within this loop we update the scheme parameters too. The resulting is sketched in Algorithm 1, in which is taken, and .

Remark 1

For the solution of both nonlinear subproblems iter2_mod and iter1_mod, we use a monotonicity-based Newton method (details see e.g., in Wi17_SISC ) with the tolerance . Inside Newton’s method, we solve the linear systems with a direct solver.

4 Numerical tests

We consider two test examples. Details for the first test van be found in MieWelHof10a . The setup of the second test can be found for instance in MesBouKhon15 . Both examples were already computed in BrWiBeNoRa19 and the results therein are compared to the ones obtained here. The scheme is implemented in a code based on the deal.II library dealII91 .

Figure 1: Examples 1 and 2. The following conditions are prescribed: on the left and right boundaries, = and traction-free in -direction. On the bottom part, = . On , = and is as stated in diri_ex_2. Finally, the lower part of the slit is fixed in -direction, i.e., = . Right: Asymmetric notched three point bending test. The three holes have each a diameter of . All units are in .
Figure 2: Examples 1 and 2. Numerical solutions on the finest meshes and at the end time. The cracks are displayed in dark blue color.

Single edge notched shear test:  The configuration is shown in Figure 1. Specifically, we use = , = , and = . The crack growth is driven by a non-homogeneous Dirichlet condition for the displacement field on , the top boundary of . We increase the displacement on over time, namely we apply non-homogeneous Dirichlet conditions: u_x &= t ¯u,  ¯u = , where denotes the current loading time. Furthermore, we set [mm] and [mm]. We evaluate the surface load vector on the as


with normal vector , and we are particularly interested in the shear force . Three different meshes with (Ref. 4), (Ref. 5) and (Ref. 6) elements are observed in order to show the robustness of the proposed schemes. The results are shown in Figure 6.

Our findings are summarized in Figure 3. The numerical solutions for all four different strategies for choosing are practically identical, only the number of iterations being different. Here, and denote tests in which are taken constant throughout the entire computation. The newly proposed dynamic versions are denoted by L dynamic and L dyn. weighted. We observe a significant reduction in the computational cost when using the dynamic -schemes. The maximum number of iterations is for both the weighted version and the spatially-constant -scheme. This number is reduced to iterations using while the accuracy only slightly changes.

Figure 3: Example 1. Comparison of dynamic updates, the weighted version, and constant . Left: number of iterations. Right: load-displacement curves.
Figure 4: Example 1. Comparison of different for the dynamic scheme.

Asymmetrically notched three point bending test:  The configuration is shown in Figure 1 (right). The initial mesh is and times uniformly refined, yielding and mesh elements with the minimal mesh size parameter and . As material parameters, we use = , = , and = . Furthermore, we set [mm] and .

Figure 5 presents the number of iterations and the load-displacement curves. The number of iterations is decreasing from (in the figures cut to ) for the classical L-scheme, to a maximum of when using the dynamic updates. The choice of weighting does not seem to have a significant influence on the number of iterations though. The crack starts growing a bit later when using the dynamic updates, which can be inferred from the right plot in Figure 5. Thus, the stabilization parameters have a slight influence on the physical solution. This can be explained in the following way. In regions where the solution component is not uniquely defined. This leads to a sub-optimal convergence behaviour of the L-scheme. With the dynamic L-scheme we regain uniqueness, but at the cost of a slightly modified physical problem.

Figure 5: Example 2: Left: The number of iterations for the different schemes; the results for and are taken from BrWiBeNoRa19 . Right: The load-displacement curves; a slight difference can be observed in the results, indicating that the dynamic updates lead to a slight delay in the prediction of the starting time for the fracture growth.
Figure 6: Examples 1 and 2 for the dynamic scheme using ; three different mesh levels are used in order to verify the robustness of the proposed scheme. The results indicate that the mesh size does not influence the number of the iterations.
Remark 2

Noteworthy, the number of iterations for the dynamic L-scheme is robust with respect to the mesh refinement, as shown in Figure 6. This is in line with the analysis in BrWiBeNoRa19 ; list2016study ; MR2079503 , where it is proved that the convergence rate does not depend on the spatial discretization.


TW is supported by the German Research Foundation, Priority Program 1748 (DFG SPP 1748) under the project No. 392587580. CE is supported by the German Research Foundation, via Priority Program 1648 (DFG SPP 1648) under the grant No. EN-1042/2-2 and via EXC 2044-390685587, Mathematics Münster: Dynamics-Geometry-Structure. ISP is supported by the Research Foundation-Flanders (FWO), Belgium through the Odysseus programme (project G0G1316N).


  • [1] E. L. Allgower and K. Georg. Numerical continuation methods: an introduction. Springer, 1990.
  • [2] D. Arndt, W. Bangerth, T. C. Clevenger, D. Davydov, M. Fehling, D. Garcia-Sanchez, G. Harper, T. Heister, L. Heltai, M. Kronbichler, R. M. Kynch, M. Maier, J.-P. Pelteret, B. Turcksin, and D. Wells. The deal.II library, version 9.1. J. Numer. Math., 2019. accepted.
  • [3] M. K. Brun, T. Wick, I. Berre, J. M. Nordbotten, and F. A. Radu. An iterative staggered scheme for phase field brittle fracture propagation with stabilizing parameters. arXiv preprint arXiv:1903.08717, accepted for publication in Comp. Meth. Appl. Mech. Engrg., 2019.
  • [4] T. Gerasimov and L. D. Lorenzis. A line search assisted monolithic approach for phase-field computing of brittle fracture. Comp. Meth. Appl. Mech. Engrg., 312:276 – 303, 2016.
  • [5] F. List and F. A. Radu. A study on iterative methods for solving Richards’ equation. Comput. Geosci., 20(2):341–353, 2016.
  • [6] A. Mesgarnejad, B. Bourdin, and M. Khonsari. Validation simulations for the variational approach to fracture. Comp. Meth. Appl. Mech. Engrg., 290:420 – 437, 2015.
  • [7] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations. Int. J. Numer. Methods Engrg., 83:1273–1311, 2010.
  • [8] I. S. Pop, F. Radu, and P. Knabner. Mixed finite elements for the Richards’ equation: linearization procedure. J. Comput. Appl. Math., 168(1-2):365–373, 2004.
  • [9] M. Wheeler, T. Wick, and W. Wollner. An augmented-Lagangrian method for the phase-field approach for pressurized fractures. Comp. Meth. Appl. Mech. Engrg., 271:69–85, 2014.
  • [10] T. Wick. An error-oriented Newton/inexact augmented Lagrangian approach for fully monolithic phase-field fracture propagation. SIAM J. Sci. Comput., 39(4):B589–B617, 2017.