Convergence Analysis of a Local Stationarity Scheme for Rate-Independent Systems and Application to Damage

This paper is concerned with an approximation scheme for rate-independent systems governed by a non-smooth dissipation and a possibly non-convex energy functional. The scheme is based on the local minimization scheme introduced in [EM06], but relies on local stationarity of the underlying minimization problem. Under the assumption of Mosco-convergence for the dissipation functional, we show that accumulation points exist and are so-called parametrized solutions of the rate-independent system. In particular, this guarantees the existence of parametrized solutions for a rather general setting. Afterwards, we apply the scheme to a model for the evolution of damage.


Introduction
The effect of rate-independence occurs in various different areas of mechanics.This concerns for example the field of elastoplasticity, damage and shape-memory, to only mention a few (see, e.g., [5,10,16,18,19]).One main characteristic of such systems is the fact that changes in the state are solely driven by an external force.What is more, as the name already suggests, the system is independent of the rate at which the loading is applied, that is to say, whenever  is a solution to some external load ℓ, then  ∘  is a solution to ℓ ∘  for every monotone increasing time rescaling .
In this paper we consider rate-independent systems that can be described by the following differential inclusion 0 ∈ ℛ( ż()) +   ℐ(, ()) a.e. in [0,  ]. (RIS) One may see this inclusion as a balance of forces, i.e., the dissipative force ℛ and the potential force −  ℐ(, ) must annihilate each other.Implicitly hidden within this formulation is the fact that the potential force as well as the dissipative force result from a (possibly non-convex) energy functional ℐ and a dissipation functional ℛ.While we postpone the exact assumptions to Section 2, let us mention at this point that the characteristic feature of the formulation in (RIS) is the positive 1-homogeneity of the dissipation ℛ.It is this property which induces that (RIS) is indeed rate-independent.However, the combination of non-convex energies and positive 1-homogeneous dissipations allow the formation of abrupt changes in the state, even if the external forces evolve smoothly.Hence, suitable notions of solutions for (RIS) need to be able to handle temporal discontinuities.One such concept are the so-called parametrized BV-solutions, whose exact definition is given in Definition 2.4.
Loosely speaking, such solutions are considered as curves in the extended phase space [0,  ]× and parametrized by arc-length.The jump path from one state to the other thus becomes an integral part of the solution itself.This idea was first applied in [7,20,21] for systems with dry friction and later on generalized in [9,25,26] for finite and infinite dimensional problems, respectively.Particularly, in [9], the authors introduced the following time-incremental local minimization scheme for the approximation of parametrized BV-solutions: = min{ −1 +  − ‖  −  −1 ‖  ,  }. (1.1b) It was moreover shown that, for  ↘ 0, subsequences of discrete solutions generated by (1.1) (weakly) converge to a parametrized BV-solution.While the authors in [9] considered a finite dimensional setting, in [28,29] and particularly [15] the results have been generalized to the infinite dimensional problems at least for semilinear energies.Furthermore, in [22] the scheme was combined with a discretization in space.In this paper we extend these result to more general energy and dissipation functionals on the one side and, additionally, build the scheme upon stationary points rather than minimizers of (1.1).The actual scheme (LISS) is presented in Section 3. Let us underline that the consideration of stationary points instead of (global) minimizers is of major importance from a numerical point of view, since optimization algorithms can in general only compute stationary points.
What is more, we incorporate unbounded dissipations into our convergence analysis, which allows us to apply our scheme to a model for the evolution of (partial) damage.Now, let us shortly outline the paper.In Section 2, we introduce our notation and state the assumptions on the energy and the dissipation functional.Moreover, we recall the precise notion of parametrized BV-solutions.Section 3 is then devoted to the presentation of the actual local stationarity scheme and its convergence analysis.Particularly, since we allow the dissipation to be approximated by some functional ℛ  , we provide suitably adapted a priori estimates for the discrete solution which still meets a discrete version of an energy-equality.Building on that, we derive our main convergence result in Theorem 3.14.In Section 4, we then focus on a rate-independent damage model and describe the algorithmic realization of the discrete stationarity scheme based on a semismooth Newton-method.Finally, we present a numerical example and compare it with results from the literature.

Basic notations and standing assumptions
Let us start with some basic notations used throughout the paper.In the following,  > 0 always stands for a generic constant.Moreover, given two normed linear spaces ,  , we denote by ⟨•, •⟩  * , the dual pairing and suppress the subscript, if there is no risk for ambiguity.By ‖ • ‖  , we denote the norm in  and ℒ(,  ) is the space of linear and bounded operators from  to  .If  is even a Hilbert space, we write   :  * →  for the Riesz isomorphism.Furthermore,   (, ) is the open ball in  around  ∈  with radius  > 0. Given a convex functional  :  → R ∪ {∞}, we denote the (convex) subdifferential of  at  by  () ⊂  * and its conjugate functional by  * :  * → R ∪ {∞}.Finally, |Ω| stands for the Lebesgue measure of a set Ω ⊂ R  ,  ∈ N and R  ≥0 describes the set of vectors in R  whose components are greater or equal to 0.

