Blaise Bourdin

Gilles A Francfort

Fracture

1 From Griffith to the variational In this section, the starting premise is Griffith's model for crack evolution, as presented in [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. Of course, continuum mechanics has evolved and it would make little sense to present fracture exactly as in [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. The reader will find below what we believe to be a very classical introduction to brittle fracture within a rational mechanical framework.

Our starting assumptions are two-fold. First we restrict our focus to quasi-static evolution, a huge restriction: At each time, the investigated sample is in static equilibrium with the loads that are applied to it at that time. We use the blanket label "loads" for both hard devices (displacement type boundary conditions) and soft devices (traction type boundary conditions and/or body forces). Then, we do not concern ourselves with changes in temperature, implicitly assuming that those will not impact upon the 1 mechanics of the evolution: in particular, thermal expansion is not covered by this model, at least to the extent that it couples thermal and mechanical effects. However, thermal stresses induced by a known temperature field fall squarely within the scope of the forthcoming analysis.

Also, we only discuss the 2d-case in this section. However, it will be clear that the resulting formulation applies as well to dimensions 1 and 3.

We consider Ω, a bounded open domain of R 2 . That domain is filled with a brittle elastic material. At this level of generality, the type of elastic behavior matters little, as long as it is represented by a bulk energy F → W (F ) which will be assumed to be a function of the gradient of the deformation field ϕ; in linearized elasticity W will become a function of e(u) := 1 2 (∇u + ∇u t ) with ϕ(x) = x + u(x). Time dependent loads are applied to Ω. We will assume that the force part of the load is given in the reference configuration (that is defined on Ω). Those are

• body forces denoted by f b (t) and defined on Ω;

• surface forces denoted by f s (t) and defined on ∂ s Ω ⊂ ∂Ω;

• boundary displacements denoted by g(t) and defined on ∂ d Ω := ∂Ω \ ∂ s Ω. Precisely, we assume throughout that g(t) is defined and smooth enough on all of R 2 and that the boundary displacement is the trace of g(t) on ∂ d Ω.

Griffith's theory

Griffith's theory is purely macroscopic. The crack or cracks are discontinuity surfaces for the deformation field of the continuum under investigation. If that continuum behaves elastically, material response under external loading will be unambiguous once the laws that preside over the onset and propagation of the crack(s) are specified. The construction of such lawsthe goal of Griffith's theory -requires three foundational ingredients, 1. A surface energy associated to the surfaces where the deformation is discontinuous; 2. A propagation criterion for those surfaces; 3. An irreversibility condition for the cracking process. The surface energy adopted by Griffith is simple. Throughout the cracking process, a(n isotropic) homogeneous material spends an energy amount which remains proportional to the area of the surface of discontinuity. We call fracture toughness of the material the proportionality factor, and denote it by k.

A simple counting argument demonstrates that, if inter-atomic bonding is ruled by a Lennard-Jones type interaction potential, then the add-energy spent in moving two atoms apart while the remaining atoms stay put is additive, which ultimately yields a total (macroscopic) energy proportional to the separation area. In the absence of contact the crack lips do not interact and cohesiveness is prohibited.

The propagation criterion is energy based. The test is a balance between the potential energy released through a virtual increase of the crack length (area) and the energy spent in creating additional length (area). The crack will extend only if the balance favors creation.

Finally a crack will form where and at the time at which the displacement field becomes discontinuous. It will then stay so forever, oblivious to the actual state of displacement at any posterior time.

We now formulate Griffith's view of the crack evolution problem in a(n isotropic) homogeneous elastic material. For now, the crack path Γ is assumed to be known a priori. We wish to include partial debonding as a possible crack behavior, so that Γ ⊂ Ω \ ∂ s Ω. The crack at time t is assumed to be a time increasing connected subset of Γ; it can thus live partially, or totally on ∂Ω. It is therefore completely determined by its length l and denoted by Γ(l).

By the quasi-static assumption, the cracked solid (see Figure 1.1.1) is, at each time, in elastic equilibrium with the loads that it supports at that time; in other words, if the crack length at that time is l, then the kinematic unknown at that time, ϕ(t, l) (the transformation, or displacement) satisfies where n denotes the appropriate normal vector.

The last relation in (1.1.1) calls for several comments. In an anti-plane shear setting, it merely states, in accord with Griffith's premise, the absence of cohesive forces along the crack lips. In a planar situation, it implicitly assumes separation of the crack lips, hence non-interpenetration.

The system (1.1.1) assumes that the crack length is known. Griffith's contribution is to propose the following criteria for the determination of that length. At time t, compute the potential energy associated to the crack of length l, that is

P(t, l) := Ω\Γ(l)
W (∇ϕ(t, l)) dx -F(t, ϕ(t, l))

(1.1.2) with

F(t, ϕ) := Ω f b (t).ϕ dx + ∂sΩ f s (t).ϕ ds. (1.1.3)
Then, l(t) must be such that it obeys

• The Griffith's criterion:

i. l t ր (the crack can only grow);

ii. -∂P ∂l (t, l(t)) ≤ k (the energy release rate is bounded from above by the fracture toughness);

iii. ∂P ∂l (t, l(t)) + k l(t) = 0 (the crack will not grow unless the energy release rate is critical).

Remark 1.1.1 From a thermodynamical viewpoint, Griffith's criterion should be interpreted as follows. The crack length is a global internal variable, and its variation induces a dissipation which must in turn satisfy Clausius-Duhem's inequality.

A convenient enforcement of Clausius-Duhem's inequality is provided through the introduction of a convex dissipation potential D( l), further satisfying D(0) = 0. Then, the inequality reduces to -∂P ∂l (t, l(t)) ∈ ∂D( l(t)).

(1.1.4)

The correct dissipation potential in Griffith's setting is denoted by D G and given by (see Figure 1.1.2)

D G ( l) := k l, l ≥ 0 ∞, l < 0, (1.1.5)
and (1.1.4) then yields precisely Griffith's criteria. So, summing up, Griffith's modeling of crack evolution reduces to (1.1.1), (1.1.4) with (1.1.5) as dissipation potential. Note that the the dissipation potential is positively 1-homogeneous. This is an essential feature if adopting a variational viewpoint.

Irreversibility forbidden region

We are now ready to explore the system (1.1.1), (1.1.4). For completeness, we should add an initial condition to (1.1.4); we will thus assume that l(0) = l 0 , (1.1.6) and denote, from now onward, any pair-solution (l(t), ϕ(t, l(t))), if it exists, by (l(t), ϕ(t)).

A variational equivalence

Assuming suitable -and unstated -smoothness of all relevant quantities, we propose to establish the equivalence between the original system (1.1.1), (1.1.4), (1.1.6) and a formulation which states that a certain energy must remain stationary at every time among all virtual admissible crackdisplacement pairs at that time, and that an energy conservation statement must be satisfied throughout the time evolution. This is the object of the following Proposition 1.2.1 The pair (l(t), ϕ(t)) (satisfying (1.1.6)) satisfies (1.1.1), (1.1.4) (with appropriate smoothness) on [0, T ] iff, for every t ∈ [0, T ], it satisfies (with that same smoothness) (Ust) (l(t), ϕ(t)) is a stationary point -in the sense of (1.2.5) below -of

E(t; ϕ, l) := Ω\Γ(l) W (∇ϕ) dx -F(t, ϕ) + kl, (1.2.1) among all l ≥ l(t) and ϕ = g(t) on ∂ d Ω \ Γ(l); (Ir) l(t) ≥ 0 ; (1.2.2) (Eb) dE dt (t) = ∂ d Ω\Γ(l(t)) ∂W ∂F (∇ϕ(t))n. ġ(t) ds -Ḟ(t, ϕ(t))
with

Ḟ(t, ϕ) := Ω ḟb (t).ϕ dx + ∂sΩ ḟs (t).ϕ ds (1.2.3) E(t) := Ω\Γ(l(t))
W (∇ϕ(t)) dx -F(t, ϕ(t)) + kl(t)

= P(t, l(t)) + kl(t).

(1.2.4)

The unilateral stationarity statement (Ust) is rather unusual because the functional E(t; •) that should be stationary at (l(t), ϕ(t)) explicitly depends on l(t); hence the label unilateral. The energy balance (Eb) can be turned, through various integration by parts in time, into what is referred to in the literature as the mechanical form of the second law of thermodynamics; see e.g. [START_REF] Gurtin | Configurational forces as basic concepts of continuum physics[END_REF].

Proof. First we should clearly articulate what is meant by (Ust). To this effect, we introduce a one-parameter family of variations of the kinematic variable ϕ(t) and of the crack length l(t) as follows. We set

l(t, ε) := l(t) + ε l; l ≥ 0; ϕ(t, ε, l) := ϕ(t, l) + εψ(t, l),
where ψ(t, l) = 0 on ∂ d Ω \ Γ(l) and ϕ(t, l(t)) = ϕ(t). Then, unilateral stationarity is meant as

d dε E(t; ϕ(t, ε, l(t, ε)), l(t, ε)) ε=0 ≥ 0. (1.2.5)
Recall the expression (1.2.1) for E. Then, the above also reads as

Ω\Γ(l(t)) ∂W ∂F (∇ϕ(t)).∇ψ dx -F(t, ψ) + ∂P ∂l (t, l(t)) l + k l ≥ 0,
where we recall that P was defined in (1.1.2). Consequently, through integration by parts, (Ust) is equivalent to

(1.1.1) and ∂P ∂l (t, l(t)) + k ≥ 0. (1.2.6)
Then, assume that (Ust), (Ir), (Eb) hold. In view of the above, (1.1.1) is satisfied, so that (Eb) reduces to

( ∂P ∂l (t, l(t)) + k) l = 0. (1.2.7) Conversely, if (1.1.1) holds true, then dE dt (t) = ∂ d Ω\Γ(l(t)) ∂W ∂F (∇ϕ(t)). ġ(t) ds -Ḟ(t, ϕ(t)) + { ∂P ∂l (t, l(t)) + k) l(t)}, (1.2.8)
and, in view of the third item in Griffith's criterion, the term in brackets in (1.2.8) cancels out and (Eb) is established.

The second item in Griffith's criterion, together with (1.1.1), implies (Ust).

Remark 1.2.2 Elimination of the kinematic field in the variational formulation leads to the sometimes more convenient equivalent formulation for (Ust). (Ust) l(t) is a stationary point of P(t, l) + kl, among all l ≥ l(t).

Modulo smoothness, Griffith's formulation and the variational formulation obtained in Proposition 1.2.1 are strictly one and the same and cannot be opposed on mechanical grounds anymore than the original formulation. Pre-assuming smoothness is universal practice in deriving a notion of weak solution, so that we feel perfectly justified in doing so, and will be quite qualify as "weak" the solutions of what we will, from now onward, label the "variational evolution".

Remark 1.2.3 Consider the case of a p > 1-homogeneous elastic energy density and of a monotonically increasing load, that is

W (tF ) = t p W (F ), F(t, tu) = t p F(1, u), g(t) = t g,
where u is the displacement field. Then, by homogeneity, u(t, l) = tū(ℓ), P(t, l) = t p P(l),

where P(l) is the potential energy associated to a crack of length l and loads corresponding to the value t = 1. We assume that P is a sufficiently smooth function of l. Then it can be shown that that energy is a strictly convex function of l on [l 0 , l 1 ], if, and only if Griffith's criterion is satisfied by a unique smooth crack propagation l(t) on [t 0 , t 1 ] given by

l(t) = ( P′ ) -1 - k t p , t 1 = p k -P′ (l 1 )
.

(1.2.9)

Then, at each time t, -t p P′ (l(t)) = k. Thus, smoothness of the propagation leads to a reinforcement of the unilateral stationarity principle (Ust). The crack length l(t) must actually be a minimizer for P(t, l) + kl, because of the necessary convexity of P. So Griffith's criterion, because it assumes smoothness of the crack evolutions, implicitly pre-supposes the global convexity of the potential energy as a function of the crack length.

Stationarity is not a very pleasant mathematical notion from the standpoint of existence and it it is tempting to somewhat strengthen (Ust). Observe that (Ust) amounts to a first order optimality condition for (l(t), ϕ(t))) to be a local unilateral minimizer -in any reasonable topology -of E(t; •).

The preceding analysis suggests the adoption of some kind of minimality principle. Consequently, we propose the following two levels of departure from Griffith's classical theory:

