Log In Sign Up

Truncated Nonsmooth Newton Multigrid for phase-field brittle-fracture problems

We propose the Truncated Nonsmooth Newton Multigrid Method (TNNMG) as a solver for the spatial problems of the small-strain brittle-fracture phase-field equations. TNNMG is a nonsmooth multigrid method that can solve biconvex, block-separably nonsmooth minimization problems in roughly the time of solving one linear system of equations. It exploits the variational structure inherent in the problem, and handles the pointwise irreversibility constraint on the damage variable directly, without penalization or the introduction of a local history field. Memory consumption is significantly lower compared to approaches based on direct solvers. In the paper we introduce the method and show how it can be applied to several established models of phase-field brittle fracture. We then prove convergence of the solver to a solution of the nonsmooth Euler-Lagrange equations of the spatial problem for any load and initial iterate.


page 1

page 2

page 3

page 4


A Damped Newton Algorithm for Generated Jacobian Equations

Generated Jacobian Equations have been introduced by Trudinger [Disc. co...

A "parallel universe" scheme for crack nucleation in the phase field approach to fracture

Crack nucleation is crucial in many industrial applications. The phase f...

Preconditioners for Saddle Point Problems on Truncated Domains in Phase Separation Modelling

The discretization of Cahn-Hilliard equation with obstacle potential lea...

Interior-point methods for the phase-field approach to brittle and ductile fracture

The governing equations of the variational approach to brittle and ducti...

An accelerated staggered scheme for phase-field modeling of brittle fracture

There is currently an increasing interest in developing efficient solver...

A Particle-in-cell Method for Plasmas with a Generalized Momentum Formulation

In this paper we formulate a new particle-in-cell method for the Vlasov-...

1. Introduction

The equations of phase-field models of brittle fracture present a number of challenges to the designers of numerical solution algorithms [1]

. Even in the small-strain case the equations are nonlinear, due to the multiplicative coupling of the mechanical stresses to the degradation function. At the same time the non-healing property introduces an inequality constraint. Finally, eigenvalue-based splittings of the energy density as in 

[21] make the equations nondifferentiable.

In this paper we focus on the spatial problems of small-strain brittle-fracture phase-field models obtained by a suitable time-discretization. The standard approach to solving these spatial problems is based on operator splitting. Algorithms based on this approach, also known as staggered schemes, alternate between solving a displacement problem with fixed damage and a damage problem with fixed displacement. Both subproblems are elliptic and well-understood, and such methods are therefore straightforward to implement. The method can be interpreted as a nonlinear Gauß–Seidel method [8], which provides a natural framework for convergence proofs. Applications of the operator splitting scheme and its extensions appear, e.g., in [4, 6, 20].

In contrast, other works propose monolithic solution schemes based on Newton’s method [9, 8, 32, 31]. For the unmodified Newton method only local convergence can be shown, and failure to converge for large load steps is readily observed in practice [34]. Therefore, various authors have proposed extensions or modifications of the Newton idea to stabilize the method. In [9] a line search strategy is applied to enlarge the domain of convergence of Newton’s method. Wu et al. [34] propose to use the BFGS quasi-Newton algorithm, claiming that it is more stable than Newton’s method and more efficient than operator splitting. Heister et al. [16] proposed a modified Newton scheme which was later improved by Wick [32] with an adaptive transition from Newton’s method to the modified Newton scheme. Finally, in [27, 19], the authors suggest an arc-length method based on the fracture surface, and an adaptive time stepping scheme to enhance the robustness.

In summary, while monolithic Newton-type methods are reported to be faster than operator-splitting algorithms, the latter ones are more robust [1, 34].

Various approaches are used by these methods to deal with the damage irreversibility. A natural approach is to regularize the constraint, as investigated in [21, 30]. This leads to an additional parameter, and to ill-conditioned tangent matrices [10]. A second approach considers the thermodynamic driving force of the fracture phase-field as a global unknown yielding a three-field formulation which results in a saddle-point principle [21]. In a third formulation one considers the Karush–Kuhn–Tucker conditions and shifts the thermodynamic driving force of the fracture phase-field into a local history field representing the maximum over time of the elastic energy [20]. The nondifferentiable maximum function is then discretized explicitly in an otherwise implicit time discretization (Section Appendix: Reformulation with a history field). Unfortunately, this approach spoils the variational structure of the spatial problems.

Augmented Lagrangian solvers as in [30] introduce extra variables. Closest in spirit to the present manuscript is the use of bound-constrained optimization solvers, used, e.g., in [2, 8, 33]. None of these approaches are fully satisfactory.

The effect of the nondifferentiable terms caused by anisotropic splits of the mechanical energy density as in [21] is rarely discussed in the literature. Hybrid formulations like the one proposed in [1] try to overcome the additional computational cost of these anisotropic models, again at the cost of sacrificing the variation structure.

All these approaches are slow in the sense that they have to solve global partial differential equations at each Newton or operator-splitting iteration. When these methods use direct sparse solvers for the linear tangent problems, memory consumption can become problematic, too. At the same time, the problem of small-strain phase-field brittle-fracture has a lot of elegant variational structure; in particular, it fits directly into the rate-independent framework of

Mielke and Roubíček [22]. As a consequence, implicit time discretization leads to a sequence of coercive minimization problems for the deformation and damage fields together. These problems are not convex, but they are biconvex, i.e., convex (even strongly convex) in each variable separately. Pointwise inequality restrictions to handle the irreversibility of the fracture process as proposed in [21] reduce the smoothness of the objective functional, but do not otherwise influence its convexity or coercivity properties. The same holds for anisotropic energy splits based on linear quantities or the eigenvalues of the mechanical strain.

In recent years, nonsmooth multigrid methods have shown to be able to solve nonsmooth problems from mechanics efficiently without the need for solving global problems [14, 26, 15]. This makes them vastly more efficient than operator-splitting or Newton-based methods. As there are no sparse matrix factorizations, memory consumption remains low. In addition, these methods can be shown to converge globally, i.e., from any initial iterate and for any load step. This works by exploiting the above-mentioned variational structure, together with certain separability properties. As one such method, TNNMG can treat the pointwise constraints of the increment problems directly, i.e., without artificial penalization or simplifications by specific time discretizations. The idea is that TNNMG only needs to handle these constraints in a series of low-dimensional subproblems, each of which is easy to solve by itself. As a consequence, solving the problems with constraints is not appreciably slower than solving the corresponding unconstrained problem.

In this paper we show how the TNNMG method can be used to solve small-strain brittle fracture problems. This involves in particular verifying that the increment functionals have the required convexity and smoothness properties. We show this for a range of different degradation functions and local crack surface densities (including both, the standard Ambrosio–Tortorelli-2 functional and the more involved Ambrosio–Tortorelli-1 functional). We also show the results for isotropic elastic energies and for energies with a spectral split, using the theory of spectral functions.

This paper is organized as follows: Chapter 2 discusses a framework of small-strain phase-field models for brittle fracture, and shows the range of applicability of the TNNMG solver. Chapter 3 introduces the natural fully implicit time discretization, and proves existence of solutions for the spatial problems. In both chapters we pay particular attention to the mathematical properties of the energy functionals. In Chapter 4, finally, we introduce the TNNMG method. We explain its construction, discuss various algorithmic options, and prove that it converges globally to stationary points of the increment energy functional.

2. Phase-field models of brittle fracture

This chapter presents a range of phase-field models for brittle fracture, and discusses its smoothness and convexity properties.

Consider a deformable -dimensional object represented by a domain . The deformation of such an object is characterized by a displacement field

. The object is supposed to exhibit small-strain deformations and elastic material behavior only, and we therefore introduce the linearized strain tensor

. Following [21], we model the fracturing by a scalar damage field , where signifies intact material, and a fully broken one. Dirichlet boundary conditions can be posed both for the displacement and for the damage field. For this we select two not necessarily equal subsets of the domain boundary, and require

where and are two given functions.

Displacement and damage field evolve together, governed by a system of coupled nonsmooth partial differential equations. Disregarding inertia effects, we obtain a rate-independent system in the sense of Mielke and Roubíček [22]. Such a system can be written using the Biot equation


where means the Gâteaux derivative with respect to the second and third arguments of , and is the convex subdifferential with respect to the second argument of the dissipation potential .

In this equation, is a potential energy, which we assume to be of the form


The term is a degraded elastic energy density, and will be discussed in detail in Section 2.1. The term models the local crack surface density, and will be discussed in Section 2.2. The number is Griffith’s critical energy release rate, a material parameter. represents time-dependent volume and surface forces, which drive the evolution. We assume that is linear and -continuous in , and differentiable in with bounded time derivative.

The last term of (2) implements the restriction that the damage field can only assume values between and . For a set we define the indicator functional

For a closed, convex, nonempty set , the functional is convex, lower semicontinuous, and proper. Adding the constraint explicitly is not always necessary, as some fracture models lead to evolutions that satisfy the constraints implicitly. However, as pointwise bounds come with practically no cost when using the TNNMG solver, we do include them to extend our range of models.

To make the potential energy well defined, we will in general consider it on the first-order Sobolev space . Incorporating the boundary conditions leads to the affine subspace

The second term of the Biot equation (1) is , where is the dissipation potential


It implements the pointwise non-healing condition as proposed by [21]. Note that is convex and lower semicontinuous, and . The fact that is positively 1-homogeneous implies the rate-independence of the system.

Remark 2.1.

In the engineering literature, the same problem is frequently formulated as

with the rate potential


This formulation is equivalent to the Biot equation (1) if the problem is sufficiently smooth.

We now in turn discuss the potential energy and the dissipation potential.

2.1. Degraded elastic energy density

We consider models that behave linearly elastic and isotropic if the material is in an undamaged state. That is, for the undamaged stored energy density we use the St. Venant–Kirchhoff material law, whose energy density is given by

with Lamé parameters and .

The undamaged energy density is split additively into a part that produces damage and another part that does not. The damage-producing part is then scaled by a so-called degradation function

and the energy density takes the form


where is the space of symmetric matrices. The residual stiffness guarantees a well-posed problem in case of fracture.

Various different degradation functions have appeared in the literature [28, 17]. While the details vary, there appears to be agreement on the following properties:

Assumption 2.1.

The degradation function is differentiable, monotone decreasing, and fulfills and .

Note that several authors require in order to ensure that the evolution does not lead to values of larger than 1. We do not need this assumption here, because the pointwise constraint is enforced explicitly by the energy term (2).

The following specific degradation functions all fulfill Assumption 2.1:

(from [5])
(from [17])
(from [17])

Note that the functions and are strictly convex, but and are not even convex. For the rest of the paper we will restrict our considerations to convex twice continuously differentiable degradation functions .

Various splittings of have been proposed in the literature. We cover four common strain-based splittings taking the form (5).111The stress-based splitting of Steinke and Kaliske [28] is left for future work. All those splitting have the property that , and we will show that all have the following essential properties:

  1. for all and for all , i.e., is differentiable with locally Lipschitz continuous derivative.

  2. The gradient is semismooth for all .

  3. The gradient is globally Lipschitz continuous uniformly in , i.e., there exists independent of such that for all matrices we have

  4. is strongly convex uniformly in , i.e., there exists independent of such that for all matrices we have

  5. is coercive uniformly in in the sense that there exists independent of such that .

We remind that the gradient is called semismooth if the limit

exists for any point and any direction . The set denotes Clarke’s generalized Jacobian of the locally Lipschitz continuous map at (cf. [25]). Notice that the strong convexity d implies strong monotonicity of , i.e.,

For the splittings considered in the following we will only prove a and b directly, and show that the simplified assumptions of the following lemma hold true. This then implies c, d, and e.

Lemma 1.

Let and be convex, non-negative, and Lipschitz continuous. Then satisfies c, d, and e.


Let and be the Lipschitz constants of and , respectively. Then is Lipschitz continuous with uniform Lipschitz constant , because .

To show strong convexity, we first note that is strongly convex on with a modulus independent of . Now consider the function

for . Since this is a weighted sum of two convex functions and with non-negative weights and , it is itself convex. Thus, as a sum of this convex function and the strongly convex functions , the function is itself strongly convex and inherits the convexity modulus of . Finally, we note that with the same and we have

Despite those strong properties of we note that is not even convex in and together for any of the splittings considered below.

2.1.1. Isotropic splitting

In this model, any strain will lead to damage. The splitting is therefore


Without proof, we note the following simple properties of the energy density defined by (5) and this splitting:

Lemma 2.

The energy density defined in (5) with the isotropic splitting (6) has the properties ae. Furthermore has the stronger property that it is in and quadratic for all .

2.1.2. Volumetric decompositions

The isotropic model is unphysical, because it produces fracturing for all kinds of strain. In [18], Lancioni and Royer-Carfagni obtained better results by letting only the deviatoric strain contribute to the degradation. They introduced the split

with the deviatoric–volumetric strain splitting

With these definition, the energies are

Lemma 3.

The energy density defined in (5) with the isotropic volumetric splitting (7) has the properties ae. Furthermore has the stronger property that it is in and quadratic for all .


-smoothness and thus a and b are straightforward. The fact that and are quadratic, convex, and non-negative allows to derive c, d, and e from Lemma 1 and implies that is also quadratic. ∎

The decomposition of Lancioni and Royer-Carfagni is still isotropic. Amor et al. [2] proposed to only degrade the expansive part of the volumetric strain. Using the ramp functions

that provide the decompositions and , they proposed the energy split


where only the tensile volumetric strain contributes to damage.

Lemma 4.

The energy density defined in (5) with the anisotropic volumetric splitting (8) has the properties ae. Furthermore is not , unless .


We first note that the squared ramp functions are convex, with derivatives having a global Lipschitz constant , and piecewise (in the sense of [29, Definition 2.19]). Hence the functions are also with globally Lipschitz gradients and piecewise , which shows a and (using Lemma 1) c. Being piecewise implies semismoothness b of [29, Proposition 2.26]. Noting that , convexity of the squared ramp functions furthermore implies that the functions are also convex and non-negative, which by Lemma 1 provides d and e.

For the functional is quadratic and thus . In the case , if would be , then the function would also be . However, this function takes the form

and is thus piecewise quadratic but not in . ∎

2.1.3. Spectral decomposition

A more elaborate nonlinear splitting separating the tensile and compressive parts of the elastic energy was introduced in [21]. To define this splitting it is convenient to introduce the ordered eigenvalue function on the space of symmetric matrices, mapping any symmetric matrix

to the vector

containing its eigenvalues in ascending order. Using the ramp functions the tensile and compressive energies and are then defined as


Note that this indeed defines a splitting . For this splitting we will make the additional assumption that .

To quantify the properties of with respect to the strain tensor we use the theory of spectral functions. To this end we note that we can write as



The functions are symmetric in the sense that does not depend on the order of the entries of . Having this form we can infer properties of the functions from properties of the symmetric functions .

Lemma 5.

Let . Then the energy density defined in (5) with the spectral splitting (9) has the properties ae. Furthermore is not , unless .


We will first show ae. An essential ingredient is that the squared ramp functions are non-negative, piecewise quadratic, and convex.

a The squared ramp functions and thus are . Now [24, Proposition 4.3] shows that the spectral functions are also . Hence the same applies to .

b The squared ramp functions are piecewise functions. Hence the gradients are piecewise functions (in the sense of [29, Definition 2.19]) and thus semismooth [29, Proposition 2.26]. Now [24, Proposition 4.5] provides semismoothness of and thus of .

c Since the functions are piecewise quadratic and the gradients are globally Lipschitz continuous. Now Corollary 43 of [3] provides global Lipschitz continuity of the gradients of the spectral functions in the more general context of Euclidean Jordan algebras (which includes the special case of symmetric matrices). In fact, the Lipschitz constant of equals the one for if is equipped with the Frobenius norm. Using Lemma 1 this implies uniform Lipschitz continuity of .

d,e Since the functions are weighted sums of convex, non-negative squared ramp functions with nonnegative weights, they are convex and non-negative themselves. Convexity of the functions then follows from [3, Theorem 41] while non-negativity of those functions is trivial. Now Lemma 1 provides d and e.

To characterize second order differentiability of we first consider . Then coincides with the quadratic function and is thus . In the case , if would be , then the function for the fixed matrix with would also be . However, this function takes the form

and is thus piecewise quadratic but not in . ∎

Remark 2.2.

One can show that decomposes into finitely many disjoint subsets such that is twice continuously differentiable in the interior in each of these sets. A matrix is in the intersection of several if it either has an eigenvalue or if . While is not differentiable at those points, there are still generalized second-order derivatives. For example, the generalized Jacobian in the sense of Clarke contains the derivatives of with respect to all the adjacent sets . Semismoothness essentially means that such generalized derivatives provide an approximation that can be exploited in a generalized Newton method.

Remark 2.3.

The additional assumption is essential for convexity of . To see this we consider for the line segment

Then, along this line segment, is quadratic and takes the form

which is strictly concave for and sufficiently small .

2.2. Crack surface density

The crack surface density function per unit volume of the solid is typically of the form [23]

with parameters , , and , and a parameter function . The internal length scale parameter controls the size of the diffusive zone between a completely intact and a completely damaged material. For the regularized crack surface yields a sharp crack topology in the sense of -convergence. For a given function , the normalization constants and must be chosen such that the integral of over the fractured domain converges to the surface measure of the crack set as .

The function models the local fracture energy. Two types of local crack density functions appear in the literature. Double well potentials (as briefly reviewed in [1]), provide an energy barrier between broken and unbroken state, but will be disregarded here. Instead, we focus on the two widely used functionals




They are referred to in the literature as Ambrosio–Tortorelli (AT) functionals of type 1 and 2, respectively.

Some authors like [17] prefer because it has a local minimizer at . Thus, in the absence of mechanical strain, the unfractured solution is a minimizer of the total energy. As a result, no additional constraints need to be applied to ensure that . However, this argument becomes void when solver technology is available that can handle the explicit constraints . In contrast, for the AT-1 functional we have in the intact state . Together with the constraint this leads to a threshold, i.e., a minimum load required to cause damage [23].

Kuhn et al. [17] proposed to regard the Ambrosio–Tortorelli functionals as special instances of the general family defined by


with . The Ambrosio–Tortorelli functionals are obtained by setting for AT-1 and for AT-2. Further choices of are proposed in [23], which also do a detailed stability analysis for one-dimensional problems.

We note the following properties of the functional in (13):

Lemma 6.

The function given in (13) has the following properties:

  1. It fulfills and .

  2. It is strictly monotone increasing on for all .

  3. It is convex for all , and strictly convex for all .

For the rest of the paper we will assume the takes the form (13) with such that is guaranteed to be convex and quadratic.

3. Discretization and the algebraic increment potential

We use a fully implicit discretization in time, and Lagrange finite elements for discretization in space. Note that time-discretizations using a local energy history field [20] are typically only semi-implicit (cf. Section Appendix: Reformulation with a history field). By using a fully implicit time-discretization we retain the variational structure of the problem. Most of this chapter is spent investigating the properties of the increment functional.

3.1. Time discretization

It is shown in [22] that there is a natural time discretization for (1) that consists of sequences of minimization problems. Let the time interval be subdivided by time points , . We obtain a time-discrete formulation by integrating the rate potential given in (4) along paths on . Given initial values at , the continuous solution on this interval minimizes in the set of all paths of sufficient smoothness with  [22]. We then set to be . Defining as the set of suitably smooth paths in from to , and using that , this can be written as

The last integral term is the length of with respect to the Finsler norm given by . Minimizing over all paths gives the distance

Then one step of the time discretization scheme is given by

In our particular case (3), does not explicitly depend on . The dissipation distance is then easily computed as [22, Example 3.2.5]

We obtain the minimization problem


with the increment potential


Note that the time step size does not appear in this functional, because the model is rate-independent. Note also that the increment potential depends on the previous time step only through the indicator functional.

Lemma 7.

Assume that is non-trivial in the sense that its -dimensional Hausdorff-measure is positive. Then the functional is coercive on .


Using the uniform coercivity e of , , and for we get

for some constant . Using Korn’s inequality for , the Poincaré inequality for , and the fact that grows at most linearly we get for another constant

where we have used that the constraint implies in the second inequality. ∎

Lemma 8.

The functional is weakly lower semicontinuous on .


Since weak lower semicontinuity of the other terms in follows from convexity and lower semicontinuity of the integrands, we only need to consider the non-convex term


To this end we note that (16) can be written as for

and the density given by

Since is a Carathéodory function, non-negative (and thus uniformly bounded from below), and convex in  for all , it satisfies the assumptions of Theorem 3.4 in [7].

Now let be a weakly convergent sequence in . Then, by the compactness of the embedding into we get

Furthermore, the -weak convergence of implies -weak convergence of

because is in for each . Now Theorem 3.4 of [7] provides

As a direct consequence of coercivity and weal lower semicontinuity we get existence of a minimizer of the increment functional:

Theorem 3.1.

There is a solution to the minimization problem (14), i.e., there exists a global minimizer of .

3.2. Finite element discretization

The increment problem (14) of the previous section is posed on the pair of spaces for the displacements and for the damage variable. Let be a conforming finite element grid for . We discretize the function spaces by standard first-order Lagrangian finite elements. In order to derive an algebraic form of the discretized increment functional we make use of the standard scalar nodal basis associated to the grid nodes . Identifying the -valued and scalar finite element functions and with their coefficient vectors and , respectively, we write

where and . For the integration we use two kinds of quadrature rules: Integrals of smooth nonlinear terms over a grid element are approximated using a higher-order quadrature rule , while the integral over the nonsmooth term is approximated using the grid nodes as quadrature point, which is ofter referred to as lumping. Using these approximations we obtain the algebraic increment functional given by