Assumptions on the data
Let us now introduce the assumptions on the quantities in (RIS).We assume that the underlying spaces  ,  and  are Banach spaces with  ˓→ ,  ˓→  , where ˓→ , stands for a dense and compact embedding.Moreover,  and  are required to be reflexive and separable.
Note that the combination of (E1)-(E3) already yields that, for all sequences   →  and   ⇀  in , it holds ℐ(, ) ≤ lim inf →∞ ℐ(  ,   ). (2.1) Moreover, we assume that ℐ satisfies the following Gårding-like inequality: With respect to the time component we require The above assumptions combined with Gronwall's lemma allow us to obtain the following estimates that hold for all ,  ∈ [0,  ],  ∈ : and An energy functional that satisfies all assumptions made here is given in Section 4.
Remark 2.1.Note that we allow for an unbounded dissipation, which is essential for the application of our method to the damage model in Appendix B.
Combining the convexity and the positive 1-homogeneity of ℛ, it is easy to verify the following triangle inequality ℛ( − ) ≤ ℛ( − ) + ℛ( − ) ∀, ,  ∈ . (2.4) In fact we allow for an approximation of the "original" dissipation in our convergence analysis.In general, this corresponds to an approximation using e.g., finite elements.However, we will keep the setting general here, that is, we assume ℛ  :  → [0, ∞] satisfies the same assumptions as ℛ, i.e., (R1)-(R3), and Mosco-converges to ℛ w.r.t. the space  in the following sense: for all sequences   ⇀  in  (for  → 0) : lim inf Remark 2.2.Note that from now on, we consider ℛ and ℛ  , respectively, as mapping from  into [0, ∞].In fact, we will subsequently always evaluate ℛ at a point in .Thus, the space  is not used in the convergence analysis.Moreover, the choice ℛ  = ℛ (i.e., no additional approximation of ℛ) is clearly possible and fulfills all these assumptions.Another example that complies with all the assumptions (R1)-(R3) and (R4) but satisfies ℛ  ̸ = ℛ is given in Section 4.2 below.
Remark 2.3.In view of the previous Remark 2.2, for any  ∈  the subdifferential ℛ  () is subsequently considered as a subset of  * , i.e., by Lemma A.1 we have  ∈ ℛ  () iff In particular, in order to ease notation, we will refrain from using   ℛ  () keeping in mind that ℛ  () ⊂  * .Clearly, by the convexity and lower semicontinuity of ℛ  , the set ℛ  (0) is also weakly closed in  * .

Initial state
The initial value  0 is supposed to satisfy  0 ∈  and the local stability −  ℐ(0,  0 ) ∈ ℛ(0).Moreover, we assume that the approximations of the initial value as well satisfy −  ℐ(0,   0 ) ∈ ℛ  (0) for all  > 0 and that   0 is bounded in  independent of .

Definition of parametrized BV-solutions
We now turn to the actual definition of the so-called parametrized BV-solutions.As indicated in the introduction, this concept takes care of possible jump paths and relies on an energy identity.Definition 2.4.Let an initial value  0 ∈  be given.We call a tuple ( t, ẑ) parametrized BV-solution of (RIS), if there exists an artificial end time  ≥  such that the following conditions are satisfied: (ii) Initial and end time condition: where dist  * {, ℛ(0)} = inf{‖ − ‖  * :  ∈ ℛ(0)}, see also Lemma A.2.
We point out that it is always possible to rescale the artificial time in order to obtain a normalized parametrized BV-solution, where t′ () + ‖ẑ ′ ()‖  = 1 f.a.a. ∈ (0, ).The key idea here is to cut out all intervals where  ′ () + ‖ ′ ()‖  = 0 and to scale the artificial time appropriately, see, e.g., Lemma A.4.3 of [30].Moreover, let us mention that the regularity conditions in our definition of parametrized BV-solutions are chosen in such a way, that all terms contained are well-defined.Depending on the actual setting, particularly the choice of ℛ and ℐ, there might exist slightly different requirements, see, e.g., Definition 4.2 of [27].

Local stationarity scheme
The ultimate goal of this section is to prove that the subsequent algorithm, which is based on the local incremental minimization scheme (1.1), provides an approximation scheme for parametrized BV-solutions.The difference compared to (1.1) is that we search for stationary points of the constrained problem rather than global minima, see (alg 1 ).Thus, for a given time-discretization parameter 0 <  ≤  , the algorithm reads as follows: Algorithm (LISS).