• Local level -(Ust) is replaced by (Ulm) (l(t), ϕ(t)
) is a local minimizer (in a topology that remains to be specified) for E(t; ϕ, l) among all l ≥ l(t) and all ϕ = g(t) on ∂ d Ω\Γ(l);

• Global level -(Ust) is replaced by (Ugm) (l(t), ϕ(t)) is a global minimizer for E(t; ϕ, l) among all l ≥ l(t)
and all ϕ = g(t) on ∂ d Ω \ Γ(l). Those criteria are common in mechanics, but can never be justified on mechanical grounds, at least when departing from a purely convex setting. Local minimality is equivalent to Lyapunov stability for systems with a finite number of degrees of freedom, while the modern treatment of finite elasticity usually resorts to global minimality of the potential energy. Of course, in our setting, the minimization criterion, be it global or local, must also accommodate irreversibility, hence the already mentioned notion of unilaterality.

We now return to the time-continuous variational evolution and recast it in a more suitable functional framework.

Functional framework-A weak variational evolution

Using a minimality criterion immediately frees the crack path. Indeed, the minimality-modified Griffith variational evolution states that the actual length l(t) of the crack is a local (or global) minimum among all lengths l greater than, or equal to l(t) along the pre-determined crack path Γ. There is then no point in restricting the future evolution precisely to the curve Γ. We may as well let the crack choose which future path it wishes to borrow, according to the minimality principle. 1 Thus, denoting by Γ(t) the crack at time t, we replace (Ulm), resp. (Ugm) by (Ulm) (Γ(t), ϕ(t)) is a local minimizer (in a topology that remains to be specified) for

E(t; ϕ, Γ) := Ω\Γ W (∇ϕ) dx -F(t, ϕ) + kH 1 (Γ), (1.3.1) among all Γ ⊃ Γ(t) and all ϕ = g(t) on ∂ d Ω \ Γ; or, resp., (Ugm) (Γ(t), ϕ(t)) is a global minimizer for E(t; ϕ, Γ) among all Γ ⊃ Γ(t
E(t) := Ω\Γ(t) W (∇ϕ(t)) dx -F(t, ϕ(t)) + kH 1 (Γ(t)) = P(t, Γ(t)) +kH 1 (Γ(t)), (1.3.2)
1 It is precisely that freedom which truly distinguishes our approach from those usually adopted in mechanics; it is also precisely that freedom which confuses the mechanician, enslaved from a very early stage to a preconceived notion of what the crack can or cannot do.

with an obvious extension of the definition (1.1.2) of the potential energy P.

We allow the test cracks Γ to be pretty much any closed set in Ω \ ∂ s Ω with finite Hausdorff measure. This allows us to envision very rough cracks, with length that coincide with the usual length when the crack is a rectifiable curve. We do not allow for the crack to lie on ∂ s Ω because the crack cannot live where soft devices are applied.

We shall refer to the above formulation, that is (Ic), (Ulm) or (Ugm), (Ir), (Eb), as the strong variational evolution.

Local minimality directly refers to a topology, whereas global minimality is topology-independent. But, even then, we need to impart upon test cracks a decent topology. A natural candidate is the Hausdorff metric, defined for two closed sets A, B as

0 = ∅, f b (0) = f s (0) = 0. Then, we should minimize Ω\Γ W (∇ϕ) dx + kH 1 (Γ)
among all pairs (Γ, ϕ) with ϕ = g(0) on ∂ d Ω \ Γ. The direct method of the calculus of variations would have us take an infimizing sequence {(Γ n , ϕ n )}. In particular, we are at liberty to assume that H 1 (Γ n ) ≤ C. Say that the sequence Γ n converges in the Hausdorff metric to some Γ; this is not a restriction, thanks to Blaschke's compactness theorem [START_REF] Rogers | Hausdorff measures[END_REF]. Then we would like to have

H 1 (Γ) ≤ lim inf n H 1 (Γ n ).
But, this is generically false, except in 2d and for, say, connected Γ n 's! That topology has been used with success to prove existence, in the global minimality framework, for the 2d variational evolution restricted to connected cracks in Dal [START_REF] Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF]. We shall come back to this point in Section 3.

In the context of image segmentation, D. Mumford and J. Shah proposed to segment image through the following algorithm: Find a pair K, compact of Ω ⊂ R 2 (the picture) representing the contours of the image in the picture, and ϕ, the true pixel intensity at each point of the picture, an element of C 1 (Ω \ K), which minimizes

Ω\K |∇ϕ| 2 dx + kH 1 (K) + Ω |ϕ -g| 2 dx, (1.3.3)
where g is the measured pixel intensity. The minimization proposed in [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] was then shown in De [START_REF] De Giorgi | Existence theorem for a minimum problem with free discontinuity set[END_REF] to be equivalent to a well-posed one-field minimization problem on a subspace SBV (Ω) of the space BV (Ω) of functions with bounded variations on Ω, namely,

Ω |∇ϕ| 2 dx + kH 1 (S(ϕ)) + Ω |ϕ -g| 2 dx, (1.3.4)
where ∇ϕ represents the absolutely continuous part of the weak derivative of ϕ (a measure), and S(ϕ) the set of jump points for ϕ.

We recall that a function ϕ : Ω → R is in BV (Ω) iff ϕ ∈ L 1 (Ω) and its distributional derivative Dϕ is a measure with bounded total variation. Then, the theory developed by E. De Giorgi (see e.g. [START_REF] Evans | Measure theory and fine properties of functions[END_REF]) implies that

Dϕ = ∇ϕ(x) dx + (ϕ + (x) -ϕ -(x))ν(x)H 1 S(ϕ) + C(ϕ),
with ∇ϕ, the approximate gradient, ∈ L 1 (Ω) (∇ϕ is no longer a gradient), S(ϕ) the complement of the set of Lebesgue points of ϕ, a H 1 σ-finite and countably 1-rectifiable set (a countable union of compacts included in C 1hypersurfaces, up to a set of 0 H 1 -measure), ν(x) the common normal to all those hypersurfaces at a point x ∈ S(ϕ), ϕ ± (x) the values of ϕ(x) "above and below" S(ϕ), and C(ϕ) a measure (the Cantor part) which is mutually singular with dx and with H 1 (it only sees sets that have 0 Lebesguemeasure and infinite H 1 -measure). The subspace SBV (Ω) is that of those ϕ ∈ BV (Ω) such that C(ϕ) ≡ 0. It enjoys good compactness properties established in Ambrosio (1990), namely

ϕ n ∈ SBV (Ω) with        ϕ n bounded in L ∞ (Ω) ∇ϕ n bounded in L q (Ω; R 2 ), q > 1 H 1 (S(ϕ n )) bounded in R ∃{ϕ k(n) } ⊂ {ϕ n }, ∃ϕ ∈ SBV (Ω) s.t.        ϕ k(n) → ϕ, strongly in L p (Ω), p < ∞ ∇ϕ k(n) ⇀ ∇ϕ, weakly in L q (Ω; R 2 ) H 1 (S(ϕ)) ≤ lim inf n H 1 (S(ϕ k(n) ) (1.3.5)
Thanks to Ambrosio's compactness result, a simple argument of the direct method applied to (1.3.4) establishes existence of a minimizer ϕ g for that functional. The further result that the pair ϕ g , (S(ϕ g )) is a minimizer for (1.3.3) is highly non-trivial and makes up the bulk of De [START_REF] De Giorgi | Existence theorem for a minimum problem with free discontinuity set[END_REF].

In De Giorgi's footstep, we thus reformulate the variational evolution in the weak functional framework of SBV , or rather of those functions that have all their components in SBV , the jump set S(ϕ) becoming the union of the jump set of each component of ϕ. To do this, it is more convenient to view the hard device g(t) as living on all of R 2 and to integrate by parts the boundary term involving ġ(t) in (Eb). So, after elementary integrations by parts, we propose to investigate

• The weak variational evolution : Find (ϕ(t), Γ(t)) satisfying

(Ic) Γ(0) = Γ 0 ; (Ulm) (Γ(t), ϕ(t)
) is a local minimizer (in a topology that remains to be specified) for Remark that, in spite of the previous considerations on SBV , we have not explicitly indicated where ϕ (or ϕ(t)) should live. This is because, when dealing with vector-valued SBV -functions (the case of plane (hyper)elasticity, for example), that space -that is the Cartesian product of SBV for each of the component -is not quite sufficient. One should really work in GSBV Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF]. Forget about that here.

E(t; ϕ, Γ) := Ω W (∇ϕ) dx -F(t, ϕ) + kH 1 (Γ), (1.3.6) among all Ω\∂ s Ω ⊃ Γ ⊃ Γ(t)
Likewise, it is not so that the crack should belong to Ω \ ∂ s Ω. Any rigorous analysis will actually require ∂ s Ω, the site of application of the surface forces, to be part of the boundary of a non-brittle piece of the material. In other words, we should single out a thin layer around ∂ s Ω with infinite fracture toughness. This also will be overlooked in the sequel. Also the test cracks Γ do not have to be even essentially, i.e., up to a set of H 1 -measure 0, closed subsets of Ω \ ∂ s Ω, but only countably 1-rectifiable curves. Whether the actual crack Γ(t) that could be produced through the weak variational evolution is closed or not will be deemed a question of regularity and briefly commented upon in Paragraph 3.4 in the setting of global minimization.

Finally, as before, the same labels have been kept. The context will clearly indicate if the relevant formulation is weak or strong.

Remark 1.3.1 The crack Γ(t) can be identified with Γ 0 ∪ s≤t S(ϕ(s)) ; see Dal [START_REF] Maso | A variational model for quasistatic crack growth in nonlinear elasticity: some qualitative properties of the solutions[END_REF].

The recasting of Griffith's evolution model in a variational framework is now complete.

Stationarity versus local or global minimality -a comparison

We wish to explore the consequences of minimality. The adopted setting, or rather settings, for such an analysis are designed so that the "crack path" is not at stake. Nor is irreversibility a concern here because the monotonicity of the loads combined with the geometry of the problems result in an increase of both the measure of the discontinuity set and the magnitude of the discontinuities on that set. The focus is squarely on minimality, although, at times energy balance (Eb) will also be invoked.

The two settings are 1. A 1d-traction experiment under a hard or a soft device; 2. A 2d-tearing experiment.

In the first setting, cracks are merely points of discontinuity along the bar; in the second setting, symmetry of the geometry and of the loads suggests a straight crack path in mode III. In both settings, we assess the potential existence of weak variational evolutions satisfying unilateral stationarity (Ust), unilateral minimality (Ulm), or still unilateral global minimality (Ugm), together with energy balance (Eb).

1d traction

A "crack-free" homogeneous linearly elastic bar of length L, cross-sectional area Σ, Young's modulus E, toughness k is clamped at x = 0 and subject to a displacement load εL, 0 ≤ ε ր (hard device), or to a force load σΣ, 0 ≤ σ ր (soft device) at x = L. The parameters σ, ε play the role of the time variable. Thus, all evolutions will be parameterized by either σ, or ε. The results are concatenated in Conclusions 2.1.1, 2.1.4.

The soft device Assume that u is an admissible displacement field for a value σ of the loading parameter; that field may have jumps S(u) ⊂ [0, L], or it may correspond to the elastic state, in which case it lies in W 1,2 (0, L).

In any case we view it as a field defined in SBV (R) and such that u ≡ 0 on (-∞, 0). Its associated energy is

E(σ, u) = 1 2 (0,L) EΣ(u ′ ) 2 dx -σΣu(L+) + kΣ#(S(u)), (2.1.1)
and that energy will only be finite if S(u) is finite and u ′ ∈ L 2 (0, L), which we assume from now onward. This in turn implies that we may as well restrict the admissible fields to be in

SBV (R) ∩ L ∞ (R). Then, Conclusion 2.1.1
In a 1d traction experiment with a soft device, the elastic evolution is the only one that satisfies the weak variational evolution with either (Ust), or (Ulm), and (Eb). There is no solution to the weak variational evolution with (Ugm) and (Eb).

Remark 2.1.2 Testing the elastic solution against non-interpenetrating jumps is easy, since it suffices to restrict test jumps to be non-negative. In this context, the elastic solution is checked to be a global minimum for σ < 0, if non-interpenetration is imposed.

Remark 2.1.3 The above result demonstrates that soft devices prohibit global minimality. This is a significant drawback of global minimality and it clearly militates for a more local criterion. This has to be somewhat tempered, since one can build a reasonable class of non linear soft devices for which global minimality works; see Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF].

