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 L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^{2}$$\end{document}-penalization in the phase field variable, driven by a viscous parameter δ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta >0$$\end{document}, 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 δ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta >0$$\end{document} and then consider the vanishing viscosity limit δ→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \rightarrow 0$$\end{document}.


Introduction
In the seminal work [16], the quasi-static propagation of brittle fractures in linearly elastic bodies is approximated in terms of equilibrium states of the Ambrosio-Tortorelli functional 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 G c is the toughness, a positive constant related to the physical properties of the material under consideration (from now on we impose G c = 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 [9,19,21,29], where the authors showed the Γ-convergence of G ε as ε → 0 to the functional where H n−1 denotes the (n − 1)-dimensional Hausdorff measure and J u 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 [12,[15][16][17][18]: 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 [39,42], the literature related to the general existence of such evolutions is very rich. The papers [25,46,47] 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 [13,[34][35][36]43,44] (see also [1,4,5,23,32,37,45,48] for further applications to fracture mechanics) 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 [40,41], 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 [2,3,6,33] (see also [31,38] for the application of alternating algorithms in different physical frameworks). Here, we describe the algorithm adopted in [33]. In dimension n = 2, given a time horizon T > 0, a timedependent Dirichlet boundary condition t → g(t), and a suitable initial condition (u 0 , z 0 ), 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 ( In the limit j → ∞, the algorithm (1.3)-(1.4) identifies a critical point (u k i , z k i ) of G ε . The analysis performed in [33] 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 [33] in a space-discrete (finite element) setting has been studied in [2] 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 [33] have been further generalized in [6] 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 [10,20]. 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 and the phase field energy becomes In the above formulas, μ and κ are two positive parameters related to the Lamé coefficients of the elastic material, and h ε , f ε : R → [0, +∞) are the degradation and the crack geometry function, respectively ZAMP Irreversibility and alternate minimization in phase field fracture Page 3 of 21 128 (see, e.g., [49]). We refer to Sect. 2 for the whole set of hypotheses and to [10] for more details about the model. Here we only mention that in [20] 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 J u and [u] stands for the amplitude of the jump of u across J u . Despite the sound mathematical results obtained in [2,6,33], 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, [12,[15][16][17][18]). 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 [6,33] out of reach. Up to our knowledge, the only existing result is [3], 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 quasi-static 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 [22,34], in Sects. 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.4 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 [22,34], the time reparametrization is based on a uniform estimate of the length of the curve t → z δ (t). Namely, in Lemma 3.6 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 [28] (see also Lemma 2.6) 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 F v : 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 s , we can rewrite the usual linear elastic energy density as where λ, μ are the Lamé coefficients and κ := λ + μ. We assume μ, κ > 0.
Following the lines of [10,20], 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 [6, Lemma 3.1] for more details.
Moreover, there exist two positive constants c, C such that for every z ∈ [0, 1] and every E 1 , E 2 ∈ 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 [27]. For every u ∈ H 1 (Ω; R 2 ) and z ∈ H 1 (Ω) ∩ L ∞ (Ω), we define the 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 crack geometry function f : R → [0, +∞) to be of class C 1,1 loc , strongly convex, and such that We notice that, by the assumptions on f , it holds However, many different functions have been extensively considered in the fracture mechanics literature (see, e.g., [7,30,49]). 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 [3, Definition 1.1]).
where the convergence is intended in the L 2 -topology.
For u ∈ H 1 (Ω; R 2 ) and z, ϕ ∈ H 1 (Ω) ∩ L ∞ (Ω) there exists finite the partial derivative of F with respect to z, i.e., 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 [  Finally, let us define, for u ∈ H 1 (Ω; R 2 ) and z ∈ H 1 (Ω) ∩ L ∞ (Ω), the functional 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 δ ): is a viscous evolution for the energy F with initial condition (u 0 , z 0 ) and boundary condition g if the following properties are satisfied: (a) Time regularity: 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. Next lemma provides a regularity property needed in our setting. The complete proof can be found in [6,Proposition 3.6]. For a more general statement, we refer to [28]. Then there exist an exponent 2 < r < p and a constant C > 0 such that for every t 1 , t 2 ∈ [0, T ] and every We collect here some technical results that will be useful in the next sections.
Proof. The lower semicontinuity of D is obvious by convexity. The semicontinuity of E follows for instance from [24,Theorem 7.5].
Now we state a semicontinuity property of the slope |∂ − z F|. The proof can be found in [6, Lemma 3.9].

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., [8]). Here we couple this standard procedure with an alternate (or staggered) minimization algorithm, which is by now typical in computational fracture mechanics [12,[15][16][17][18]. 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

2)
We notice that the solutions of (3.1) and (3.2) exist and are unique by strict convexity.
128 Page 8 of 21 S. Almi ZAMP 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 [2,6,11,31,33], 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., [12,[15][16][17][18]) the constraint in (3.2) is considered.
Again in contrast with [2,6,11,33], 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.4. A quasistatic evolution is recovered only in the limit as the viscosity parameter δ tends to 0. We refer to Sect. 4 for the full discussion. We further notice that a similar penalization has been used also in [3,[34][35][36]44]. However, only [3] 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 ).

Remark 3.3.
Notice that the function z in (3.4) satisfies z ≥ 0, as a consequence of the hypotheses on the functions h and f . Indeed, it holds for every u ∈ H 1 (Ω; R 2 ) and every z ∈ H 1 (Ω) with z ≤ z k i−1 .

Proof of Proposition 3.2.
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.9 and by the above convergences, we immediately deduce that u k i,j l →ū in H 1 (Ω; R 2 ) and thatū 2 . Therefore, passing to the liminf as l → ∞ in the previous inequality, applying Lemma 2.7, and exploiting the strong convergence of u k i,j l we get (3.4).

ZAMP
Irreversibility and alternate minimization in phase field fracture Page 9 of 21 128 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.7 we obtain 2 . 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 [3], 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.4, which would 128 Page 10 of 21 S. Almi ZAMP anyway follow from convexity and Riemann sum arguments as in [3]. 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.5 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.6 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.5. There exists a positive constant C independent of δ and k such that for every t ∈ [0, T ] 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 19) 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.6. Since 0 ≤ z k i ≤ 1, by Hölder inequality, we continue in (3.19) with for 1 ν + 2 r = 1. By Lemma 2.6, 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 Finally, by interpolation and Sobolev inequality we obtain for ε, C ε > 0. Choosing ε = c 2 in (3.22) and combining (3.16), (3.18), and (3.22) we deduce that for some C > 0 independent of k and δ.
Lemma 3.6. There exists C > 0 independent of k and δ such that Proof. We use here the same notation as in the proof of Lemma 3.5. 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 (3.31) An estimate for ż 1 H 1 follows directly from (3.24) and an application of Sobolev and interpolation inequalities, so that τ k ż 1 H 1 ≤ Cτ k ( ż 1 1 + ġ 1 W 1,p ) . 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. where R k ≥ 0 is such that R k → 0 as k → ∞.
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.4, we get (3.35) 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 128 Page 14 of 21 S. Almi ZAMP 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.6. 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.6, 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 ∈ [t k i−1 , t k i ] (3.38)