Irreversibility and alternate minimization in phase field fracture: a viscosity approach

This work is devoted to the analysis of convergence of an alternate (staggered) minimization algorithm in the framework of phase field models of fracture. The energy of the system is characterized by a nonlinear splitting of tensile and compressive strains, featuring non-interpenetration of the fracture lips. The alternating scheme is coupled with an $L^{2}$-penalization in the phase field variable, driven by a viscous parameter~$\delta>0$, and with an irreversibility constraint, forcing the monotonicity of the phase field only w.r.t.~time, but not along the whole iterative minimization. We show first the convergence of such a scheme to a viscous evolution for $\delta>0$ and then consider the vanishing viscosity limit $\delta\to 0$.


Introduction
In the seminal work [13] the quasi-static propagation of brittle fractures in linearly elastic bodies is approximated in terms of equilibrium states of the Ambrosio-Tortorelli functional (1.1) Gε(u, z) := 1 2 Ω (z 2 + ηε)σ(u) : ǫ(u) dx + Gc where Ω is an open bounded subset of R n with Lipschitz boundary ∂Ω, u ∈ H 1 (Ω; R n ) is the displacement, ǫ(u) denotes the symmetric part of the gradient of u, σ(u) := Cǫ(u) is the stress, C being the usual elasticity tensor, ε and ηε are two small positive parameters, and Gc is the toughness, a positive constant related to the physical properties of the material under consideration (from now on we impose Gc = 1). The so-called phase field function z ∈ H 1 (Ω) is supposed to take values in [0, 1], where z(x) = 1 if the material is completely sound at x, while z(x) = 0 means that the elastic body Ω has developed a crack at x. Hence, the zero level set of z represents the fracture and z can be interpreted as a regularization of the crack set. In the static framework, the connection between (1.1) and fracture mechanics has been drawn in [6,16,18,25], where the authors showed the Γ-convergence of Gε as ε → 0 to the functional (1.2) G(u) := 1 2 Ω σ(u) : ǫ(u) dx + H n−1 (Ju) for u ∈ GSBD 2 (Ω; R n ) , where H n−1 denotes the (n − 1)-dimensional Hausdorff measure and Ju is the discontinuity set of u.
From the computational point of view, the study of the functional (1.1) is very convenient in combination with the so-called alternate minimization algorithm [9,12,13,14,15]: equilibrium configurations of the energy are computed iteratively by minimizing Gε first w.r.t. u and then w.r.t. z. A fracture irreversibility is imposed by forcing z to be non-increasing in time. Exploiting the separate convexity of Gε, the above scheme guarantees the convergence to a critical point of Gε, whose direct computation would be rather time consuming because of the non-convexity of the functional.
Let us turn our attention to the problem of evolution of the phase field driven by the energy functional Gε. In the context of rate-independent processes [33,36], the literature related to the general existence of such evolutions is very rich. The papers [21,39,40] showed the existence of energetic solutions, in which the equilibrium configurations are intended to be global minimizers of Gε. Global minimality is however expected to be unphysical in this context, since it generates time discontinuities in which the solution jumps between two equilibrium states, ignoring the presence of energetic barriers between them. For this reason, the works [10,29,30,31,37,38] promoted a vanishing viscosity approach, based on a viscous regularization of the evolution problem. As the viscosity tends to zero, the solution of the regularized problem converges to a Balanced Viscosity (BV) evolution [34,35], where the equilibrium states are critical points of Gε. In all the mentioned papers, the proof of existence is very constructive and is based on a scheme that only partly resembles the alternate minimization algorithm used in computational fracture mechanics: besides a time-discretization procedure, the authors made use of a one-step scheme, that is, at each time step the energy functional Fε (or a suitable viscous perturbation) is minimized in the pair (u, z), and a time-continuous solution is obtained in the limit as the time step tends to zero.
The theoretical investigation of the relationship between alternate minimization schemes for phase field models of fracture and rate-independent processes has instead only recently started with the works [1,2,3,28] (see also [27,32] for the application of alternating algorithms in different physical frameworks). Here, we describe the algorithm adopted in [28]. In dimension n = 2, given a time horizon T > 0, a time dependent Dirichlet boundary condition t → g(t), and a suitable initial condition (u0, z0), the authors considered, as usual in the study of rate-independent systems, a time discretization procedure: for every k ∈ N set τ k := T /k the time step increment and t k i := iτ k , i = 0, . . . , k, the time nodes. A time discrete evolution is then constructed through an iterative scheme: Known the state (u The analysis performed in [28] ensures that the above discrete solutions converge to a phase field evolution as the time increment τ k tends to zero. In particular, the limit solutions are characterized in terms of parametrized BV evolutions. We mention that a discrete version of [28] in a space-discrete (finite element) setting has been studied in [1] together with the limit of the solutions in a discrete to continuum sense, i.e., as the mesh becomes finer and finer. The results of [28] have been further generalized in [3] to not separately quadratic energy functionals. An example of such energies, that will be considered also in this paper, is inspired by the phase field model introduced in [7,17]. The underlying idea is that an elastic material behaves differently when under tension or compression. Moreover, fracture propagation is allowed only as a result of tensile or shear stresses, while compression does not lead to inelastic behaviors. Thus, in contrast with (1.1), the factor z 2 + ηε shall not affect the whole stress σ(u). Instead, a splitting of the strain ǫ(u) into its volumetric ǫv(u) := 1 2 (trǫ(u))I and deviatoric ǫ d (u) := ǫ(u) − ǫv(u) parts is considered, where tr denotes the trace of a matrix and I is the identity matrix. Introducing also tensile and compressive strains ǫ ± v (u) := 1 2 (trǫ(u))±I, the elastic energy density writes as In the above formulas, µ and κ are two positive parameters related to the Lamé coefficients of the elastic material, and hε, fε : R → [0, +∞) are two degradation functions. We refer to Section 2 for the whole set of hypotheses and to [7] for more details about the model. Here we only mention that in [17] it has been proven that, in dimension n = 2, Fε Γ-converges to G in (1.2) for u ∈ SBD 2 (Ω; R 2 ) satisfying the non-interpenetration constraint [u]·νu ≥ 0, where νu is the approximate unit normal to Ju and [u] stands for the amplitude of the jump of u across Ju. Despite the sound mathematical results obtained in [1,3,28], one of the drawback of the alternating scheme (1.3)-(1.4) lies in the irreversibility condition z ≤ z k i,j−1 , which forces the phase field variable to be non-increasing along the whole algorithm. This requirement, indeed, could lead to an accumulation of numerical error, making the fracture simulation very inaccurate. In order to bypass such a problem, the weaker irreversibility constraint z ≤ z k i−1 is usually numerically imposed (see, for instance, [9,12,13,14,15]). From a theoretical viewpoint, very little is known about convergence of the scheme (1.3)-(1.4) with this new irreversibility condition, mainly because the lack of monotonicity of z along the algorithm prevents from deducing any time regularity of the discrete evolutions, such as BV in time, and makes the analysis of [3,28] out of reach. Up to our knowledge, the only existing result is [2], where the minimization (1.4) is replaced by In particular, we notice that the functional Gε is perturbed with an L 2 -penalization, which makes the time discrete evolution regular in time. Furthermore, the minimum problem in (1.7) is unconstrained, while the irreversibility is a posteriori imposed by truncating the minimizerz k i,j with z k i−1 . In the limit as τ k → 0 we have been able to show the convergence to an L 2 -gradient flow of Gε, while the presence of the pointwise minimization prevented us from studying the vanishing viscosity limit, and hence the convergence to a BV evolution.
The scope of this note is to give a first result of convergence of an alternate minimization scheme to a quasistatic evolution, in the presence of the (weaker) irreversibility constraint z ≤ z k i−1 . Precisely, with the notation introduced in (1.5)-(1.6), for every δ > 0 and k ∈ N we consider the iterative algorithm . Indeed, we have now explicitly imposed in the minimization the monotonicity constraint z ≤ z k i−1 , which only ensures an irreversibility w.r.t. time but not along the whole scheme. As in (1.7), we have regularized in time the evolution of the phase field z by adding to Fε an L 2 -penalization driven by a small viscous parameter δ.
Following the lines of [19,29], in Sections 3 and 4 we study the limit of (1.8)-(1.9) as k → ∞ and δ → 0, in the given order. Hence, in Theorem 4.1 we show for δ > 0 the convergence of the iterative scheme to a viscous evolution (u δ , z δ ) satisfying the following displacement equilibrium and energy balance: where |∂ − z Fε| denotes the unilateral slope of Fε and the dot indicates the derivative w.r.t. time. We refer to Definitions 2.2 and 2.5 for the full details.
Finally, in Theorem 4.2 we consider the vanishing viscosity limit δ → 0 and prove that the pair (u δ , z δ ) converges, in a time reparametrized setting, to a BV evolution, now represented by a triple (t, u, z), where t is a suitable Lipschitz parametrization of the time interval [0, T ]. Besides the equilibrium condition (1.10), the triple (t, u, z) also satisfies the energy balance where the index ′ denotes the derivative w.r.t. the new time variable s. As in [19,29], the time reparametrization is based on a uniform estimate of the length of the curve t → z δ (t). Namely, in Lemma 3.5 we prove that the arc length of the algorithm (1.8)-(1.9), measured in terms of the distance between two consecutive states of the iterative minimization, is bounded uniformly w.r.t. k and δ. Eventually, this allows us to perform an arc length reparametrization s → t δ (s) of time which makes (u δ , z δ ) 1-Lipschitz continuous, and thus compact, in the new time variable. We notice that the uniform bound of the arc length is a consequence of the combination of a general regularity and continuity result [24] (see also Lemma 2.7) for PDEs with non-constant coefficients and of Sobolev embeddings, valid only in dimension two.

Notation and setting of the problem
In order to explain the phase field model considered in this work, we have to introduce some notation. Let M 2 denote the space of squared matrices of order 2 and M 2 s be the subspace of symmetric matrices. For every F ∈ M 2 , we consider its splitting in volumetric and deviatoric part Note that Fv : F d = 0, where the symbol : denotes the scalar product between matrices. As a consequence, we have that where | · | denotes Frobenius norm. Moreover, we set F ± v := 1 2 (trF)±I , where (·)+ and (·)− stand for positive and negative part, respectively. It is clear that s , we can rewrite the usual linear elastic energy density as where λ, µ are the Lamé coefficients and κ := λ + µ We assume µ, κ > 0.
Following the lines of [7,17], we consider a phase field model that does not allow for fracture under compression, that is, when (trE)− = 0. To model such a behavior, the phase field variable z ∈ [0, 1] is assumed to affect only the energetic contribution of the tensile strain E + v and of the deviatoric strain E d . Hence, the elastic energy density reads Notice that, under these assumptions, h is non-decreasing in [0, +∞).
For z fixed, the function W (z, ·) is differentiable w.r.t. E and We collect here some useful properties of the energy density W . We refer to [3, Lemma 3.1] for more details.
is of class C 1,1 loc . Moreover, there exist two positive constants c, C such that for every z ∈ [0, 1] and every E1, E2 ∈ M 2 s the following holds: Let Ω be an open bounded subset of R 2 with Lipschitz boundary ∂Ω. For later use, we also fix ∂DΩ ⊆ ∂Ω regular in the sense of Gröger [23]. For every u ∈ H 1 (Ω; R 2 ) and z ∈ H 1 (Ω) ∩ L ∞ (Ω) we define the elastic energy where ǫ(u) = 1 2 (∇u + (∇u) T ) denotes the strain. We introduce the dissipation potential associated to the phase field variable z ∈ H 1 (Ω) ∩ L ∞ (Ω) given by Here, we assume the degradation function f : R → [0, +∞) to be of class C 1,1 loc , strongly convex, and such that However, many different degradation functions have been extensively considered in the fracture mechanics literature (see, e.g., [4,26,41]).
The total energy F : of the system is given by the sum of elastic energy (2.1) and dissipation potential (2.2) Notice that, in comparison with (1.6), we have fixed ε = 1 2 . An important role in the definition of evolution we consider in this work is played by the following notion of unilateral L 2 -slope (see also [2, Definition 1.1]).
where the convergence is intended in the L 2 -topology.
Remark 2.3. The minus sign appearing in the notation |∂ − z F| reminds that only negative variations are allowed and it should not be confused with a similar notation for the relaxed slope (see, e.g., [5, Section 2.3]). (2.5) The natural relationship between partial derivatives (2.5) and slope (2.4) is stated in the next lemma, whose proof can be found, for instance, in [ We are now in a position to give the precise definition of viscous and vanishing viscosity evolutions we consider in this paper.
We say that a pair (u δ , z δ ) : [0, T ] → H 1 (Ω; R 2 ) × H 1 (Ω) is a viscous evolution for the energy F with initial condition (u0, z0) and boundary condition g if the following properties are satisfied: (c) Displacement equilibrium: for every t ∈ [0, T ] we have u δ (t) = g(t) on ∂DΩ and (d) Energy balance: the map t → F(u δ (t), z δ (t)) is absolutely continuous and for a.e. t ∈ [0, T ] it holds Our first goal is to prove the convergence to a unilateral L 2 -gradient flow of the time discrete solutions obtained by an alternating minimization scheme (see (3.1)-(3.2) for the algorithm and Theorem 4.1 for the convergence result). Our second aim is to characterize the limit δ → 0, of the above evolution. The limit solutions are described in the following definition.
For S ∈ (0, +∞) we say that a triple (t, u, z) is a vanishing viscosity evolution for the energy F with initial condition (u0, z0) and boundary condition g if the following properties are satisfied: (e) Energy balance: the map s → F(u(s), z(s)) is absolutely continuous and for a.e. s ∈ [0, S] it holds Next lemma provides a regularity property needed in our setting. The complete proof can be found in [3, Proposition 3.6]. For a more general statement we refer to [24].

Alternate minimization algorithm and a priori estimates
Let us fix a time horizon T > 0 and a Dirichlet boundary datum g ∈ H 1 ([0, T ]; W 1,p (Ω; R 2 )), with p ∈ (2, +∞), to be applied on the Dirichlet part of the boundary ∂DΩ ⊆ ∂Ω. For fixed δ > 0, the construction of a gradient flow of F is done by time-discretization (see, e.g., [5]). Here we couple this standard procedure with an alternate (or staggered) minimization algorithm, which is by now typical in computational fracture mechanics [9,12,13,14,15].
Let us describe the minimization scheme. For every k ∈ N we define the time step τ k := T /k and the time nodes t k i := iτ k for i ∈ {0, . . . , k}. For i = 1, . . . , k assume that we know the configuration (u k i−1 , z k i−1 ) at time t k i−1 . The new state (u k i , z k i ) at time t k i is constructed as limit of an alternate minimization procedure: for j = 0 we set u k i,0 := u k i−1 and z k i,0 := z k i−1 , while for j ≥ 1 we define We notice that the solutions of (3.1) and (3.2) exist and are unique by strict convexity.
Remark 3.1. We stress that the irreversibility constraint, stating that the damage (or the crack) can only grow and healing is not allowed, is expressed in (3.2) by the inequality z ≤ z k i−1 , which only involves the configuration at time t k i−1 . This is in contrast with the recent theoretical literature developed in [1,3,8,27,28], where the irreversibility has been implemented in a stronger way by imposing z ≤ z k i,j−1 . The latter choice, even if mathematically correct, could lead to the accumulation of numerical error in the simulation of the fracture process. For this reason, in most of the applications in computational fracture mechanics (see, e.g., [9,12,13,14,15]) the constraint in (3.2) is considered.
Again in contrast with [1,3,8,28], in (3.2) we have perturbed the functional F with an L 2 -penalization of the distance from z k i−1 . This is necessary in order to carry out our analysis, and leads in the limit as k → ∞ to the construction of a gradient flow of the energy F according to Definition 2.5. A quasi-static evolution is recovered only in the limit as the viscosity parameter δ tends to 0. We refer to Section 4 for the full discussion. We further notice that a similar penalization has been used also in [2,29,30,31,38]. However, only [2] deals with an alternate minimization procedure, where the minimization in (3.2) is unconstrained and the irreversibility is imposed a posteriori by a simple pointwise minimization, which made the study of the viscosity limit δ → 0 not accessible.
In the following proposition we prove the convergence, up to subsequence, of the sequence (u k i,j , z k i,j ). Proof. By definition of u k i,j and of z k i,j we have that Hence, the sequences u k i,j and z k i,j are bounded in H 1 , uniformly w.r.t. j. Thus, there exist a subsequence j l ,ū ∈ H 1 (Ω; R 2 ), andz ∈ H 1 (Ω; [0, 1]) such that u k i,j l ⇀ū weakly in H 1 (Ω; R 2 ) and z k i,j l ⇀z weakly in H 1 (Ω) as l → ∞. Up to a further (not relabeled) subsequence, we may assume that the above convergences are strong in L 2 (Ω) and that z k i,j l −1 ⇀ẑ weakly in H 1 (Ω), for someẑ ∈ H 1 (Ω; [0, 1]). By Lemma 2.10 and by the above convergences, we immediately deduce that u k i,j l →ū in H 1 (Ω; R 2 ) and that u = arg min {F(u,ẑ) : u ∈ H 1 (Ω; R 2 ), u = g(t k i ) on ∂DΩ} . Moreover, we have that (3.4) is satisfied. Indeed, by (3.2) we have that for every z ∈ H 1 (Ω) with z ≤ z k 2 . Therefore, passing to the liminf as l → ∞ in the previous inequality, applying Lemma 2.8, and exploiting the strong convergence of u k i,j l we get (3.4). It remains to prove thatẑ =z. To do this, we notice that Hence, passing to the liminf in the previous chain of inequalities and applying again Lemma 2.8 we obtain By uniqueness of minimizer, this implies thatẑ =z, and the proof is thus concluded.
In view of Proposition 3.2, we are allowed to define where the limits are intended to be up to a subsequence and strong in H 1 (Ω; R 2 ) and weak in H 1 (Ω), respectively. In particular, u k i and z k i solve In the following two propositions we prove the boundedness of u k i and z k i together with a discrete energy inequality.
We now define the following interpolation functions: We notice that at this point, arguing as in [2], we could already show the convergence of the sequence (ū δ k , z δ k ) to a viscous evolution (u δ , z δ ), without passing through a priori estimates showing time regularity of z δ k . The most delicate point would indeed be the energy balance (d) of Definition 2.5, which would anyway follow from convexity and Riemann sum arguments as in [2]. However, such an analysis would preclude the study of the limit δ → 0.
Therefore, in the next two lemmas we provide some a priori estimates for the discrete evolutions defined in (3.9)-(3.11). Precisely, we show in Lemma 3.4 that z δ k is bounded in H 1 ([0, T ]; H 1 (Ω)) w.r.t. k. Since the bound is not uniform in δ > 0, we further prove in Lemma 3.5 that the length of the curve t → z δ k (t), namely is bounded uniformly w.r.t. k and δ. This will allow us to consider the limit δ → 0 in Theorem 4.2.
Lemma 3.4. There exists a positive constant C independent of δ and k such that for every t ∈ [0, T ] Proof. We follow here the lines of [29,Propositions 4.6 and 5.7] and [19,Proposition 2.8]. For δ and k fixed, for simplicity of notation we setżi := In view of (3.6) we have that Let us consider i ≥ 2. Insertingżi as a test function in (3.14) at time t k i−1 we have that Subtracting (3.15) from the previous inequality we get (3.16) where, in the last step, we have used the inequality 2a We now estimate the left-hand side of (3.16). First, we split the difference using the definition (2.3) of F, so that

(3.17)
Since f is strongly convex, there exists c ∈ (0, 1) such that As for the first two terms on the right-hand side of (3.17), we write where, in the last inequality, we have used the convexity of h, which yieldsżi(h ′ (z k i−1 ) − h ′ (z k i )) ≤ 0. Let r, ℓ ∈ (2, +∞) be as in Lemma 2.7. Since 0 ≤ z k i ≤ 1, by Hölder inequality we continue in (3.19) with for 1 ν + 2 r = 1. By Lemma 2.7 we know that u k i W 1,r is bounded uniformly w.r.t. k, i, and δ. Furthermore, we deduce from (3.20) that for some positive constant C independent of k, i, and δ. Setting λ := max{ν, ℓ} and applying Young inequality to (3.21) we infer that żi 2 H 1 ≤ Cτ k ( żi 2 2 + ġi 2 W 1,p ) for some C > 0 independent of k and δ.
Lemma 3.5. There exists C > 0 independent of k and δ such that Proof. We use here the same notation as in the proof of Lemma 3.4. Arguing as in (3.16), for i ≥ 2 we get that where the last step is due to Cauchy and triangle inequality. In order to estimate the left-hand side of (3.27) we proceed as in (3.17)-(3.21), obtaining for some positive c, C independent of i, k, and δ and some λ ∈ (2, +∞). By interpolation and Sobolev inequality we deduce from (3.28) that for ε, Cε > 0. Rewriting (3.29) for ε < c/2 and using the inequality żi 2 ≤ żi H 1 we get Multiplying (3.30) by 2 and dividing by δ we finally obtain an inequality of the form for ai := żi 2, bi := cτ k 2δ 1 2 żi H 1 , ci := 2Cτ k δ 1 2 ġi W 1,p , di := Cτ k δ żi 1, and γ := cτ k 4δ . Applying now the discrete Gronwall inequality of [29,Lemma 5.9] and performing the computation contained in [19,Proposition 2.8] we eventually deduce that An estimate for ż1 H 1 follows directly from (3.24) and an application of Sobolev and interpolation inequalities, so that Hence, summing up (3.31) and (3.32) and applying Hölder inequality toġi we obtain (3.33) We conclude by observing that in view of (3.12) for t = t k 1 = τ k we have that δ ż1 2 ≤ Ce C δ τ k for some C > 0 independent of k and δ. Hence, for τ k ≤ δ we infer that δ ż1 2 is bounded and, in view of (3.33), which implies the thesis.
In the next proposition we provide a discrete energy inequality.
Proposition 3.6. Let δ > 0 be fixed. Then, the following facts hold: ) on ∂DΩ} ; (c) There exists a constant C > 0 such that for every k ∈ N and every t ∈ [0, T ] Proof. We notice that condition (b) is true by definition ofū δ k andz δ k and by (3.5). Let us prove the energy estimate (c). For k ∈ N, t ∈ (0, T ] fixed, let i ∈ {1, . . . , k} be such that t ∈ (t k i−1 , t k i ]. By convexity of z → F(u k i , z) and by applying Proposition 3.3 we get We now have to pass from u k i to u k i−1 in the left-hand side of (3.35). We perform this passage in two steps, first moving from u k i to u k i,1 and then from u k i,1 to u k i−1 . For the second step we can simply rely on the minimality (3.1) of u k i,1 . For the first step, instead, we make use of the regularity estimate of Lemma 2.7. Hence, by convexity of u → F(u, z) we have that (3.36) We estimate the right-hand side of (3.36). By minimality of u k i we have that Hence, choosing r > 2 as in Lemma 2.7, from the regularity of h and Sobolev embedding we deduce that where C is a positive constant independent of k and δ and 1 ν + 2 r = 1. In a similar way, by Lemma 2.1, we obtain that for every s In particular, R k → 0 as k → ∞ thanks to the bound (3.13) of Lemma 3.4 and to the regularity of g.

Viscous and quasistatic evolutions
In this section we show the existence of a viscous evolution in the sense of Definition 2.5 for every δ > 0. Furthermore, we study the limit as δ → 0 of the above evolutions proving that, in a suitable time-reparametrized setting, they converge to a vanishing viscosity evolution in the sense of Definition 2.6.
The energy equality (d) of Definition 2.5 follows by Young inequality. Finally, we notice that (4.4) and the energy balance imply (4.2), and the proof is thus concluded.
We can now conclude showing that in the limit as δ → 0 the viscous evolutions determined in Theorem 4.1 converge to a vanishing viscosity evolution (see Definition 2.6). In order to do this, we first introduce an arc length reparametrization of time which makes z δ Lipschitz continuous. For every t ∈ [0, T ] we define σ δ (t) := t + t 0 ż δ (τ ) H 1 dτ .