The hard device The admissible deformations are still in SBV (R) and they satisfy u ≡ 0 on (-∞, 0) and u ≡ εL on (L, ∞). The associated energy is

E(ε, u) = 1 2 (0,L) EΣ(u ′ ) 2 dx + kΣ#(S(u)), (2.1.2)
and, once again it is only finite if #(S(u)) is finite and u ′ ∈ L 2 (0, L), which we assume. Then, Conclusion 2.1.4 In a 1d traction experiment with a hard device, the elastic evolution, and all admissible evolutions with a set finite number of jumps satisfy the weak variational evolution with (Ulm) -and also (Ust) -and (Eb). Only u g (ε) defined as

u g (ε) = εx if 0 < ε ≤ 2k/EL 0, x ≤ a; εL, x > a if ε ≥ 2k/EL
satisfies the weak variational evolution with (Ugm) and (Eb). Also, all evolutions that are elastic, up to ε = 2ik/EL, then have i jumps satisfy (Ulm) -and also (Ust) -and (Eb).

A tearing experiment

Consider a thin semi-infinite homogeneous, linearly elastic slab of thickness 2H, Ω = (0, +∞) × (-H, +H). Its shear modulus is µ and its toughness k. Tearing amounts to a displacement load tHe 3 on {0} × (0, +H) and -tHe 3 on {0} × (-H, 0). The upper and lower edges are traction free and no forces are applied.

We assume throughout that all solutions respect geometric symmetry, although doing so cannot be justified; see in this respect the numerical experiment in Subsection 4.3. The symmetry assumption permits to look for anti-plane shear solution, anti-symmetric with respect to y = 0 and for a crack along that axis. We seek a displacement solution field of the form

u(x, y, t) = sign(y)u(t, x)e 3 with u(t, 0) = tH (2.2.1)
and note that such a displacement cannot be the exact solution, because it fails to ensure the continuity of the normal stress at the points (l(t), y), y = 0 (see (2.2.3)). The true symmetric solution can only be evaluated numerically, but it will be close to the proposed approximate solution as H becomes large.

The field u(t) will be discontinuous at the points x on the y = 0 -axis where u(x, t) = 0, that is S(u(t)) = {x ≥ 0 : u(t, x) = 0}. Then, the energy reads as

E(ϕ) = ∞ 0 µH(u ′ (x)) 2 dx + k ∞ 0 sg + (|u(x)|) dx, with sg + (z) := 0, z ≤ 0 1, z > 0 .
The kinematically admissible test fields u at time t are elements of W 1,2 (0, +∞) and satisfy u(t, 0) = tH. A global minimum for E exists for each t by elementary lower semi-continuity properties. We propose to show that (Ust), (Eb) has a unique solution, which identifies with the global minimum for E at t, which is thus unique.

Fix t. First, if u(t) is solution to (Ust), then it is such that, once it reaches 0, it stays equal to 0. Indeed, assume that a, with 0 ≤ a, is such that u(t, a) = 0. Take v with v = -u(t, x) in (a, ∞) and v = 0 otherwise. For h ∈ (0, 1), u + hv is an admissible test and

E(u + hv) -E(u) = (-2h + h 2 ) ∞ a u ′ (t, x) 2 dx +k ∞ a sg + (|(1 -h)u(t, x)|) -sg + (|u(t, x)|) dx = (-2h + h 2 ) ∞ a u ′ (t, x) 2 dx.
Thus, invoking (Ust),

0 ≤ d d h E(u + h v) h=0+ ≤ -2 ∞ a u ′ (t, x) 2 dx ≤ 0, so that u(t) = 0 in (a, ∞). But u(t) is continuous in x; thus, there exists ∞ ≥ l(t) > 0 such that S(u(t)) = [0, l(t)) with u(t, 0) = tH and u(t, l(t)) = 0.
We now perform an inner variation in E. Take v be in

C ∞ 0 (0, ∞). When |h| is sufficiently small, x → φ h (x) = x + hv(x) is a direct diffeomorphism onto R + . Moreover, if u(0) = tH, u h (0) = tH and u h converges to u when h goes to 0. The change of variables y = (φ h ) -1 (x) in the energy yields E(u(t) • φ -1 h ) = ∞ 0 µH u ′ (t, x) 2 φ ′ h (x) + φ ′ h (x)k sg + (u(t, x)) dx,
which in turn leads to

d dh E(u(t) • φ -1 h ) h=0 = ∞ 0 -µHu ′ (t, x) 2 + k sg + (u(t, x)) v ′ (x)dx.
Thus, provided that (Ust) holds for this kind of variations,

µHu ′ (t, x) 2 -k sg + (u(t, x)) = c, (2.2.2)
for some constant c. Now, take v ≥ 0 be in C ∞ 0 (0, l(t)). Then,

d dh E(u(t) + hv) h=0 = -2µH l(t) 0 u ′′ (t, x)v(x)dx.
Thus, invoking (Ust) again, we get that, on (0, l(t)), u ′′ ≤ 0, that is that u ′ is monotonically decreasing there. Now, if l(t) were infinite, we would have from (2.2.2) that u ′ ≡ d, some negative constant, which contradicts the convergence of u to 0 at infinity. Thus l(t) is finite and u ′ ≡ 0 on (l(t), ∞), hence c = 0, and

u ′ (t, x) = - k sg + (u(t, x)) µH , x > 0, u(0) = tH. (2.2.3)
We then conclude that the solution u(t) to (Ust) is unique and that it is given by

S(u) = [0, l(t)) with l(t) = tH µH k , while u(t, x) = tH 1 - x l(t) + .
Also note that l(t) and u(t) increase with t, so that irreversibility is automatic, while energy balance is guaranteed by the evoked smoothness of u(t).

Here, in contrast with the setting of Subsection 2.1, unilateral stationarity, unilateral local, or unilateral global minimality are indistinguishable, at least for an increasing load.

Remark 2.2.1 Note the surreptitious assumption that inner variations are valid tests for stationarity. In the presentation of Section 1, stationarity was introduced in the form of a combination of outer and inner variation (see (1.2.5)). It is in that sense that the re-formulated problem of Proposition 1.2.1 was equivalent to the original problem (1.1.1), (1.1.4) and an investigation of possible additional constraints on that problem resulting from the introduction of inner variations should be undertaken.

For a given length l of the tear (crack), the total energy at time t is immediately seen to be µH 3 t 2 /l + kl, hence strictly convex in l, so that, according to Remark 1.2.3, the smoothness of the evolution l(t) is hardly surprising.

Global minimality

In a Griffith setting, irreversibility is a simple notion: the crack can only extend with time, Γ(t) ⊃ Γ(s), s < t.

With that notion in mind, we now discuss the variational evolution in a global minimality setting, noting that existence in such a context will automatically provide existence of that evolution for any kind of local minimality criterion. Once again, the argument put forth in Paragraph 2.1 (see Remark 2.1.3) essentially prohibits force loads. We thus assume throughout this subsection that the only load is a displacement g(t) defined on ∂ d Ω, or rather, as we saw earlier in Subsection 1.3, on R 2 \ Ω.

Discrete evolution

As mentioned in the Introduction, the basic tool is also the natural computational tool: time discretization over the interval [0, T ]. We thus consider

t 0 = 0 < t n 1 < ...... < t n k(n) = T with k(n) n ր ∞, ∆ n := t n i+1 -t n i n ց 0.
Time-stepping the strong or weak minimality condition (Ugm), we obtain

(Sde) The strong discrete evolution: Find (Γ n i+1 , ϕ n i+1 ) a minimizer for min ϕ,Γ Ω\Γ W (∇ϕ) dx + kH 1 (Γ \ ∂ s Ω) : ϕ = g(t n i+1 ) on ∂ d Ω \ Γ; Γ ⊃ Γ n i ;
resp.

(Wde) The weak discrete evolution: Find ϕ n i+1 a minimizer for

min ϕ Ω W (∇ϕ) dx + kH 1 (S(ϕ) \ (Γ n i ∪ ∂ s Ω)) : ϕ = g(t n i+1 ) on ∂ d Ω \ S(ϕ) ; then, Γ n i+1 = Γ n i ∪ (S(ϕ n i+1 ) \ ∂ s Ω).
The balance (Eb) seems to have been forgotten all together in the discrete evolution, yet it will reappear in the time-continuous limit of those evolutions.

The first mathematical issue to confront is the existence of a solution to those discrete evolutions. As we mentioned before in Subsection 1.3, we cannot expect, even in 2d, a direct existence proof for the strong discrete evolution without imposing further restrictions on the class of admissible cracks. This is easily understood through the Neumann sieve example Murat (1985).

A Neumann sieve situation occurs when boundaries close up at a critical speed that creates channels of non-zero capacity in the domain. For example, consider Ω = (-1, 1) 2 and assume, in a linear anti-plane shear setting,

that the crack Γ n is given as {0} × [-1, 1] \ p=-n,...,n p n -e -n , p n +e -n n 0 n 1 n Figure 3.1.1. The Neumann sieve.
with

ϕ n = 0 on {x 1 = -1}, 1 on {x 1 = 1}.
Then ϕ n satisfies

-∆ϕ n = 0 on Ω n := (-1, 1) 2 \ Γ n ,
with ∂ϕ n ∂ν = 0 on all boundaries of Ω \ Γ n , except {x 1 = ∓1}. According to the results in Murat (1985) ϕ n → ϕ strongly in L 2 (Ω), with Ω = [(-1, 0) ∪ (0, 1)] × (-1, 1) and ϕ is the solution, for some µ = 0 of

-∆ϕ = 0 on Ω, with              ∂ϕ ∂x 2 = 0 on ∂Ω ∩ {x 2 = ±1} ϕ = 0, resp. 1 on {-1} × (-1, 1), resp. {1} × (-1, 1) ∂ϕ ∂x 1 = µ[ϕ] on {0} × (-1, 1).
Hence ϕ n does not converge to the solution φ = 0 on (-1, 0) × (-1, 1); 1 on (0, 1) × (-1, 1)

of the Neumann problem on Ω \ Γ, with Γ = {0} × (-1, 1). The Neumann sieve must thus be prevented so as to ensure the very existence of a pair-solution to the strong discrete evolution at each time step. A possible exit strategy consists in "prohibiting" disconnected cracks. A result of Chambolle and F. Doveri [START_REF] Chambolle | Continuity of Neumann linear elliptic problems on varying two-dimensional bounded open sets[END_REF] (see also [START_REF] Bucur | Boundary variation for a Neumann problem[END_REF]) states that, if Ω is a Lipschitz two dimensional domain and {Γ n } is a sequence of compact connected sets with H 1 (Γ n ) ≤ C and such that it converges -for the Hausdorff metric -to Γ, the solution to a Neumann problem of the form

     -∆ϕ n + ϕ n = g in Ω \ Γ n ∂ϕ n ∂ν = 0 on ∂[Ω \ Γ n ] is such that ϕ n , ∇ϕ n n -→ ϕ, ∇ϕ, strongly in L 2 (Ω), with ϕ solution to      -∆ϕ + ϕ = g in Ω \ Γ ∂ϕ ∂ν = 0 on ∂[Ω \ Γ]
An adaptation of that result by Dal Maso and Toader Dal [START_REF] Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF] proves the existence of a minimizer to the strong discrete evolution at each time step under the restriction that the cracks have an a priori set number of connected components. In turn, [START_REF] Chambolle | A density result in two-dimensional linearized elasticity, and applications[END_REF] proves an analogous result for plane elasticity. We will not discuss the strong evolution any further in these notes.

Note that the connectedness restriction can be weakened to include cracks with an a priori set number of connected components Dal [START_REF] Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF].

The discrete weak evolution behaves better as far as existence is concerned. Indeed, existence is a direct consequence of Ambrosio's compactness result (1.3.5), together with the following lower semi-continuity result which applies to the kind of elastic energy under consideration and to the sequence ϕ n in (1.3.5) (see Ambrosio (1994))