3:
Compute a stationary point  ,  , i.e., with the indicator function   (see (3.1)), which, additionally, satisfies Time update: Set  →  + 1. 6: end while Note that the constraint in (1.1a) is incorporated into the algorithm through (alg 1 ) and the indicator function In comparison to (1.1), however, we do not use the "min" from (1.1b) in the time-update for purely technical reasons.In addition, the notation   is used here just once more to highlight that the subdifferential is in fact calculated in terms of the space , see also Remark 2.3.Nevertheless, the proposed method is closely related to (1.1), since a local minimizer of necessarily satisfies (alg 1 ).Moreover, thanks to the assumptions on ℐ and ℛ  , in particular weak lower semicontinuity, the existence of a global minimum of (3.2) and therefore also the existence of a stationary point fulfilling (alg 2 ) is guaranteed by the direct method in the calculus of variations.The reason for investigating (LISS) instead of (1.1), is the fact that a numerical algorithm for solving (1.1a) or rather (3.2) naturally provides a stationary point  ,  that satisfies ℐ )︁ but, in case of a nonconvex energy, is not guaranteed to be a global optimum of (1.1a) and (3.2), respectively.Moreover, since the concept of parametrized BV-solutions is based on a local stability condition, it is consistent to look for locally stable points, which are exactly the stationary points of (1.1a).Despite its necessity for the convergence analysis, the inequality in (alg 2 ) is also physically meaningful since it enforces the system to look for energetically preferable states, i.e., states with a lower energy cost.Concerning the exploration of this algorithm, particularly with a view to convergence, we proceed as follows: We start with characterizing properties of the stationary points.Afterwards, we turn to the essential a priori estimates that will allow a passage to the limit in the discrete version of the energy identity in (2.9), which is deduced in Lemma 3.10.The limit procedure itself is elaborated in the final Section 3.4.

Approximate discrete parametrized solution
The foundation for both, the a priori estimates and the discrete version of the energy identity, is given by the following Lemma 3.1.It provides various properties of a stationary point  ,  in (alg 1 ) and shows some similarities with the complementarity in (2.8).Indeed, we will see that one can interpret the stationarity condition as a discrete version of (2.8).Lemma 3.1 (Discrete optimality system).Let  ≥ 1 and  ,  be an arbitrary stationary point in the sense of (alg 1 ) with associated  , −1 given by (alg 3 ).Then the following properties are satisfied: There exists a subgradient  ,  ∈   ( ,  −  , −1 ) such that Herein, dist  * { •, ℛ  (0)} denotes the extended distance as defined in Lemma A.2.
Proof.The proof is given, e.g., in [22,30] under a slightly different setting.For convenience of the reader, we thus repeat the main steps.To shorten the notation we also set ℛ , = ℛ  +   , cf. (A.3), and suppress the superscripts ,  for the iterates throughout the proof.Thanks to a classical result of convex analysis, (alg 1 ) is equivalent to Moreover, from Lemma A.2, we infer Inserting this together with (3.5) in (3.4) gives (3.3c).
To prove (3.3a), we consider (alg 1 ) once more.Since 0 ∈ dom(ℛ  ) ∩ dom(  ) and   is continuous in 0, the sum rule for convex subdifferentials is applicable giving the existence of a   ∈   (  −  −1 ), such that and thereby A comparison with (3.3c), shows that Now, the fact that   ∈   (  −  −1 ), and the characterization in Lemma A.3 immediately yields (3.3a).
Remark 3.2.In fact, since (alg 1 ) is equivalent to the properties (3.3a)-(3.3d) it might be practical to exploit the characterization via (3.3a)-(3.3d)for the actual numerical realization of (LISS) instead of (alg 1 ) in order to calculate a stationary point.Moreover, we will solely build upon this discrete optimality system (and the inequality (alg 2 )) for the convergence analysis.
Let us take a further look at (3.3d).Exploiting the properties of  ,  from (A.6) it is easy to see that Combining this with (3.3a) and the time-update which is a discrete version of the complementarity condition in the definition of parametrized BV-solutions, cf. (2.8).

A priori estimates
Based on the previous Lemma 3.1, we subsequently provide several a priori estimates that will allow a passage to the limit in the discrete energy identity in Sections 3.3 and 3.4, respectively.Furthermore, we show that the discrete physical time  ,  given the time update in (alg 3 ) reaches the final time  in a finite number of iterations, see Proposition 3.7 below.We start with the following collection of results, whose proofs are basic so that we refer to [15,22] here.Let us, nevertheless, remark that these a priori estimates are the only point where one needs to exploit that  ,  is energetically preferred, that means (alg 2 ) holds, and not only a stationary point satisfying (3.3a)-(3.3d).
Lemma 3.3 (Boundedness for energy and dissipation).For all ,  > 0 and all  ∈ N, it holds sup where  and  are the components from Section 2.
Proof.The proof mainly relies on the estimates in (2.2) and (2.3) and the coercivity of the energy from assumption (E2), cf.[15,22].
Remark 3.4.As a consequence of Lemma 3.3 and the boundedness of   0 by assumption we have that  ,  ∈   (0, ) for some  > 0 independent of  and .
The estimate (3.12) will, on the one hand, provide us with a uniform  ∞ -bound for the linear interpolants of  and, on the other hand, allows us to obtain a bound for the derivative   ℐ.In preparation for that, we derive the following: Lemma 3.5.For every  > 0, there exists  1 (),  2 () > 0, such that for all ,  ∈   (0, ) and ,  ∈ [0,  ].
Proof.According to Ehrling's lemma, for every  > 0, there exists a constant   such that Combining the Gårding-like inequality from (E5) with (E6) we find To proceed, we consider each of the two last terms separately.For the first one, we exploit (3.14 where we used the lower bound for ℛ  from (R3) and the embedding  ˓→  in the last line.Next, we turn to the last term in (3.15).For this, we take advantage of Young's inequality which gives (3.17) Inserting (3.16) and (3.17) in (3.15) we eventually arrive at (3.13).