Ω W (∇ψ)dx ≤ lim inf n Ω W (∇ϕ n ) dx. (3.1.1)
(There is a slight modification in (1.3.5) which consists in replacing

H 1 by H 1 (Γ n i ∪ ∂ s Ω) c .
) To be precise, existence is established in • The anti-plane shear case: ϕ is scalar-valued and W is convex and has p-growth, p > 1; • The non-linear elasticity case: ϕ is vector-valued and W is quasiconvex with p-growth, p > 1. We refer the reader to the abundant literature on quasi-convexity (see e.g. [START_REF] Ball | W 1,p -quasiconvexity and variational problems for multiple integrals[END_REF]) for details on that notion; for our purpose, it suffices to remark that quasiconvexity, plus growth implies sequential weak lower semi-continuity on the Sobolev space W 1,p (Ω; R 2 ) [START_REF] Ball | W 1,p -quasiconvexity and variational problems for multiple integrals[END_REF], but also, see Ambrosio (1994), on

L ∞ (Ω; R 2 ) ∩ ϕ ∈ SBV (Ω; R 2 ) : ∇ϕ ∈ L p (Ω; R 2×2 ) .
It should be noted that the growth assumption prevents the energy density W (F ) from blowing up as det F ց 0, a desirable feature in hyperelasticity; and, most recently, • The case of finite elasticity: ϕ is vector-valued and W is poly-convex with blow-up as det F ց 0 and a weak form of non-interpenetration is imposed Dal Maso and Lazzaroni (2010). We will not consider this case in the sequel. Existence will not however be achieved in the setting of linearized elasticity which thus seems confined, for the time being, to the strong formulation.

Consider any setting for which the discrete evolution is meaningful. Then, for a given n (a given time step), we define the piecewise in time fields

       ϕ n (t) := ϕ n i Γ n (t) := Γ n i g n (t) = g(t n i ) on [t n i , t n i+1 ), and, for i = -1, Γ n -1 := Γ 0 .
Remark that irreversibility is guaranteed at the discrete level because of the definition of Γ n i in terms of its predecessors. In other words, Γ n (t) ր with t. Summing up, we have constructed, for each time t ∈ [0, T ], a pair (Γ n (t), ϕ n (t)) such that (Wde) The weak discrete evolution: ϕ n (t) is a minimizer for

min ϕ Ω W (∇ϕ) dx + kH 1 (S(ϕ)\(Γ n (t-∆ n ) ∪ ∂ s Ω)) : ϕ = g n (t) on ∂ d Ω\S(ϕ) and Γ n (t) = Γ n (t-∆ n ) ∪ S(ϕ n (t)).
Here again, the functional dependence of ϕ n (t) is not specified because it depends upon the scalar/vectorial nature of the specific problem, as well as on the coercivity/growth properties of the bulk energy density W .

At time t = 0, generically, it is not true that Γ n 0 (independent of n) ≡ Γ 0 , but merely that Γ n 0 ⊃ Γ 0 . There is an increase in the initial condition. The goal is to pass to the limit in n and hope that the limit fields will be solutions to the strong/weak variational evolutions. As will be seen below, this is not a straightforward proposition.

Global minimality in the limit

A usual first step in a limit process is to obtain n-independent a priori estimates on the fields. This will be obtained here upon testing the weak discrete evolution (Wde) at each time by appropriate test fields. The two choice test fields are g n (t) and ϕ n (t-∆ n )+g n (t)-g n (t-∆ n ) (the addition of the terms involving g n are so that the test deformations satisfy the boundary conditions at time t).

Then, provided we impose decent regularity on g, namely

g ∈ W 1,1 (0, T ; W 1,p (Ω(; R 2 )) ∩ L ∞ ((0, T ) × Ω(; R 2 )), (3.2.1)
for an energy with p > 1-growth, we obtain the following a priori bounds:

∇ϕ n (t) L p (Ω(;R 2 )) ≤ C H 1 (S(ϕ n (t)) ≤ C, (3.2.2) and H 1 (Γ n (t)) ≤ C, (3.2.3)
together with the following upper bound on the total energy

E n (t) := Ω W (∇ϕ n (t)) dx + kH 1 (Γ n (t)\ ∂ s Ω) ≤ E n (0) + τ n (t) 0 Ω ∂W ∂F (∇ϕ n (s)).∇ ġ(s) dx ds, (3.2.4)
where τ n (t) := sup{t n i ≤ t}. Remark that the derivation of (3.2.4) actually requires a bit of care; see Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF], Section 6.

It remains to pass to the n-limit in the minimality statement (Wde) under the above convergences and to avoid a Neumann sieve phenomenon as more and more crack components accumulate at a given time when n ր.

To this effect, remark that the circumstances that presided over the appearance of the Neumann sieve phenomenon in Paragraph 3.1 were deceiving, for they failed to account for the role played by the surface energy. Indeed, consider n large enough; the pair ϕ n , Γ n considered in that example cannot be a joint minimizer of

1 2 Ω\Γ |∇ϕ| 2 dx + H 1 (Γ), Γ ⊃ Γ n
with the same boundary conditions. By lower semi-continuity,

lim inf n 1 2 Ω\Γ |∇ϕ n | 2 dx + H 1 (Γ n ) ≥ 1 2 Ω\Γ |∇ϕ| 2 dx + 1,
with ϕ, the solution to the Neumann sieve. Now, ϕ has non zero bulk energy 1 2 Ω\Γ |∇ϕ| 2 dx, say C, so that, for n large enough,

1 2 Ω\Γ |∇ϕ n | 2 dx + H 1 (Γ n ) ≥ 1 + C 2 .
But the energy associated to the pair ({0}×[-1, 1], φ) is exactly 1, a strictly smaller value, while {0}×[-1, 1] ⊃ Γ n . For n large enough, closing the holes of the sieve and taking the crack to be {0}×[-1, 1] is the energetically sound choice. There is thus hope for a derivation of the global minimality condition (Ugm) in the weak variational evolution from (Wde) under refinement of the time step. That this is by no means a trivial endeavor can be illustrated as follows. We note first that, since H 1 (B \ A) ≥ H 1 (B) -H 1 (A), (Wde) implies in particular that ϕ n (t) is a minimizer for its own jump set, that is

1 2 Ω W (∇ϕ n (t)) dx ≤ 1 2 Ω W (∇ϕ) dx + kH 1 (S(ϕ) \ (S(ϕ n (t)) ∪ ∂Ω s )).
(3.2.5) If (Ugm) is to be obtained in the limit, then ϕ(t) should also in particular be a minimizer for its own jump set. In view of (3.2.2) and of the already quoted lower semi-continuity result of Ambrosio (1994), the left hand side of (3.2.5) is well behaved and the result would follow easily, provided that

lim sup n H 1 (S(ϕ) \ S(ϕ n (t))) ≤ H 1 (S(ϕ) \ S(ϕ(t))).
Consider however ϕ such that S(ϕ) ⊂ S(ϕ(t)), while the jump set of ϕ n (t) does not intersect that of ϕ(t) (which would surely happen if S(ϕ n (t)) ⊂ K n , with K n ∩K = ∅ and the Hausdorff distance from K n to K goes to 0). Then H 1 (S(ϕ)) must be 0!

The stability of the own jump set minimality condition cannot be established without a modification of the test fields ϕ. This is the essence of the jump transfer Theorem [START_REF] Francfort | Existence and convergence for quasi-static evolution in brittle fracture[END_REF], Section 2. We now quote it without proof in its simplest version.

Theorem 3.2.1 Let ϕ n , ϕ ∈ SBV (Ω) with H 1 (S(ϕ)) < ∞, be such that

• |∇ϕ n | weakly converges in L 1 (Ω);and

• ϕ n → ϕ in L 1 (Ω). Then, for every ζ ∈ SBV (Ω) with ∇ζ ∈ L p (Ω), 1 ≤ p < ∞, and H 1 (S(ζ)) < ∞, there exists a sequence {ζ n } ⊂ SBV (Ω) with ∇ζ n ∈ L p (Ω), such that

• ζ n → ζ strongly in L 1 (Ω); • ∇ζ n → ∇ζ strongly in L q (Ω); and • lim sup n H 1 A (S(ζ n ) \ S(ϕ n )) ≤ H 1 A (S(ζ) \ S(ϕ)) , for any Borel set A.
We fix a time t and recall (3.2.2). Ambrosio's compactness result permits to assert the existence of a t-dependent subsequence {ϕ nt (t)} of {ϕ n (t)} and of ϕ(t) such that the assumptions of Theorem 3.2.1 -or rather of a corollary of Theorem 3.2.1 which takes into account the boundary conditions on the test fields at t, namely ϕ nt (t) = g nt (t) on R 2 \ Ω -are met. The conclusion of that theorem then allows for a corresponding sequence {ζ nt } that is an admissible test in (W de), so that

Ω W (∇ϕ nt (t)) dx ≤ Ω W (∇ζ nt ) dx + kH 1 (S(ζ nt ) \ (S(ϕ nt (t)) ∪ ∂ s Ω)) ,
and then, from the convergences obtained in the theorem, together with the assumed p-growth of the energy, we pass to the limit in n t and obtain that the limit ϕ(t) is a minimizer for its own jump set, that is

Ω W (∇ϕ(t)) dx ≤ Ω W (∇ζ) dx + kH 1 (S(ζ) \ (S(ϕ(t)) ∪ ∂ s Ω)) .
We are inching ever closer to the global minimality statement (Ugm) in the weak variational evolution, but are not quite there yet, because we would like to remove not only S(ϕ(t))∪∂ s Ω but Γ(t)∪∂ s Ω in the minimality statement above. To do this, we need to define the limit crack Γ(t). There are various setting-dependent paths to a meaningful definition of the limit crack. An encompassing view of that issue is provided by the notion of σ p -convergence introduced in Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF], Section 4, a kind of set convergence for lower dimensional sets.

Definition 3.2.2 Γ n σ p -converges to Γ if H 1 (Γ n
) is bounded uniformly with respect to n, and

(1) Whenever

ϕ j , ϕ ∈ SBV (R 2 ) are such that            ϕ j weak- * ⇀ ϕ, in L ∞ (R) ∇ϕ j weak ⇀ ∇ϕ, in L p (R 2 ) S(ϕ j ) ⊂ Γ nj
for some sequence n j ր ∞, then S(ϕ) ⊂ Γ;

(2) there exist a function ϕ ∈ SBV p (R 2 ) with S(ϕ) = Γ and a sequence ϕ n with the properties of item (1).

The following compactness result proved in Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF], Section 4.2, holds true: Theorem 3.2.3 Let Γ n (t) be a sequence of increasing sets defined on [0, T ] and contained in a bounded set B. Assume that H 1 (Γ n (t)) is bounded uniformly with respect to n and t. Then there exist a subsequence Γ nj and Γ(t) defined on [0, T ] such that

Γ nj (t) σ p -converges to Γ(t), ∀t ∈ [0, T ].
The estimate (3.2.3) permits to apply the theorem above and thus to define a meaningful crack Γ(t) such that, for a subsequence still labeled Γ n (t), Γ n (t) σ p -converges to Γ(t), hence also Γ nt (t). Thanks to item (2) in Definition 3.2.2, we can construct ϕ with S(ϕ) = Γ(t) and ϕ nt satisfying the assumptions of Theorem 3.2.1 with S(ϕ nt ) ⊂ Γ nt (t). But (Wde) implies in particular that

Ω W (∇ϕ nt (t)) dx ≤ Ω W (∇ζ) dx + kH 1 (S(ζ)\(Γ nt (t) ∪ ∂ s Ω)) ≤ Ω W (∇ζ) dx + kH 1 (S(ζ)\(S(ϕ nt ) ∪ ∂ s Ω)).
and the jump transfer Theorem 3.2.1 delivers the minimality property (Ugm).

Having obtained global minimality, we still have to derive energy conservation (Eb). This is the object of the next paragraph.

Energy balance in the limit

Inequality (3.2.4) derived at the onset of Paragraph 3.2 hints at the possibility of an energy inequality. To obtain such an inequality in the limit, one should ensure that, as

n t ր ∞, Ω W (∇ϕ nt (t)) dx → Ω W (∇ϕ(t)) dx (3.
E(t) := Ω W (∇ϕ(t)) dx + kH 1 (Γ(t)\ ∂ s Ω) ≤ E(0) + t 0 Ω ∂W ∂F (∇ϕ(s)).∇ ġ(s) dx ds, (3.3.3)
The other energy inequality is a byproduct of the minimality statement (Ugm). Simply test global minimality at time s by ϕ(t) + g(s)g(t), t > s. Then, since S(ϕ(t) ⊂ Γ(t),

Ω W (∇ϕ(s)) dx ≤ Ω W (∇ϕ(t) + g(s) -g(t)) dx +H 1 (S(ϕ(t)) \ (Γ(s) ∪ ∂ s Ω)) ≤ Ω W (∇ϕ(t) + g(s) -g(t)) dx +H 1 (Γ(t) \ (Γ(s) ∪ ∂ s Ω)) = Ω W (∇ϕ(t)) dx + H 1 (Γ(t) \ (Γ(s) ∪ ∂ s Ω)) - Ω ∂W ∂F ∇ϕ(t) + ρ(s, t) t s ∇ ġ(τ ) dτ ) . t s ∇ ġ(τ ) dτ dx,
for some ρ(s, t) ∈ [0, 1]. Hence

E(t) -E(s) ≥ Ω ∂W ∂F ∇ϕ(t) + ρ(s, t) t s ∇ ġ(τ ) dτ ) . t s ∇ ġ(τ ) dτ dx.
We then choose a partition 0 < s n 1 < .... < s n k(n) = t of [0, t], with ∆ ′ n := s n i+1s n i ց 0; summing the contributions, we get

E(t) -E(0) ≥ k(n) i=0 Ω ∂W ∂F ∇ϕ(s n i+1 ) + ρ(s n i ,s n i+1 ) s n i+1 s n i ∇ ġ(τ ) dτ ) . s n i+1 s n i ∇ ġ(τ )dτ dx.
A uniform continuity type result -already implicitly used in the derivation of (3.2.4) -permits to drop the term depending on ρ(s n i ,s n i+1 ) in the previous inequality in the limit ∆ ′ n ց 0; see Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF], Section 4.3. Thus

E(t) -E(0) ≥ lim sup n    k(n) i=0 s n i+1 s n i Ω ∂W ∂F ∇ϕ(s n i+1 ) .∇ ġ(τ ) dx dτ    .
The expression on the right hand-side of the previous inequality looks very much like a Riemann sum. A not so well-known result in integration asserts that Riemann sums of a Lebesgue integrable function do converge to the integral of that function, but only for carefully chosen partitions [START_REF] Hahn | Über Annäherung an Lebesgue'sche integrale durch Riemann'sche summen[END_REF]. Since we enjoy complete liberty in our choice of the partition {s n j } of [0, t], we conclude that

E(t) -E(0) ≥ t 0 Ω ∂W ∂F (∇ϕ(s)) .∇ ġ(τ ) dx dτ,
which, together with (3.3.3), provides the required equality (Eb).

The time-continuous evolution

Here, the results obtained in the previous paragraphs are coalesced into an existence statement to the weak variational evolution. The result is expressed in a 2d setting, but it applies equally in a 3d setting, upon replacing H 1 by H 2 . We also recall similar existence results obtained in Dal [START_REF] Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF], [START_REF] Chambolle | A density result in two-dimensional linearized elasticity, and applications[END_REF] in the 2d connected case.

In what follows, the energy density W is a nonnegative convex -in the anti-plane shear setting -or quasiconvex -in the plane setting -

C 1 function on R 2 with (1/C)|F | p -C ≤ W (F ) ≤ C|F | p + C, ∀F, 1 < p < ∞.
Note that the assumptions on W immediately imply that (see, e.g., [START_REF] Dacorogna | Direct Methods in the Calculus of Variations[END_REF])

|DW (F )| ≤ C(1 + |F | p-1 ).
The domain Ω under consideration is assumed throughout to be Lipschitz and bounded , and the function g, which appears in the boundary condition on ∂Ω d , is assumed to be defined on all of R 2 ; actually, each of its components is taken to be in W 1,1 loc ([0, ∞); W 1,p (R 2 )). The traction-free part ∂ s Ω of the boundary is assumed to be closed. Finally, the pre-existing crack Γ 0 is a closed set in Ω, with H 1 (Γ 0 ) < ∞.

Theorem 3.4.1 ∃Γ(t) ⊂ Ω and ϕ such that

• each component of ϕ(t) ∈ SBV (R 2 ), with ∇ϕ p-integrable;

• Γ(t) ⊃ Γ 0 increases with t and H 1 (Γ(t)) < +∞;

• S(ϕ(t)) ⊂ Γ(t) ∪ ∂ s Ω and ϕ(t) = g(t) a.e. on R 2 \ Ω;

• For every t ≥ 0 the pair (ϕ(t), Γ(t)) minimizes

Ω W (∇ϕ) dx + kH 1 (Γ \ ∂ s Ω)
among all Γ ⊃ Γ(t) and ϕ with components in SBV (R 2 ) s.t. S(ϕ) ⊂ Γ and ϕ = g(t) a.e. on R 2 \ Ω; • the total energy

E(t) := Ω W (∇ϕ(t)) dx + kH 1 (Γ(t) \ ∂ s Ω) is absolutely continuous, DW (∇ϕ) • ∇ ġ ∈ L 1 loc ([0, ∞); L 1 (R 2 )), and
E(t) = E(0) + t 0 Ω DW (∇ϕ(s)) • ∇ ġ(s) dx ds.
We could incorporate body or surface loads, provided they belong to a certain class of soft devices Dal [START_REF] Maso | Quasistatic crack growth in nonlinear elasticity[END_REF], Section 3.

Note that, in the vector-valued setting, it is assumed that somehow, the deformations are always capped in sup-norm by some set number. This is an a-priori assumption which can be verified for certain classes of quasi-convex energies [START_REF] Leonetti | Maximum principle for vector valued minimizers[END_REF]. There is no need for such an assumption in the anti-plane shear case, provided that the displacement load g is also bounded in sup-norm.

In 2d only and in the case where the cracks are assumed a priori to be connected -or to have a pre-set number of connected components -then the same existence result for the strong variational evolution is obtained in Dal [START_REF] Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF] in the quadratic case and in [START_REF] Chambolle | A density result in two-dimensional linearized elasticity, and applications[END_REF] in the case of linearized elasticity. The statement is identical to that of Theorem 3.4.1 at the expense of replacing Ω by Ω\Γ , and considering ϕ's with components in

L 1,2 (Ω \ Γ) := {f ∈ L 2 loc (Ω \ Γ) : ∇f ∈ L 2 (Ω \ Γ)}, resp. ϕ ∈ LD(Ω \ Γ) := {ζ ∈ L 2 loc (Ω \ Γ; R 2 ) : e(ζ) ∈ L 2 (Ω \ Γ; R 4 )}.
in the case of linear elasticity.

This existence result calls for comments. First and foremost, it is an existence result, not a uniqueness result. As in other non-convex problems in mechanics, uniqueness should not be expected.

Then the lack of regularity of the field ϕ(t) indicates that time jumps could appear in the various fields. Indeed, still referring to that same example, we witness there a brutal decrease to 0 at time t i of the bulk energy with a corresponding increase of the surface energy. The total energy (Eb) will remain impervious to those jumps.

Third, an implicit change of initial conditions may occur, since it might happen that Γ(0) contains, but does not equal Γ 0 .

Finally, the weak evolution might just turn out to be a strong evolution, as was the case for image segmentation thanks to De [START_REF] De Giorgi | Existence theorem for a minimum problem with free discontinuity set[END_REF], in which case there would be no need for the strong variational evolution. Recent results of J. F. Babadjian and A. Giacomini seem to confirm this in a 2d setting.

Remark 3.4.2 The unilateral global minimality condition (item 4. in Theorem 3.4.1) can actually be strengthened as follows: For every t ≥ 0 the pair (ϕ(t), Γ(t)) minimizes

Ω W (∇ϕ) dx + kH 1 (Γ \ ∂ s Ω) among all Γ ⊃ ∪ s<t Γ(s) and ϕ with components in SBV (R 2 ) s.t. S(ϕ) ⊂ Γ and ϕ = g(t) a.e. on R 2 \ Ω.
This states that the admissible test cracks do not have to contain the current crack, but only those up to, but not including the current time, a clearly stronger minimality condition. The two conditions are actually equivalent because, for s < t, unilateral global minimality implies in particular that

Ω W (∇ϕ(s)) dx + kH 1 (Γ(s) \ ∂ s Ω) ≤ Ω W (∇ϕ + ∇g(s) -∇g(t)) dx +kH 1 (Γ \ ∂ s Ω),
for any ϕ with components in SBV (R 2 ) s.t. ϕ = g(t) a.e. on R 2 \Ω, S(ϕ) ⊂ Γ and any Γ ⊃ ∪ s<t Γ(s). Let s ր t and use item 5. (the continuity of the total energy) to pass to the limit in the left hand-side of the inequality above. The stronger minimality result is then obtained by dominated convergence (since W has p-growth).

Numerics

At first glance, numerical implementation of the variational approach advocated in these notes is hopeless because the "classical" numerical methods dealing with discontinuous displacement fields rely on some non-negligible amount of a priori knowledge of that path. This includes the extended finite element method and other enrichment-based variants. A proper discretization scheme for the total energy needs to both approximate potentially discontinuous displacement fields -and thus the position of their discontinuity sets -and to lead to an accurate and isotropic approximation of the surface energy. Such a scheme does not easily accommodate cohesive finite element methods or discontinuous Galerkin methods. Note that this is partially addressed by a careful estimate of the anisotropy induced by the mesh in [START_REF] Negri | The anisotropy introduced by the mesh in the finite element approximation of the Mumford-Shah functional[END_REF], [START_REF] Negri | A finite element approximation of the Griffith model in fracture mechanics[END_REF] or still through the use of adaptive finite element methods [START_REF] Bourdin | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF].

Further, if the variational framework contends that it addresses crack initiation and crack propagation in a unified framework, the same should be true of the numerical method. In particular, methods based on considering energy restitution caused by small increments of existing cracks are ruled out. It can actually be shown (see [START_REF] Chambolle | Crack initiation in brittle materials[END_REF]) that "small" cracks will never lead to descent directions for the global minimization of the total energy in the absence of strong singularities in the elastic field.

Non-convexity of the total energy is yet another major obstacle to overcome. The typical size of the discrete problems prohibits appeal to global or non-deterministic optimization techniques. Global minimization of the energy is an arguable postulate, but it is at present the only one theoretically suitable for a thorough investigation of any numerical implementation. So far, the only exception to that would be the 1d setting where one can establish convergence of the critical points of the approximating functional discussed later on in this section [START_REF] Francfort | Critical points of ambrosio-tortorelli converge to critical points of mumford-shah in the one-dimensional dirichlet case[END_REF].

Since, as already noted at the onset of Section 3, global minimality deals badly with force loads, the only loads considered throughout this section are displacement loads.

The numerical method that will be described below finds, once again, its inspiration in the Mumford-Shah functional for image segmentation (see Subsection 1.3). The main ingredients were first introduced in the latter context in Ambrosio and Tortorelli (1990), [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], [START_REF] Bellettini | Discrete approximation of a free discontinuity problem[END_REF], [START_REF] Bourdin | Une méthode variationnelle en mécanique de la rupture[END_REF], [START_REF] Bourdin | Image segmentation with a finite element method[END_REF], [START_REF] Negri | Numerical minimization of the Mumford-Shah functional[END_REF] and later adapted to fracture in Bourdin et al. (2000), [START_REF] Giacomini | A discontinuous finite element approximation of quasi-static growth of brittle fractures[END_REF], Chambolle (2004), [START_REF] Chambolle | Addendum to "an approximation result for special functions with bounded deformation[END_REF], [START_REF] Giacomini | Ambrosio-Tortorelli approximation of quasi-static evolution of brittle fractures[END_REF], [START_REF] Giacomini | Discontinuous finite element approximation of quasistatic crack growth in nonlinear elasticity[END_REF].

The method allows for an isotropic and mesh independent approximation of the total energy. It copes rather successfully with both initiation and propagation as seen through the various numerical experiments presented in Subsection 4.3. Like the actual variational model, it applies to the one, two, or three dimensional cases without alteration.

Finally, time dependence will be approached through time discretization, and all computations will be performed for a sequence of times

t 0 = 0 < t n 1 < ...... < t n k(n) = T with k(n) n ր ∞, ∆ n := t n i+1 -t n i n
ց 0. We will mostly drop the n-dependence, unless explicitly referring to the putative convergence of the time-discrete evolution to the time-continuous evolution.

Numerical approximation of the energy

The essence of the numerical implementation is to be found in the concept of variational convergence. Specifically, the first step consists in devising a good approximation of the total energy in the sense of Γ-convergence. We refer the reader to Dal [START_REF] Maso | An introduction to Γ-convergence[END_REF]), Braides (2002) for a complete exposition of the underlying theory.