1231
Clearly, from the uniform boundedness of the iterates and ‖ +1 −   ‖  ≤  we can conclude the following result.
One major issue in the convergence analysis for parametrized BV-solutions concerns the boundedness of the artificial time, even in the continuous setting, see, e.g., the discussion in [24], p. 218.For the discrete counterpart, the artificial time reads In order to bound this term, we need to estimate ∑︀  =1 ‖  −  −1 ‖  , which is purpose of the next proposition.Moreover, we will show that the physical end time  is reached after a finite number of iterations, which guarantees that the algorithm finishes in a finite number of steps.Proposition 3.7 (Bound on artificial time).For every parameter ,  > 0 there exists an index  (, ) ∈ N such that  ,  (,) ≥  .Moreover, there are constants  1 ,  2 ,  3 > 0 independent of ,  such that for all ,  > 0 it holds (,) and dist Proof.The arguments are similar to [15,22].However, since there are some significant differences, particularly the estimate (3.20), we present the arguments in detail.Let  ∈ N be arbitrary.For convenience, we again suppress the superscript ,  throughout the proof, except for  , 0 in order to avoid confusion with the initial data.We start by testing (3.3d) with  =  +1 −   to obtain Inserting (3.3b) into (3.3c) and rewriting this identity for the index Thanks to the constraint Consequently, it holds Now, inserting the estimate from Corollary 3.6 gives Rearranging terms and summing up the resulting estimate with respect to , we arrive at Thanks to (3.11) it now suffices to estimate ⟨ 1 ,  1 −  0 ⟩  * , to prove (3.20), which is shown next.To this end, we again insert (3.3b) into (3.3c) to obtain for  = 1: .
Adding a zero, and rearranging terms yields By assumption, we have by the characterization in Lemma A.2, so that (3.24) implies For the first term on the left-hand side we take advantage of the assumption in (E5) and the boundedness of the iterates from Lemma 3.3 which results in Hence, using again the constraint ‖ +1 −   ‖  ≤  we obtain By adding (3.27) to (3.23) and applying (3.11), we arrive at where we used that   ≤  +  ≤ 2 by the time update in (alg 3 ) and  ≤  for the last estimate.Clearly, by the uniform boundedness of  , 0 by assumption and the continuity of ℐ(0, •), the term ℐ Hence, there must exist a finite index  (, ), possibly depending on  and , so that   (,) ≥  .Lastly, since In what follows we will abbreviate the index  (, ) simply by  having in mind that the number  of time steps always depends on  and .

Discrete energy-equality
In the following section, we aim at deriving a discrete analogue to the energy identity (2.9).To this end, we introduce the piecewise affine as well as the left-and right-continuous piecewise constant interpolants associated with the iterates  ,  .As indicated in the introduction, potential discontinuities of the parametrized BV-solution are resolved by introducing an artificial time.The physical time is accordingly interpreted as a function of the very same and jumps are characterized by the plateaus of this function.This is also reflected by the timeincremental stationarity scheme (LISS), where, loosely speaking, the artificial time is divided into equidistant subintervals with step size  and the approximation of the parametrized BV-solution is implicitly defined through the optimization in (LISS).To be more precise, we set  ,  :=  , so that Moreover, we define the artificial end time  , as that point where t reaches the end time  , i.e., it holds (see also Fig. 1) t, ( , ) = ,  Proof.The statements are a direct consequence of the properties in (3.10).
Proof.While the first three bounds are an immediate consequence of the results in Lemmas 3.8 and 3.3, the last one requires some slightly more explanation.

⃦ ⃦ ⃦
is missing in front of the distance.It would be possible to reformulate the optimality conditions in Lemma 3.1 in a way such that this coefficient would arise in (3.39).This, however, would complicate the passage to the limit in the next section.As we will see at the end of the proof of Theorem 3.14, equation (3.39) is sufficient to obtain the desired energy identity in (2.9).