Consider a R-valued functional F defined over, say a metrizable topological space X, and a sequence F ε of the same type. Then, F ε Γ-converges to F as ε ց 0 iff the following two conditions are satisfied for any u ∈ X:

1. lower-inequality: for any sequence (u ε ) ε ∈ X converging to u,

lim inf ε→0 F ε (u ε ) ≥ F(u); (4.1.1)
2. existence of a recovery-sequence: there exists a sequence (u

ε ) ε ∈ X converging to u, such that lim sup ε→0 F ε (u ε ) ≤ F(u). (4.1.2)
The interest of Γ-convergence from the standpoint of numerics lies in the following elementary theorem in Γ-convergence:

Theorem 4.1.1 If F ε Γ-converges to F and u * ε is a minimizer of F ε and if, further, the sequence u * ε is compact in X, then there exists u * ∈ X such that u * ε → u, u * is a global minimizer for F, and

F ε (u * ε ) → F(u * ).
Stability of global minimizers under Γ-convergence is indeed a powerful numerical tool. Rather than attempting to minimize the total energythus having to reconcile discretization and discontinuous functions -we propose to construct, at each time step t i , a family of regularized energies E i ε that Γ-converge to E i , the energy for the weak variational evolution at that time step (see (1.3.6)). In the footstep of Ambrosio and Tortorelli (1990), [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], we will approximate the potentially discontinuous field ϕ i and its crack set Γ i by two smooth functions. The implementation of the first time step, which is very close to that of the original approximation in the context of the Mumford-Shah functional, is presented in Paragraph 4.1. while Paragraph 4.1 shows how to account for irreversibility and approximate the weak discrete time evolution (Wde).

The first time step Consider the first time step of the weak discrete evolution under the unilateral global minimality condition (Ugm). The irreversibility condition is trivially satisfied, so that it suffices to minimize the total energy

E(ϕ) = Ω W (∇ϕ)dx + kH 1 (S(ϕ))
with respect to any kinematically admissible ϕ. In all that follows, Ω denotes a "large enough" open bounded set such that Ω ⊂ Ω, and the Dirichlet boundary conditions are enforced on Ω \ Ω, not on R 2 \ Ω because, as will be seen below, the computations are performed on that larger domain, and not only on Ω.

Following Ambrosio and Tortorelli (1990), [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], we introduce a secondary variable v ∈ W 1,2 ( Ω \ ∂ s Ω) and two small positive parameters ε, and η ε = o(ε), and define, for any kinematically admissible ϕ,

F(ϕ, v) =      Ω W (∇ϕ) dx + kH N -1 (S(ϕ) \ ∂ S Ω) if v = 1 a.e. +∞ otherwise, (4.1.3)
and

F ε (ϕ, v) = Ω (v 2 + η ε )W (∇ϕ) dx + k Ω\∂sΩ (1 -v) 2 4ε + ε|∇v| 2 dx.
(4.1.4) In the anti-plane case, proving the Γ-convergence of F ε to F is a simple adaptation of Ambrosio and Tortorelli's result (see [START_REF] Bourdin | Une méthode variationnelle en mécanique de la rupture[END_REF]) while it is more involved in that of linearized elasticity Chambolle (2004). We limit the analysis to the former case assuming that the energy is quadratic in the field, i.e., W (F ) := 1/2µ|F | 2 . The proof of the lower inequality of Theorem 4.1.1 is technical and does not shed much light on the proposed numerical method. Instead, we present a simpler but weaker inequality i.e. a version of (4.1.1) with a "wrong" constant that still highlights the link between the regularized and weak energies. The construction of the recovery sequence in Theorem 4.1.1 provides valuable insight and we propose to detail it, at least when the target is a mildly regular kinematically admissible field for E. Actually, deriving the lim-sup inequality (4.1.2) for minimizers can easily be seen to be no restriction. But for those, the mild regularity assumption below holds true, at least in anti-plane shear and energy densities of the form |F | p with p > 1.

We start with the compactness of minimizing sequences. Let (ϕ ε , v ε ) be a sequence of minimizers for F ε . We show that v ε → 1 almost everywhere and that there exists ϕ ∈ SBV (Ω) such that ϕ ε → φ in L 2 .

A simple truncation argument shows that we are at liberty to assume that

ϕ ǫ L ∞ ≤ C (4.1.5) and 0 ≤ v ǫ ≤ 1. (4.1.6)
Using the classical inequality a 2 + b 2 ≥ 2ab, we also have that

(1 -v ε ) 2 4ε + ε 2 |∇v ε | 2 ≥ (1 -v ε )|∇v ε |. (4.1.7)
Finally, testing with v = 1 and ϕ kinematically admissible, we get

F ε (ϕ ε , v ε ) ≤ C. (4.1.8) From (4.1.8), v ε → 1 a.e.
in Ω. In order to obtain the compactness of the sequence ϕ ε , we consider the function

ω ε := (2v ε -v 2 ε )ϕ ε and note that ω ε is uniformly bounded in L ∞ (Ω). We have that ∇ω ε = (2v ε -v ε ) 2 ∇ϕ ε + 2(1 -v ε )∇v ε ϕ ε .
From (4.1.6) and (4.1.8), we easily obtain that (2v εv ε ) 2 ∇ϕ ε is uniformly bounded in L 1 and from (4.1.7), (4.1.5), and (4.1.8), that 2(1v ε )∇v ε ϕ ε is also bounded in L 1 . Thus ω ε is uniformly bounded in L ∞ (Ω) ∩ BV (Ω), so that there exists ω such that ω ǫ → ω a.e. in Ω.

Finally, remark that

ϕ ε = ω ε /2v ε -v 2 ε and that, since ω ε → ω and v ε → 1 a.e. in Ω, ϕ := ω is such that ϕ ε → ϕ in L 2 (Ω).
We now focus on the lower inequality of Γ-convergence, or rather on a weaker inequality. Take any

(ϕ ε , v ε ) such that (ϕ ε , v ε ) → (ϕ, 1) in L 2 (Ω) 2 .
Then,

lim inf ε→0 F ε (ϕ ε , v ε ) ≥ Ω W (∇ϕ) dx + k 2 H N -1 (S(ϕ)). (4.1.9)
We may as well assume that (4.1.5), (4.1.6), (4.1.8) hold true, in which case ϕ ∈ BV (Ω).

Using (4.1.7) again,

F ε (ϕ ε , v ε ) ≥ Ω v 2 ε |∇ϕ ε | 2 dx + k Ω\∂sΩ (1 -v ε )|∇v ε | dx.
By increasing rearrangement,

Ω v 2 ε |∇ϕ ε | 2 dx = 1 0 2s {vε>s} |∇ϕ ε | 2 dx ds,
and using the co-area formula for BV -functions (see for instance [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]; [START_REF] Evans | Measure theory and fine properties of functions[END_REF]),

Ω\∂sΩ (1 -v ε )|∇v ε | dx = 1 0 (1 -s)H N -1 (∂ * {v ε > s} \ ∂ s Ω) ds.
Putting the above expressions together, we obtain

F ε (ϕ ε , v ε ) ≥ 1 0 2s {vε>s} |∇ϕ ε | 2 dx ds + k 1 0 (1 -s)H N -1 (∂ * {v ε > s} \ ∂ s Ω) ds. (4.1.10) Consider any s = 1 and define ω s ε := ϕ ε χ {vε>s} . Since ϕ ε ∈ W 1,2 and {v ε > s} is a set of finite perimeter, ω s ε ∈ SBV , ∂ * {v ε > s} = S(ω s ε ),
the jump set (set of non Lebesgue points) of ω s ε , and

{vε>s} |∇ϕ ε | 2 dx = Ω |∇ω s ε | 2 dx.
Since v ε → 1 a.e., we also obtain that ω s ε → ϕ in L 2 . Combining (4.1.10) and the two identities above and using Fatou's lemma, we obtain

lim inf ε→0 F ε (ϕ ε , v ε ) ≥ 1 0 2s lim inf ε→0 Ω |∇ω s ε | 2 dx ds + k 1 0 (1 -s) lim inf ε→0 H N -1 (S(ω s ε ) \ ∂ s Ω) ds.
By Ambrosio's Compactness and lower semi-continuity theorems in SBV , we obtain

lim inf ε→0 F ε (ϕ ε , v ε ) ≥ 1 0 2s Ω |∇ϕ| 2 dx ds + k 1 0 (1 -s)H N -1 (S(ϕ) \ ∂ s Ω) ds.
Integrating in s yields (4.1.9).

We finally turn to the construction of the recovery sequence (4.1.2). The following construction does not account for the Dirichlet boundary condition and the interested reader is referred to [START_REF] Bourdin | Une méthode variationnelle en mécanique de la rupture[END_REF] for the corresponding technicalities. As a corollary, we may as well take Ω ≡ Ω in the construction of the attainment sequence that follows. In the case of interest to us, i.e., that with Dirichlet boundary conditions on a part ∂ d Ω = ∂Ω \ ∂ s Ω it will be enough to re-introduce Ω \ ∂Ω s in lieu of Ω in the second integral in (4.1.4).

We also assume that ϕ is a solution to the minimization of E that satisfies

H 1 (S(ϕ)) = H 1 (S(ϕ)). (4.1.11)
For minimizers of the Mumford-Shah functional, this mild, albeit difficult regularity property was established in De [START_REF] De Giorgi | Existence theorem for a minimum problem with free discontinuity set[END_REF]. In the scalarvalued setting, the case of a certain class of convex bulk energies which includes p > 1-homogeneous energies was investigated in [START_REF] Fonseca | Regularity results for anisotropic image segmentation models[END_REF]. The regularity result was generalized to our setting in [START_REF] Bourdin | Une méthode variationnelle en mécanique de la rupture[END_REF], at least for minimizers in anti-plane shear with a quadratic elastic energy density. The closure property (4.1.11) is not so clearly true in more general settings, and different approximation processes must be used in such cases; the interested reader is invited to consult e.g. Braides (2002). Consider a kinematically admissible field ϕ -an element of SBV (Ω)satisfying (4.1.11). Define

d(x) := dist(x, S(ϕ)).
The volume of the area bounded by the s-level set of d is

ℓ(s) := x ∈ R 2 ; d(x) ≤ s .
The distance function is 1-Lipschitz, i.e., |∇d(x)| = 1 a.e., while, by the co-area formula for Lipschitz functions (see e.g. [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]),

ℓ(s) = s 0 H 1 ({x ; d(x) = t}) dt, so that, in particular, ℓ ′ (s) = H 1 ({x ; d(x) = s}) .
(4.1.12) Also, see [START_REF] Federer | Geometric measure theory[END_REF]-3.2.39, lim s→0 ℓ(s) 2s = H 1 (S(ϕ)).

We choose α ε such that

α ε = o(ε), η ε = o(α ε ), which is possible since η ε = o(ε)
, and define the functions (4.1.13) and

v ε (x) :=        0 if d(x) ≤ α ε 1 -exp - d(x) -α ε 2ε otherwise,
ϕ ε (x) :=      d(x) α ε ϕ(x) if 0 ≤ d(x) ≤ α ε ϕ(x) otherwise.
Note that it is easily seen that

ϕ ε ∈ W 1,2 (Ω). Further, ϕ ε → ϕ in L 2 (Ω), and v ε → 1 almost everywhere. Since v ε ≤ 1, Ω v 2 ε + η ε |∇ϕ ε | 2 dx ≤ d(x)≤αε η ε |∇ϕ ε | 2 dx + d(x)≥αε (1 + η ε )|∇ϕ| 2 dx.
Observe now that, for d(x) ≤ α ε , ∇ϕ ε = d(x)/α ε ∇ϕ + 1/α ε ϕ∇d, so, in view of the 1-Lipschitz character of d and of the L ∞ -bound on ϕ,

Ω (v 2 ε + η ε )|∇ϕ ε | 2 dx ≤ 2 η ε d(x)≤αε |∇ϕ| 2 dx + M 2 η ε α 2 ε ℓ(α ε ) + d(x)≥αε (1 + η ε )|∇ϕ| 2 dx.
Since Ω |∇ϕ| 2 dx < ∞, the first term in the parenthesis on the right hand side above converges to 0 as ε → 0. Recalling that ℓ(α ε )/α ε = O(1), while η ε /α ε = o(1) permits one to conclude that the limit of the second term in that parenthesis also converges to 0 with ε. We conclude that lim sup

ε→0 Ω v 2 ε + η ε |∇ϕ ε | 2 dx ≤ Ω |∇ϕ| 2 dx. (4.1.14)
Let us examine the surface energy term. Using once again the 1-Lipschitz character of d, together with the co-area formula, we get

Ω ε|∇v ε | 2 + (1 -v ε ) 2 4ε dx ≤ ℓ(α ε ) 4ε + d(x)≥αε 1 2ε exp(- d(x) -α ε ε ) ≤ ℓ(α ε ) 4ε + 1 2ε ∞ αε exp(- s -α ε ε )H 1 ({d(x) =
ε|∇v ε | 2 + (1 -v ε ) 2 4ε dx ≤ H 1 (S(ϕ)).
(4.1.17)

Collecting (4.1.14), (4.1.17) gives the upper Γ-limit inequality.

Remark 4.1.2 The form of the field v ε in (4.1.13) may seem somewhat ad-hoc. It is not. The choice of the profile for the field v ε is derived from the solution of an "optimal profile" problem (see Alberti (2000)). Consider, in e.g. 2d, a point x on the crack and a line orthogonal to the crack and passing through x, parameterized by the variable s. Consider the restriction of the regularized surface energy to this line

F ε,x (s) = k ∞ 0 (1 -v(s)) 2 4ε + ε|v ′ (s)| 2 ds.
Then the profile

v ε (s) = 1 -exp - (s -α ε ) 2ε
corresponds to the minimizer of F ε,x under the following boundary conditions:

v ε (α ε ) = 0; lim s→∞ v ε (s) = 1.
Indeed, it is also possible to construct the field v ε for the upper Γ-limit along lines intersecting the crack set at 90 0 angles, using the solution to the optimal profile problem on each of those. Integration of the result along the crack set will also permit one to recover the upper Γ-limit.

The Γ-convergence result above can be extended to the restriction F ε,h of F ε to a linear finite element approximation, provided that the discretization parameter h is such that h = o(ε) (see [START_REF] Bellettini | Discrete approximation of a free discontinuity problem[END_REF], [START_REF] Bourdin | Image segmentation with a finite element method[END_REF]). A closer look at the construction for the upper Γ-limit and at its adaptation to F ε,h provides some useful insight into possible error estimates.

The construction of the sequence (ϕ ε,h , v ε,h ) for the upper Γ-limit for F ε,h can be obtained from that above. Let T h be a conforming mesh of Ω \ ∂Ω s and S h be the set of all elements in T h intersecting S(ϕ). Let π h be a linear finite element projection operator associated with T h , and consider

v ε,h (x) := 0 if x ∈ S h ; π h (v ε ) otherwise, (4.1.18) and ϕ ε,h (x) := π h (ϕ ε ) . (4.1.19)
Following a path similar that developed in the computation of the upper Γ-limit above, the first term ℓ(α ε )/4ε on the right hand-side of inequality (4.1.15) becomes |S h |/4ε ≃ H 1 (S(ϕ))h/4ε, which converges to 0 only if h = o(ε). The consideration of quadratic finite elements in lieu of linear ones would still induce an error on the surface energy of the order of h/ε, albeit with a different constant. This is why the proposed implementation only resorts to piecewise linear finite elements for ϕ and v.

In a different direction, this term links the anisotropy of the mesh to the quality of the approximation of the surface energy. In [START_REF] Negri | The anisotropy introduced by the mesh in the finite element approximation of the Mumford-Shah functional[END_REF], M. Negri studied the effect of various types of structured meshes on the surface energy for the Mumford-Shah problem. In the numerical experiments, the isotropy of the surface term is ensured through the use of "almost" isotropic Delaunay meshes.

From the construction above, it is deduced that the relation h = o(ε) only needs to be satisfied "close" to S(ϕ). Of course, barring prior knowledge of S(ϕ), uniformly homogeneous fine meshes are a must. However, a posteriori re-meshing the domain will then improve the accuracy of the energy estimate. However, a priori mesh adaption is not advisable because the local size of the mesh affects the quality of the approximation of the surface energy and can potentially create spurious local minimizers. So, a posteriori mesh refinement around the cracks shields the computations from artificial cracks that would correspond to local minima created by a priori mesh refinement! Note that some recent computations in [START_REF] Burke | An adaptive finite element approximation of a variational model of brittle fracture[END_REF] seem to deal rather successfully with a priori mesh adaption.

The sequence for the upper Γ-limit is also admissible for the lower Γlimit, so that, if ϕ if a minimizer for the total energy, the sequence (ϕ ε , v ε ) constructed above approximates a minimizing sequence for F ε,h and this asymptotically in h, that is

Ω\∂Ωs (1 -v ε ) 2 4ε + ε|∇v ε | 2 dx = 1 + h 4ε H 1 (S(ϕ)). (4.1.20)
In practice, it is as if the fracture toughness had been amplified by a factor 1 + h/4ε, which has to be accounted for when interpreting the results.

The experiments in Section 4.3 highlight the effect of mesh isotropy on the results, and show how the fracture toughness is overestimated.

Quasi-static evolution The approximation scheme devised in Subsection 4.1 should now be reconciled with the evolutionary character of the weak discrete formulation. Irreversibility of the crack growth is enforced at the time-discrete level in the manner described below. Consider a fixed ε and a fixed conforming mesh T h of Ω \ ∂Ω s with characteristic element size h. Introduce a small parameter η > 0, and at each step t i , the set of vertices

K i ε,h,η := s ∈ T h ; v i ε,h (s) ≤ η , i > 0; K 0 ε,h,η := ∅.
In the light of the Γ-convergence properties of F ε,h , the crack growth condition translates into a growth condition on the sets K i ε,h,η and leads to the following fully spatially and temporally discrete evolution scheme: Recently, Giacomini conducted a rigorous analysis of a slightly different approach to the time evolution for F ε . In [START_REF] Giacomini | Ambrosio-Tortorelli approximation of quasi-static evolution of brittle fractures[END_REF], crack growth is enforced through the monotonicity of v in time, i.e., by successively minimizing F ε among all (ϕ, v) such that ϕ = g(t i+1 ) on Ω \ Ω, and v ≤ v ε i almost everywhere on Ω. In that setting, as both the time discretization parameter (∆ n ) and ε go to 0 (in a carefully ordered fashion), the discrete evolution converges to a continuous evolution satisfying the conclusions of Theorem 3.4.1.

In the forthcoming numerical experiments, crack growth is enforced as described in (4.1.21). Implementing monotonicity would not generate additional difficulties, but only slightly increase the computational cost as equality constraints would have to be replaced by inequality ("box") constraints.

Remark 4.1.3 The Γ-convergence based approach to minimization is not so easily amenable to the treatment of local minimization. If (ϕ, 1) is an isolated L 1 -local minimizer for F (see (4.1.3)), then Theorem 2.1 in [START_REF] Kohn | Local minimizers and singular perturbations[END_REF] can be adapted to the current setting to prove the existence of a sequence of L 1 -local minimizers (ϕ ε , v ε ) for F ε converging to (ϕ, 1) in L 1 . Unfortunately, the isolation hypothesis is generically false: see for instance the 1d-traction experiment with a hard device in Paragraph 2.1.

Even when the isolation hypothesis applies, the above-mentioned theorem grants the existence of a sequence of local minimizers for F ε converging to a local minimizer of F, but does not however guarantee that a converging sequence of local minimizers for F ε converges to a local minimizer for F.

The only positive results in this direction concern the numerically uninteresting 1d case [START_REF] Francfort | Critical points of ambrosio-tortorelli converge to critical points of mumford-shah in the one-dimensional dirichlet case[END_REF].

Minimization algorithm

Recall that F ε,h is the restriction of F ε defined in (4.1.4) to a linear finite element approximation. Also note that, although F ε is separately convex in its arguments ϕ and v, it is not convex in the pair (ϕ, v).

In the numerical experiments below, we fix the regularization parameter ε and generate a mesh with characteristic size h. We do not try to adapt the values on ε and h during the numerical minimization of F ε,h . Thus, the numerical implementation reduces to a sequence of minimizations for F ε,h , each corresponding to a separate time step. All presented experiments have been tested on meshes of various size and with different values of the parameter ε and/or of the time discretization length; the results seem impervious to such changes, at least for reasonably small choices of the parameters ε, h, ∆ n .

Because of the lack of convexity of F ε,h , the minimization scheme is purely heuristic. As per Section 4.1, we should choose a mesh size h which remains "small" compared to the regularization parameter, which in turn needs to be "small". In a 2d setting, this typically results in meshes with (10) 5 elements, while in 3d, meshes will consist of over a million elements.

Although the analysis of such large problems can be tackled thanks to the wider availability of massively parallel computers, there are, to our knowledge, no global minimization algorithms capable of handling them. At best, the algorithms will satisfy necessary optimality conditions for minimality.

The alternate minimization algorithm The first building block in the numerical implementation is an alternate minimization algorithm, leading to evolutions satisfying a first set of necessary conditions for optimality.

The functional F ε -and therefore F ε,h -is Gateaux-differentiable around any (ϕ, v). We compute the first order variation of F ε,h around any kinematically admissible (ϕ, v) in the directions ( φ, 0) and (0, ṽ), where φ and ṽ are admissible variations ( φ = 0 on Ω \ Ω and ṽ = 0 on K i ε,h,η ) and obtain that the solution (ϕ i+1 ε,h , v i+1 ε,h ) of the fully discrete evolution at time step

t i+1 satisfies                    Ω (v i+1 ε,h ) 2 + η ε DW (∇ϕ i+1 ε,h ).∇ φ dx = 0 Ω v i+1 ε,h ṽ W (∇ϕ i+1 ε,h ) dx + k Ω\∂Ωs v i+1 ε,h ṽ 4ε + ε∇v i+1 ε,h • ∇ṽ dx = k Ω\∂Ωs ṽ 4ε dx.
(4.2.1) This leads to the following algorithm (δ is a fixed tolerance parameter):

Algorithm 1 The alternate minimization algorithm:

1: let p = 0 and v (0) := v i ε,h . 2: repeat 3: p ← p + 1 4: compute ϕ (p) := arg min ϕ F ε,h (ϕ, v (p-1) ) under the constraint ϕ (p) = g(t i+1 ) on Ω \ Ω. 5: compute v (p) := arg min v F ε,h (ϕ (p) , v) under the constraint v (p) = 0 on K i ε,h,η 6: until v (p) -v (p-1) ∞ ≤ δ 7: set ϕ i+1 ε,h := ϕ (p) and v i+1 ε,h := v (p)
Since F ε,h is separately convex in each of its arguments, the algorithm constructs at each time step a sequence with decreasing total energy; it is therefore unconditionally convergent in energy. A more detailed analysis conducted in [START_REF] Bourdin | Numerical implementation of the variational formulation for quasi-static brittle fracture[END_REF] proves that, whenever the cracks are a priori known to propagate smoothly, the alternate minimization algorithm converges to the global minimizer of F ε,h for fine enough time discretization steps. In cases where cracks propagate brutally, this algorithm can only be proved to converge to critical points of F ε , which may be a local (or global) minimizers, but also saddle points for F ε . As per Remark 4.1.3, local minimizers of F ε can sometimes be proved to converge to local minimizers of F, but similar results are lacking at present in the case of saddle points, except in 1d. The detection of saddle points require a detailed stability study. Because of the typical size of the problems, this is a difficult task which has yet to be implemented. Instead, we investigate additional necessary conditions for minimality and propose to devise compatible algorithms.

The backtracking algorithm, a tool for global minimization. When cracks propagate brutally, the alternate minimization algorithm, or any other descent-based algorithm for that matter, cannot be expected to converge to the global minimizer of F ε,h . Indeed, a numerical method that relies solely on (4.2.1) will lead to evolutions whose total energy E(t) is not an absolutely continuous (or even continuous) function (see Figure 11 in [START_REF] Negri | A finite element approximation of the Griffith model in fracture mechanics[END_REF] or Figure 3(b) in Bourdin et al. (2000)). This is incompatible with Theorem 3.4.1. So, since (4.2.1) is satisfied at each time step, those evolutions have to correspond to local minimizers or saddle points of the regularized energy. Such solutions -spurious from the standpoint of global minimization -can actually be eliminated by enforcing an additional optimality condition.

Consider a monotonically increasing load, as in Section 1.2.3, and suppose the elastic energy density W to be 2-homogeneous (adapting this argument to p-homogenous W is trivial). If (ϕ i ε,h , v i ε,h ) is admissible for a time step t i , then t j /t i ϕ i ε,h , v i ε,h is admissible for all time steps t j with 0 ≤ j ≤ i, and

F ε,h t j t i ϕ i ε,h , v i ε,h = t 2 j t 2 i F b ε,h (ϕ i ε,h , v i ε,h ) + F s ε,h (v i ε,h ),
F b ε,h and F s ε,h denoting respectively the bulk and surface terms in F ε,h . But if the sequence {(ϕ i ε,h , v i ε,h )} is a solution of the fully discrete evolution, (ϕ j ε,h , v j ε,h ) must minimize F ε,h among all admissible pairs (ϕ, v), and in particular, for 0 ≤ j ≤ i ≤ n,

F b ε,h ϕ j ε,h , v j ε,h + F s ε,h v j ε,h ≤ t 2 j t 2 i F b ε,h (ϕ i ε,h , v i ε,h ) + F s ε,h (v i ε,h ). (4.2.2)
In establishing (4.2.2), the global minimality of the evolution {(ϕ i ε,h , v i ε,h )} was used, so that (4.2.2) is a necessary condition for global minimality but it is neither necessary, nor sufficient for local minimality. Since t j ≤ t i , the total energy, that is {F ε,h (ϕ i ε,h , v i ε,h )}, associated with an evolution satisfying (4.2.2) is monotonically increasing.

Algorithmically, we check condition (4.2.2) against all previous time steps t j , with j varying from 0 to i. If for some t j , (4.2.2) is not satisfied, then (ϕ j ε,h , v j ε,h ) cannot be the global minimizer for the time step t j , and (t j /t i ϕ i ε,h , v i ε,h ) provides an admissible field with a strictly smaller energy at time t j . In this case, we backtrack to time step t j , and restart the alternate minimization process, initializing the field v with v i ε,h . Because the alternate minimization algorithm constructs sequences with monotonically decreasing energy (at a given time step), repeated backtracking will converge to a solution such that (4.2.2) is satisfied for this particular choice of i and j.

The Tearing experiment

In order to illustrate the numerical method presented in the Sections above, we revisit the tearing experiment investigated in Section 2.2. Consider this time a rectangular domain Ω = (0, L) × (-H, H) (with traction free boundary conditions at x = L). The analysis in Subsection 2.2 still applies and the field constructed there under assumption (2.2.1) is an admissible test field for this problem, provided of course that 0 ≪ l(t) ≤ L.

However, when the domain has finite length, a crack splitting the whole domain is a minimizing competitor. Let ϕ c represent that solution. Following the notation in Section 2.2, we set S(ϕ c ) = (0, L) × {0} u c (t, x) = tH, so that E(ϕ c ) = kL.

A comparison of the energy of both types of evolutions demonstrates that, under assumption (2.2.1), the global minimizer for the tearing problem is such that u(x, y, t) = sign(y)u(t, x)e 3 and S(ϕ) = [0, l(t)) × {0}, with an evolution satisfies (necessary conditions for) (Ulm). We will comment further on this evolution in the discussion of That relation is used to compute the load t for which the crack length is l, and thereafter l(t).

u(x, t) =      tH 1 - x l(t) + if t ≤ L 2H
Once again, a crack splitting the whole domain along the x-axis is a minimizing competitor. Consider t c and l c := l(t c ) such that E b (t c , l c ) + kl c = kL. For t > t c , splitting the domain is energetically preferable. The value of t c can be estimated from the computations of E b (1, l). Using the finite element computations described above, we get t c ≃ .47. The critical length l c is such that

E b (1, l c ) = -(L -l c ) ∂E b ∂l (1, l c ).
Numerically, we obtain l c ≃ 2.28. That value is strictly less than the length l * c for which the constraint G ≤ k can no longer be met, as expected when global energy minimization presides. Indeed, the energetic landscape is explored in its entirety through global minimization, allowing the crack to decrease its energy at l c , rather than waiting for G to reach k at l * c .2 Notice the sudden jump introduced in Griffith's evolution -that satisfying (4.3.4) -at t c . Classically, such a jump would not be allowed to take place and Griffith's evolution would cease to hold at l * c . Figure 4.3.3 represents the numerically computed globally minimizing evolution of the bulk, surface, and total energies (thin lines), together with the analytically computed energies of the 1d solution -see (4.3.1), (4.3.2) -obtained in Section 2.2 and above (thick lines).

The computed evolution has the crack propagating smoothly for 0 ≤ t < t c , until it reaches the critical length l c , then cutting brutally through the domain. For small loads, the one-dimensional analysis overestimates the crack length; note that as l → 0, G(1, l) → ∞, and that the accuracy of our finite element computations cannot be guaranteed. When t, and therefore l, become large enough, the values of dE b (t)/dt and dE s (t)/dt become very close to those obtained in Section 2.2. Numerically we obtain dE b (t)/dt ≃ 7.43 (10) -2 and dE s (t)/dt ≃ 6.83 (10) -2 while the 1d result is dE b (t)/dt = dE s (t)/dt = H √ kµH ≃ 7.22 (10) -2 . Next, a numerical experiment that uses the algorithms developed in this section is conducted. So as to favor symmetric solution, we use a gave in Section 4.1 (k(1 + h/4ε)L = 6.25 (10) -2 ).

As noted before, the first evolution computed above using finite element analysis -that is that following Griffith until it jumps at l * c ≃ 4.19 -can be argued to be one satisfying (necessary conditions for) (Ulm). It propagates smoothly until it reaches l * c ≃ 4.19 at t = t * c ≃ .75, then brutally to the right end-side of the domain. Figure 4.3.6, represent the bulk, surface and total energies of this solution, compared to an experiment using the variational approximation and the alternate minimization, but without backtracking. Following the analysis in [START_REF] Bourdin | Numerical implementation of the variational formulation for quasi-static brittle fracture[END_REF], we expect that, as long as the crack propagates smoothly following local minimizers, the alternate minimization will provide the right evolution. When the crack propagates brutally, nothing can be said. However, once again, the agreement between our experiments is striking. Using the variational approximation, we obtain t * c ≃ .82 (instead of .75 using the finite element analysis). The estimate for the critical length is l * c ≃ 4.08 (vs. 4.19 for the finite element computations). The symmetry assumption about the x-axis was instrumental in deriving the theoretical results in Subsection 2.2; it was also imposed as a meshing restriction in the previous computation. In its absence, a bona fide theoretical prediction is difficult to make, but an educated guess may provide insight into the possible crack path. We thus introduce a third class of so-

  Figure 1.1.1. The cracked solid.

Figure 1

 1 Figure 1.1.2. Griffith dissipation potential.

  d H (A, B); = max{sup a∈A d(a, B), sup b∈B d(b, A)}. Examine for instance the initial time in the global minimality context with Γ

  Figure 2.2.1. Tearing.

  e -t ℓ ′ (tε) dt.

  (4.1.16) Since ℓ ′ (0) = lim s→0 ℓ(s)/s = 2H 1 (S(ϕ)), α ε = o(ε) and ∞ 0 e -t dt= 1, insertion of (4.1.16) into (4.1.15) and application of Lebesgue's dominated convergence theorem yields lim sup ε→0 Ω

(

  Fde) Find a sequence ϕ i+1 ε,h , v i+1 ε,h i=0,...,n of global minimizers for F ε,h under the constraints ϕ = g(t i+1 ) on Ω \ Ω and v = 0 on K i ε,h,η . (4.1.21)

Figure 4

 4 Figure 4.3.1. Tearing experiment: l → E b (1, l).

  Figure 4.3.6. The numerical values of E b (1, l) lead to an estimate of the position of the crack tip as a function of the load. Let ϕ(t, l) be the equilibrium deformation associated with the load t, and E b (t, l) := t 2 E b (1, l) the associated bulk

Figure 4

 4 Figure 4.3.2. Tearing experiment: l → G(1, l).

Figure 4

 4 Figure 4.3.6. Evolution of the bulk surface and total energies following (Ulm), as a function of the load t, computed using a classical finite element analysis. Comparison to the variational approximation without backtracking.

As an aside, note that the critical length lc does not depend upon the fracture toughness k!

where

This corresponds to a crack that propagates at constant speed dl dt = H µH k along the symmetry axis, until its length reaches L/2, and then jumps along the x-axis until the end point of that axis in the domain. Note that, during the smooth propagation phase, the bulk and surface energies of the sample are equal, and that, throughout the evolution, the total energy of the solution is E(t) = min 2tH µHk, kL .

(4.3.3)

Now, the anti-plane tearing problem is numerically solved by a method developed in [START_REF] Destuynder | Sur une interprétation mathématique de l'intégrale de rice en théorie de la rupture fragile[END_REF], then compared to the crack evolution analytically obtained above.

We consider a domain with dimensions H = 1, L = 5. The material properties are E = 1, ν = .2 (corresponding to µ ≃ .4167), k = 1.125 (10) -2 . Following Subsection 2.2, the analysis is restricted at first to symmetric solutions consisting of a single crack of length l(t) propagating along the x-axis, starting from the left edge of the domain, with l(0) = 0. In order to estimate l(t), we compute the equilibrium deformation ϕ(1, l) corresponding to a unit load and a crack of length l, using finite element meshes consisting of approximately 70,000 nodes, automatically refined around the crack tip. For various choices of l ∈ [0, L], we estimate the elastic energy E b (1, l) associated with ϕ(1, l), as well as the energy release rate G(1, l) = -∂E b /∂l(1, l), using classical formulae for the derivative of E b with respect to the domain shape. From now onward, we refer to the analytical solution as the "1d solution" in all figures, as well as in the text.

A quick analysis of the numerical results shows that G(1, l) is strictly decreasing (and therefore that E b is strictly convex) for 0 ≤ l < l * c , with l * c ≃ 4.19. For l * c ≤ l ≤ 5, G is an increasing function of l. Recalling Remark 1.2.3 in Section 1, we deduce that the crack will first propagate smoothly, following Griffith's criterion. When it reaches the length l * c , it will then jump brutally to the right edge of the domain because not doing so would violate the constraint that G ≤ k. It could be argued that such The mesh size is (10) -2 ; the regularization parameters are ε = (10) -2 and η ε = (10) -9 . We consider 100 equi-distributed time steps between 0 and 1. Recalling (4.1.20), the effective fracture toughness in the computations is k eff = (1 + h/4ε) k = .0125. Figure 4.3.4 represents the computed bulk, surface and total energy, as well as their values obtained via the proposed algorithm, as a function of t. Once again, the backtracking algorithm leads to an evolution with a monotonically increasing and continuous total energy.

Figure 4.3.5 represents the v field, representing the crack for t = .49 and t = .5. The values v = 0 are coded in red and v = 1 in blue.

The agreement with the classical solution is remarkable. The bulk energies are within 1% of each others, and the surface energies within 10%. For long enough cracks, the surface and bulk energies grow at a constant rate, and dE b (t)/dt ≃ 6.95 (10) -2 and dE s (t)/dt ≃ 7.03 (10) -2 . The critical load upon which the crack propagates brutally is .49 ≤ t c ≤ .5 (vs. a estimated value of .47), and the critical length is l c := l(.49) ≃ 2.46 which, again, is in agreement with the finite element analysis presented above (l c ≃ 2.28). The final surface energy is 6.38 (10) -2 , which is consistent with the estimate we lutions: a crack propagating along the symmetry axis with length l(t) until some critical t c at which it brutally bifurcates, reaching one of the sides of the domain. The crack for t ≥ t c is assumed L-shaped, i.e., of the form (0, l(t c )) × {0} ∪ {l(t c )} × (0, -H) or its mirror image with respect to the x-axis. It then remains to minimize in t c . Appealing to (4.3.3), (4.3.2) and comparing the energy associated with the straight crack, i.e., 2tH √ µHk, to that associated with the bifurcated crack, i.e., k(tH

The total energy of this branch of solution as a function of the loading parameter t is E(ϕ) = min 2tH µHk, 2kH .

If L > 2H, this asymmetric solution has a lower energy than its symmetric counterpart as soon as t ≥ k/µH. We propose a second set of experiments that use a non-symmetric Delaunay-Voronoi mesh. The mesh size is still h = (10) -2 , and the other parameters are those of the previous experiment. Once again, the position of the crack tip lags behind its theoretical position and the comparison between the numerical and theoretical energies is difficult.