Main convergence theorem
Before we come to the main result, i.e., the passage to the limit in the discrete energy identity and therewith ultimately the existence of parametrized BV-solutions, we need one last preparatory result, which guarantees the weak lower semicontinuity of the distance term in (3.39).Lemma 3.12.Let   ∈  * with   ⇀  in  * for  → 0. Suppose, moreover, that the distance is uniformly bounded, i.e., dist  * {−  , ℛ  (0)} ≤  with  independent of .Then the following weak lower semicontinuity result holds true: Since this holds for all subsequence of   , we ultimately arrive at the desired lower semicontinuity in (3.44).
We now have everything at hand to prove our main convergence result.Theorem 3.14 (Convergence towards parametrized BV-solutions).Assume that   0 converges to the initial state  0 in  for  → 0. Then there exists a sequence of parameters {  ,   } ∈N ⊂ R + × R + converging to zero so that the affine interpolants generated by the fully discrete local stationarity scheme (LISS) and the artificial end time defined in ( Proof.The arguments are analog to the ones in [22,30].For convenience of the reader, we briefly repeat the main steps. The existence of a (sub-)sequence satisfying (3.48)-(3.50) is an immediate consequence of the uniform estimates in Lemmas 3.3, 3.9, and (3.34).The pointwise convergence in (3.51) follows from the Aubin-Lions lemma, i.e.,  1,∞ (0, ; ) ∩  ∞ (0, ; ) ˓→  (0, ; ), the density of  in  and the fact that for every  ∈ [0, ], {ẑ , ()} ∈N is bounded in  by Lemma 3.3.
It remains to show that every (weak) limit is a parametrized BV-solution.For this purpose, let {  ,   } be an arbitrary null sequence and assume that the convergences in (3.48)-(3.51)hold.In order to simplify the notation, we indicate by {•}  the sequence of {•} , corresponding to {  ,   }.Analogously, we abbreviate the index   simply by .We proceed in several steps and start with the following: Convergence of piecewise constant interpolants.One easily verifies using the estimate which holds for all  ∈ {1, ...,  } and all  ∈ [︀   −1 ,    )︀ that the piecewise constant interpolants converge pointwise to the same limit as the affine interpolants.We therefore have Initial and end time conditions.By assumption we have ẑ (0) =   0 →  0 in , so that the pointwise convergence in (3.51) implies ẑ(0) =  0 as desired.Moreover, thanks to (3.49), t converges uniformly to t so that 0 = t (0) → t(0) and  = t (  ) → t(), where we also used (3.48).
Energy identity.The energy identity is a direct consequence of its discrete version in Lemma 3.10.Indeed, we find by the weak lower semicontinuity of ℐ(, •) from (2.1) and Corollary 4.5 of [31] which gives due to the assumptions on the space  and condition (R4) on ℛ  .Moreover, we used ‖ẑ ′ ()‖  ≤ 1 as well as Fatou's lemma together with (3.53) for the distance term.Now, exploiting the estimate in (3.42), assumption (E4) for   ℐ and the strong convergence of ẑ (0) to  0 in  we finally end up with which is the desired energy inequality.Taking into account that ẑ ∈  1 (0, ; ), it is well-known that the sole inequality (3.55) is already equivalent to the energy identity (2.9), see Lemma 6.6 of [16] or Lemma 2.4.6 of [30], which completes the proof.
Unfortunately, we do not obtain the nondegeneracy let alone normalization of the limit (︀ t, ẑ)︀ here.The main problem is the fact that the weak convergence of ẑ in  1 (0, ; ) from (3.50) is not sufficient in order to pass to the limit in (3.36), that is, t′  () + ‖ẑ ′  ()‖  = 1, and still obtain equality in the end.In [9,25], the authors therefore provide sufficient conditions, which guarantee the nondegeneracy of the limit function.Moreover, in [9], a condition is given, which also preserves the normalization.Nevertheless it is always possible to reparametrize a parametrized BV-solution in order to normalize it, see Lemma A.4.3 of [30].Regardless of this fact, we note that the above theorem, while dedicated to the convergence analysis of the fully discrete local stationarity scheme, also provides an existence result for parametrized BV-solutions in case of an unbounded dissipation ℛ (choose

Application to a damage model
We now aim at applying the local stationarity scheme to model the evolution of damage within a workpiece during a time interval [0,  ].For this, we let Ω ⊂ R 2 be a bounded domain that corresponds to an elastic body and satisfies Ω ⊂ R 2 has a Lipschitz boundary Ω = Γ  ∪ Γ  with Dirichlet boundary Γ  such that ℋ 1 (Γ  ) > 0 and Neumann boundary Γ  .Moreover, Γ  and Γ  are supposed to be regular in the sense of Grger, see [12].
Note that, although we focus on the two-dimensional case here, it is also possible to consider the threedimensional case as well provided the spaces and the energy are adapted appropriately; compare with the elaborations in [16,23] which we also follow with regard to notation.During the time [0,  ], time dependent boundary conditions   as well as external boundary and volume forces ℓ may be applied, which lead to a certain displacement  and possibly even to a damage, represented by the variable , of the body.Usually,  is supposed to take values in [0, 1] whereby (, ) = 1 means the body is completely sound and, correspondingly, (, ) = 0 means the body is completely damaged.With a view to the energy functional, we define where  > 2 is chosen as in Lemma 4.1 below.Note that we will use a scaled version of the  2 -norm, that is This choice has been shown to be advantageous in the numerical experiments particularly with a view to iteration numbers.Now, we set ℰ : [0,  ] ×  ×  → R as where C is the usual elasticity tensor with ∃ 0 > 0 such that for all  ∈ R 2×2  and for almost all  ∈ Ω : and is the linearized strain tensor.The nonlinearity  in the energy is chosen such that the term is continuously differentiable in .Since we will neglect this term in our numerical examples, we do not particularize the exact assumptions and refer to [16].However, a typical choice for  in the context of Ambrosio-Tortorelli approximation of brittle fracture is  () = (1 − ) 2 , see [2,11], which is clearly sufficiently smooth.Note at this point that the meaning of the damage variable  is in some sense inverse to the one used in [16], i.e., (, ) = 0 here means the material is completely sound, while (, ) → ∞ if the material gets completely damaged.Furthermore, the function , which somehow represents the preservation of the elasticity of the material depending on the state of damage, is supposed to fulfill: In particular, the lower bound  ≥  1 > 0 is to be noted here.It implies that even if the material is completely damaged, it does not lose all its rigidity.This is often referred to as partial damage model.Finally, the dissipation with the so-called fracture toughness  > 0. Now, in order to bring this model into the setting of Section 2, it is convenient to reduce the system to the damage variable .This means that we require the displacement () to minimize the energy ℰ(, •, ()) at every time point  ∈ [0,  ], i.e., () ∈ arg min{ℰ(, , ) :  ∈ }.
It is, in fact, easy to see that this problem has a unique minimizer for every  ∈ [0,  ] and  ∈ .Hence, we define

Properties of the energy functional
We now want to verify that the model from above fits into the setting of Section 2. Thereby, we rely on the results from [16,23].We start with the following observation, which was first proven in [13] and states that the minimization with respect to  is well-defined and provides a unique solution.
Therefore, the reduced energy ℐ 2 (, ) = inf ∈ ℰ 2 (, , ) is also well-defined and we can focus on the properties of this part in the overall energy ℐ.Lemma 4.2.Let  > 2 and the assumptions (4.1), (4.2) hold.Then there exist constants  1 ,  2 ,  3 > 0 such that as well as ) , where  > 2 is as in Lemma 4.1.Moreover, for any sequences   →  and   ⇀  in , it holds Proof.This is a combination of Lemma 2.4, 2.6 and 2.8 as well as Corollary 2.9 from [16].
Proof.The conditions (E1)-(E4) and (E6), (E7) follow immediately from the above Lemma 4.2.In addition, the Gårding-like inequality (E5) is an easy consequence of the properties of  and the inequality in (4.7).Thus, we see that ℐ fulfills all assumptions from Section 2 so that applying Theorem 3.14 proofs the existence of a parametrized BV-solution.
Theorem 3.14 guarantees the existence of at least one parametrized BV-solution for (RIS) in the setting of (partial) damage here.This particularly applies to the two dimensional case with the standard Laplacian as regularizer for the damage variable which as also been investigated in [16].In comparison, however, we slightly reduced the assumption on the initial state, that is, we do not require that −  ℐ(0,  0 ) ∈  but −  ℐ(0,  0 ) ∈ ℛ(0) ⊂  * .What is more, we incorporated the possibility of using a space-discretization via the approximation of the dissipation ℛ.Hence, we may approximate a parametrized BV-solution for the damage model by the local incremental stationarity scheme (LISS), which is purpose of the following subsections.

Finite element discretization
As the convergence analysis from Theorem 3.14 allows us to use an approximation ℛ  of the dissipation ℛ, we may use some finite element discretization to approximate a parametrized BV-solution.Hence, we assume that a family { ℎ } ℎ>0 of shape-regular triangulations of the domain Ω be given.Herein, ℎ denotes the mesh size defined by ℎ := max  ∈ ℎ diam( ).To keep the discussion concise, we also assume that Ω is a polygon and that the triangulations exactly fit the boundary.For the discrete space we choose the space of piecewise linear and continuous test functions, i.e., In addition, we set By standard arguments the lower inequality in (R4) is satisfied for  = ℎ → 0. For the upper inequality assume that ℛ() < ∞, i.e.,  ≥ 0 a.e. in Ω (otherwise there is nothing to show).Since Ω has a Lipschitz-boundary we can extend  beyond Ω (cf.[3], A8.12) and use standard convolution in order to obtain an approximation . By the construction of the extension and the non-negativity of the convolution kernel,   is also non-negative.Moreover, we have ‖  − ‖  → 0 for  → ∞.For   ∈  ∞ (︀ Ω )︀ the classical, pointwise Lagrangeinterpolation  ℎ is well-defined so that ‖  −  ℎ (  )‖  → 0 for ℎ → 0. Hence, for any  ∈ N, there exists ℎ  > 0 with ‖  −  ℎ (  )‖  ≤ 1/ for all ℎ ≤ ℎ  .W.l.o.g.we may assume that {ℎ  } ∈N is strictly monotonic decreasing.Therewith, we define  ℎ :=  ℎ (  ) for ℎ  ≥ ℎ > ℎ +1 .Combining the above properties, we find that  ℎ ≥ 0 a.e. in Ω with  ℎ →  in  for ℎ → 0 and by the continuity of ℛ on the set of non-negative functions in .

Discrete energy functional
In analogy to  we denote by  = ( With this operator at hand, we can also reduce the discrete energy Ẽ to the discrete damage variable  ℎ .Hence, we define Ĩ :  ℎ → R as Ĩ(,  ℎ ) = Ẽ(,  ℎ ,  ℎ ( ℎ )).Altogether, with a little abuse of notation, we denote the reduced energy functional considered as mapping acting on the coefficient vector by  : R  ℎ → R.

Numerical solution of the local minimization problems
With all the notations above, particularly the description of ℛ ℎ in (4.10), the stationary equation (alg 1 ) is equivalent to the following problem for the coefficient vector   : Here and for the rest of this section, we abbreviate  , −1 simply by  −1 as well as  = 1 |Ω| .Inserting the characterizations of   and   and taking )︁ , we therefore find that which can be equivalently formulated as

Numerical results
For our numerical tests, we use two benchmark tests from [8].In any of these cases we set the external volume and surface forces to zero, so that ℓ ≡ 0, and the softening function  to () = exp(−) +  with  = 0.01.The numerical computations are performed with Matlab c ○ and the linear systems of equations arising in each semismooth Newton step (cf.Appendix B) are solved by Matlab's inbuilt direct solver based on UMFPACK.The implementation for the elasticity part of the energy relies on the Matlab code from [1].
Example I: Pre-cracked brick The geometry for this example is shown in Figure 2. Due to its symmetry the computation is performed using only one quarter of the whole system.Therefore, the symmetry axes become parts of the boundary of Ω and we impose the following symmetry boundary conditions  1 = 0 on Γ 1 = {0} × [16,40] and  2 = 0 on Γ 2 = [0, 100] × {40}.
During the time interval [0, 16] the workpiece is stretched at its ends, which we realize by Dirichlet conditions on Γ  , i.e., Furthermore, we use the parameters given in the table of Figure 2 (right).We initiate the evolution with  0 ≡ 0 and choose  = 0.1 for the time step size.The state of damage at several points during the time interval [0, 16] are shown in Figure 3.
Obviously, damage occurs at first at the tip of the crack and evolves along the symmetry axis Γ 1 afterwards.As expected, it also concentrates on regions with large stresses.However, the sharpness of interfaces between damaged and undamaged areas highly depends on the choice of the functional ℐ 1 , see [6].In our case, this interface is rather diffuse and cannot be sharpened by a refinement of the mesh, see Figure 4. Clearly, one may reduce the factor  included in the operator .However, this leads to instabilities in the semismooth Newton method and globalization strategies might be necessary, which is subject of future research.Nevertheless, the results are stable with respect to mesh refinement, which can be observed in Figure 5 that shows the forcedisplacement-diagram for three different mesh sizes.The reaction force is calculated by integrating the stress  in normal direction along the boundary Γ  .In comparison with the results from [8], there is a high degree of conformity, except for the larger values of the reaction force here after reaching its maximum at approximately 0.08 mm of displacement, cf. also with Example II.With regard to the time function t that is depicted in Figure 6 we see that the spreading of the damage area starting at  ≈ 8 is slightly faster than the rest of the evolution but does not cause a jump.
Let us finally compare the results obtained with an alternative approach.Instead of incorporating the constraint ‖  −  −1 ‖  ≤  which corresponds to an  2 -ball of size  it is convenient to take ‖  −  −1 ‖  ∞ ≤  .+  as a constraint since this leads to pointwise box constraints that are particularly well suited to semismooth Newton-methods, see e.g., [14].Unfortunately, this choice is not covered by our convergence analysis for the local stationarity scheme (LISS).Nevertheless, we provide the numerical results that can be obtained by using this version.Thus, while we change (4.15) to       we keep all parameters unaltered.Clearly, the Newton-matrix   also changes in the obvious way.The combined force-displacement-diagram for both solutions using LISS is given in Figure 7 while Figure 8 shows the timefunction t for the solution using box-constraints.At first glance the plot of the time-function seems to indicates that there is a larger deviation for the two solutions around the time  ≈ 8.However, it is easy to see that at least the reaction forces for both solutions are very similar.What is more, the  2 -distance at the end time It is easy to see that both solutions, using  2and  ∞ -constraints, are very similar.Indeed, the  2 and  ∞ distance at the end time  = 16 calculates to Since the incorporation of box constraints is natural in this context and its implementation is also easier to realize it provides an interesting topic for further research.
A possible approach for this is to consider the -Laplacian with  >  for the regularization instead of  = −∆ as proposed in [17].In this case, we have which may open up the opportunity to use the  ∞ -norm as constraint in (3.2).

Example II: Brick with a circular hole
As proof-of-concept we consider the situation depicted in Figure 9. Again, due to its symmetry the computation is performed using only one quarter of the whole system.This also implies certain symmetry conditions for the boundary of the domain Ω.For the situation at hand, we impose  1 = 0 on Γ 1 = {0} × [50, 100] and  2 = 0 on Γ 2 = [50, 100] × {0}.
The parameters used are given in the table of Figure 9 (right).Finally, we initiate the evolution with  0 ≡ 0 and choose  = 0.1 for the time step size.The numerical solution of the damage variable for an intermediate time point  ≈ 7.8 is shown in Figure 10 and the corresponding force-displacement-diagram is depicted in Figure 11.We observe a strong similarity with the results in [8].However, the reaction force in the end phase of the evolution is significantly higher in our case.This may be a result of the additional regularization by the introduction of a "local field"  in [8].Certainly, this requires further investigation.

Conclusions
We presented a numerical scheme for the approximation of parametrized BV-solutions for rate-independent systems including a fairly general setting for the energy and dissipation functionals.The scheme itself is based on the local minimization scheme introduced in [9], but relies on stationary points rather than local minima, making it very accessible for numerical optimization algorithms (limit points are, in general, stationary).Moreover, by adapting the convergence analysis of the recent contributions in [15,22] and using arguments from [16], we proved that the scheme provides parametrized BV-solutions of the original rate-independent system under Mosco-convergence of approximations for the dissipation ℛ.While this is at first glance a result that verifies the consistency of the local incremental stationarity scheme, we, moreover, gain an existence result for parametrized BV-solutions in the case of a nonconvex energy and unbounded dissipation.We then focused on the realization of our scheme for a model of the evolution of damage within a workpiece.We employed a finite element discretization in space and used a semismooth Newton method for solving the discrete stationary system arising in each step of the scheme.The resulting algorithm behaves efficient and robust in our numerical tests.Eventually there are several topics for future research.This concerns for example the usage of an  ∞ -norm in the indicator functional in (alg 1 ) which leads to an easy to implement algorithm for the damage model considered in this paper.In the same context, it might be interesting to relax the assumptions for the energy in order to incorporate functionals that allow for a sharper resolution of the interface between damaged and undamaged regions.Moreover, considering dissipation functionals which are also depending on the state , i.e., ℛ = ℛ(,  ′ ) should be noted here.As there are only few results in this direction for rate-independent systems this does not only concern parametrized BV-solutions.with the iterate   = (  ,   ,   ) and a Newton-derivative    .However, the matrix   contains second order information of  which, by (B.3) necessitates the determination of  ′ ℎ ( ℎ ).In order to keep track of this, we blow up the whole system so that in every semismooth Newton-step we actually solve Eventually, we update  +1 =   + ∆  ,  +1 =   + ∆  and  +1 =   + ∆  .With this choice, all matrices H appearing in our numerical test have shown to be invertible and the semismoothNewton method performed well with respect to both, robustness and efficiency.In particular, no globalization efforts are needed to ensure convergence of the method and, moreover, no line search was necessary in order to guarantee condition (alg 2 ) in LISS.A rigorous convergence analysis of the method however would go beyond the scope of this paper and is subject to future research.

(3. 31 )Figure 1 .
Figure 1.Qualitative illustration of the affine interpolant t, the choice of the artificial end time  , via the equality t( , ) =  and the upper bound S.

Figure 2 .
Figure 2. Geometry of the domain (left); Table of parameters (right).

Figure 3 .
Figure 3. State of the dimensionless damage variable  ℎ at different points in time.Note that the scale of the colormap varies with the time.(A)  ≈ 4.9.(B)  ≈ 7.88.(C)  ≈ 8.74.(D)  ≈ 16.

Figure 4 .
Figure 4. Refined and deformed mesh combined with the state of the dimensionless damage variable at the final time  = 16.Note that the displacement is magnified by a factor of 20.

Figure 5 .
Figure 5. Reaction Force on the boundary Γ  depending on the displacement   for three different mesh sizes.

Figure 6 .
Figure 6.Function t in dependence of the artificial time .

Figure 7 .
Figure 7.Comparison of the reaction force for the solution obtained by using the  2constraints (red) and box-constraints (blue).

Figure 8 .
Figure 8. Function t in dependence of the artificial time  for the solution obtained by using box-constraints.

Figure 9 .
Figure 9. Geometry of the domain (left); Table of parameters (right).

Figure 10 .
Figure 10.State of the dimensionless damage function  ℎ at time point  ≈ 7.8

Figure 11 .
Figure 11.Reaction force on the boundary Γ  depending on the displacement   .