# Comptes Rendus
## Mécanique

Jean-Jacques Marigo

L'approche variationnelle de la rupture : un exemple de collaboration fructueuse entre mécaniciens et mathématiciens

Volume 351, Numéro spécial S1 (2023), p. 591-613

En ligne depuis le 28 février 2023

Numéro spécial : L'héritage scientifique de Roland Glowinski

Rédacteurs en chef invités : Grégoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France), Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université) et Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

https://doi.org/10.5802/crmeca.170

<figure><img src="image_1.png" /></figure>

Cet article est publié sous la license

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL.

http://creativecommons.org/licenses/by/4.0/
where $k$ is a positive constant which corresponds to the volume dissipated energy during the damaging process. Furthermore the strain work $W$ becomes the following function of state:

$$W(\chi, \boldsymbol{\epsilon}) = \frac{1-\chi}{2} A_0 \boldsymbol{\epsilon} \cdot \boldsymbol{\epsilon} + \frac{\chi}{2} A_1 \boldsymbol{\epsilon} \cdot \boldsymbol{\epsilon} + k\chi,$$

where $\chi \in \{0, 1\}$. The consequence is that the evolution of the damage in a body during a loading process can be formulated in a variational form. Specifically, if $\mathbf{u}_t$ denotes the displacement field and $\chi_t$ the characteristic function of the damaged zone at time $t$, then $(\mathbf{u}_t, \chi_t)$ must satisfy

$$\mathcal{E}_t(\chi_t, \mathbf{u}_t) \le \mathcal{E}_t(\chi, \mathbf{u}_t), \quad \forall \chi \ge \chi_t, \quad \mathcal{E}_t(\chi_t, \mathbf{u}_t) \le \mathcal{E}_t(\chi_t, \mathbf{u}_t), \quad \forall \mathbf{u} \in \mathcal{V}$$

where $\mathcal{V}$ is the set of kinematically admissible displacement fields and $\mathcal{E}_t(\mathbf{u}, \chi)$ denotes the energy of the body $\Omega$ in the virtual state $(\mathbf{u}, \chi)$ at time $t$

$$\mathcal{E}_t(\chi, \mathbf{u}) = \int_{\Omega} W(\chi, \boldsymbol{\epsilon}(\mathbf{u})) dx - \mathbf{f}_t(\mathbf{u}),$$

$\boldsymbol{\epsilon}(\mathbf{u})$ being the symmetrized gradient of $\mathbf{u}$ and $\mathbf{f}_t$ denoting the current applied forces. But these partial minimization properties are not sufficient to determine $(\chi_t, \mathbf{u}_t)$. The idea developed in [3] was to reinforce the minimization principles and to introduce the so-called global stability condition

$$\mathcal{E}_t(\chi_t, \mathbf{u}_t) \le \mathcal{E}_t(\chi, \mathbf{u}) \quad \forall \chi \ge \chi_t, \quad \forall \mathbf{u} \in \mathcal{V}.$$

(The formulation must be complemented by the irreversibility condition on damage and an energy balance.) But it turns out that the problem has now no more solution and must be relaxed because in order to minimize the energy it can be more interesting to consider locally fine mixtures of the healthy and the damaged states of the material. Like in shape optimization [4], the relaxation procedure leads to use the optimal bounds for the effective elastic tensor coming from the mixture of the two states in a given volume fraction $D \in [0, 1]$. Finally, the relaxed strain work reads as

$$W(\boldsymbol{\epsilon}, D) = \varphi_*(\boldsymbol{\epsilon}, D) + kD,$$

where the relaxed elastic potential $\varphi_*(\boldsymbol{\epsilon}, D)$ can be obtained in a closed form in the case of an isotropic material, see [3]. So it appears that starting from a model where the damage variable can only take the values 0 or 1, one obtains finally a model where the damage variable can continuously vary from 0 to 1. Moreover the resulting continuous damage law comes from an optimal microstructuration of the material. It has not to be postulated, in contrast with the usual procedure of construction of damage laws.

In Section 2, Griffith's theory of brittle fracture is revisited. Adopting its (too) simple hypothesis that the surface energy is proportional to the area of the cracks, although it leads to unrealistic size effects and cannot account for the nucleation of cracks in a healthy body, one focusses on its key idea that cracks propagate on a principle of least energy [5]. That idea will be the mainstay of the revisited theory (the so-called variational approach to fracture [6]), providing fundamental tools both from a theoretical and numerical point of view. Specifically, in the restricted context of quasi-static processes and small deformations, at time $t$ of a loading process, the following energy $\mathcal{E}_t(\mathbf{u}, \Gamma)$ is associated to the displacement field $\mathbf{u}$ and the crack state of a body $\Omega$:

$$\mathcal{E}_t(\mathbf{u}, \Gamma) = \mathcal{P}_t(\mathbf{u}, \Gamma) + \mathcal{S}(\Gamma),$$

where $\mathcal{P}_t(\mathbf{u}, \Gamma)$ and $\mathcal{S}(\Gamma)$ are the potential energy and the surface energy respectively, that is

$$\mathcal{P}_t(\mathbf{u}, \Gamma) = \int_{\Omega \cap \Gamma} \frac{1}{2} A (\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}_t^0) \cdot (\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}_t^0) dx - \mathbf{f}_t(\mathbf{u}), \quad \mathcal{S}(\Gamma) = \int_{\Gamma} G_c dS$$
where A is the stiffness tensor, $\boldsymbol{\epsilon}(\mathbf{u})$ the symmetric gradient of $\mathbf{u}$, $\boldsymbol{\epsilon}_t^0$ and $\mathbf{f}_t$ being the prestrain field and the applied forces at time $t$. After solving the elastic problem giving the displacement field, the energy becomes a function of $\Gamma$ alone at given $t$:

$$ \hat{\mathcal{E}}_t(\Gamma) = \mathcal{E}_t(\mathbf{u}_t(\Gamma), \Gamma) = \min_{\mathbf{u} \in \mathcal{V}} \mathcal{P}_t(\mathbf{u}, \Gamma) + \mathcal{S}(\Gamma). $$

Then, following Griffith's idea, it was proposed in [7] that the current crack state $\Gamma_t$ satisfies the following global stability criterion:

$$ \hat{\mathcal{E}}_t(\Gamma_t) \le \hat{\mathcal{E}}_t(\Gamma), \forall \Gamma \supset \Gamma_t $$

where the crack states $\Gamma$ which test the stability of $\Gamma_t$ must be greater than it because of the irreversibility of the crack path. (The formulation must be completed by an energy balance but it is not necessary to give its precise statement here, see [8] for the details.) This global stability criterion is in general too strong and must be replaced by a local one. (To be short, the crack states $\Gamma$ testing the stability of the state $\Gamma_t$ must be “close” to $\Gamma_t$ and of course the difficulty is to make precise that definition of neighborhood.) However the global stability criterion remains useful because the most important rigorous results [9, 10] have been established with it.

In any case, from a practical point of view the real issue to apply that principle was to generate arbitrary crack paths. The solution came from Ambrosio and Tortorelli who proved in [11, 12] that, in the context of image segmentation, the Mumford-Shah functional, which is similar to the mechanical energy in Griffith theory, can be approximated by a functional which can be interpreted as an energy of a damage model with an additional term containing the damage gradient. This fundamental result is established by Gamma-convergence [13] (and hence once more by global minimization) where the small parameter corresponds to a length which fixes the size of the damage localisation. It was possible to actually transpose this idea in damage and fracture mechanics. Finally, a typical energy functional which is used is the so-called AT1 model. It reads as

$$ \mathcal{E}_t(\mathbf{u}, D) = \int_{\Omega} W_t(\boldsymbol{\epsilon}, D, \nabla D) dx - \mathbf{f}_t(\mathbf{u}) $$

where $W_t$ is the energy density given by

$$ W_t(\boldsymbol{\epsilon}, D, \nabla D) = \frac{1}{2}(1-D)^2 A_0 (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_t^0) \cdot (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_t^0) + \frac{3G_c}{8} \left( \frac{D}{\ell} + \ell \nabla D \cdot \nabla D \right) $$

where $A_0$ is the stiffness tensor of the healthy material and $\ell$ is a length. From a theoretical point of view, starting from this form of the energy, the damage evolution process is governed by an irreversibility condition, a stability criterion and an energy balance [8, 14, 15]. In particular the local stability criterion is the heir of Griffith's energy minimisation principle and reads as

$$ \mathcal{E}_t(\mathbf{u}_t, D_t) \le \mathcal{E}(\mathbf{u}^*, D^*), \quad \forall (\mathbf{u}^*, D^*) \text{ close to } (\mathbf{u}_t, D_t), 1 \ge D^* \ge D_t $$

where the notion of neighbourhood is given by the $H^1(\Omega)$ norm. By standard arguments of the calculus of variations, one deduces that the local damage criterion reads

$$ \frac{1}{(1 - D_t)^3} S_0 \boldsymbol{\sigma}_t \cdot \boldsymbol{\sigma}_t + \frac{3G_c \ell}{4} \Delta D_t \le \frac{3G_c}{8\ell}, $$

where $S_0$ is the healthy compliance tensor and $\boldsymbol{\sigma}_t$ is the current stress. So it contains a term in laplacian $D$ which comes from the damage gradient terms of the energy. Therefore, the maximal uniaxial stress that the material can sustain is given by $\sigma_c = \sqrt{\frac{3G_c E_0}{8\ell}}$, and it can be shown [16] that the optimal damage profile to create a crack is parabolic with a width equal to $4\ell$.

The variational formulation is also very useful from a numerical point of view. Specifically, after time discretisation, the displacement and damage states at time step $i$ are obtained by an alternate minimisation algorithm: the energy at time step $i$ is minimised with respect to $\mathbf{u}$ at given $D$, then with respect to $D$ (under the constraint of irreversibility) at fixed $\mathbf{u}$ up to convergence [6].
Moreover each minimisation admits a unique solution and is easy to implement by virtue of the convexity and the quadratic character of the energy with respect to each variable.

In a first step, the gradient damage model was used as an approximation of Griffith model, in the spirit of Ambrosio-Tortorelli. In that case, the length has a purely numerical meaning and can be chosen as small as desired. The advantage with respect to the original Griffith formulation is that it becomes possible to manage arbitrary crack paths in space and time with classical finite elements. Several examples have illustrated those qualities [6, 17, 18].

But by considering the length $\ell$ as an arbitrary small parameter one falls in the defect of the Griffith theory and one cannot explain the nucleation of a crack because the critical stress tends to infinity as $1/\sqrt{\ell}$. Therefore it becomes necessary to consider that the length $\ell$ is a characteristic length of the material which can be identified from the measurement of $E_0$, $\sigma_c$ and $G_c$ by virtue of the relation coming from the optimal damage profile:

$$ \ell = \frac{3G_cE_0}{8\sigma_c^2}. $$

And even if the Gamma-convergence result can no more be used because one does no more pass to the limit in $\ell$, one practically recovers the Griffith theory as soon as the characteristic material length is small by comparison with the length of the crack [19]. The thermal shock experiment considered in [20] and whose more important aspects are shown in Section 3 is an illuminating example of application of the variational approach to fracture, see also [21] for the details. The problem of the nucleation of a crack at the tip of a notched specimen is a fundamental test for the ability of gradient damage models to predict with a good precision the critical load at which the phenomenon occurs. That has been done in [22] by using the AT1 model where the numerical predictions are compared with experiments made for several materials and a range of notch angles. Here also the agreement is very good.

In Section 3, we also consider the tearing of a strip cut from a sheet of anisotropic paper according to the “trousers test” shown in Figure 9. From the experimental point of view, the goal of the experiment is to identify the function $G_c(\theta)$ giving the Griffith surface energy density as a function of the orientation of the crack in the sheet of paper. From the theoretical point of view this leads to solving the issue of the kinking of a preexisting crack and to find what is the criterion giving the angle of the kink. Assuming that the strip is perfectly flexible and inextensible, there is no more elastic energy in the system and the potential energy reduces to the opposite of the work of the prescribed force at point A. Accordingly, there is no other strategy else than to follow the variational approach and to use the stability criterion that it contains. It is what is made in [23,24]. It turns out that the local stability criterion imposes that the orientation of the crack after the kink is given by the angle which maximises the ratio $G/G_c$ between the potential energy release rate and the surface energy density. Adopting this criterion allows the authors to obtain all or a part of the function $G_c(\theta)$ from the measurement of the critical force and the direction of the kink for various orientations of the strip with respect to the paper sheet. A complete discussion can be found in [23].

## Remarque introductive

Le texte de cet article ne suit pas le standard des textes scientifiques usuellement publiés dans les revues internationales, le style formel de rigueur cédant parfois le pas à un discours à caractère plus personnel et narratif. C'est l'objectif même de cet article, qui se propose de retracer l'historique de l'approche variationnelle de la rupture, qui le demande. En conséquence, dans ces parties plus personnelles et donc plus subjectives, le discours se fera à la première personne du singulier ou du pluriel suivant qu'il implique son auteur seul ou la compagnie de proches. Dans les parties plus classiques, le ton impersonnel conventionnel sera adopté.
Cet article est aussi l'occasion pour moi de rendre hommage à Roland Glowinski et à travers lui à l'école française de mathématiques appliquées qui a joué un rôle essentiel dans l'éclosion de l'approche variationnelle de la rupture. Au début des années 80, le jeune ingénieur-chercheur d'EDF que j'étais a profité de la liberté totale que lui offrait sa hiérarchie pour affiner sa culture mathématique, besoin qu'il ressentait pour mieux appréhender les questions de modélisation mécanique qui lui étaient posées. Le Duvaut-Lions, le Glowinski-Lions-Trémolières, l'Ekeland-Temam ou le Bensoussan-Lions-Papanicolaou le disputaient au Germain vert ou au Truesdell version Euvrard. Ils venaient enrichir la formation de mécanique théorique que j'avais reçue via le cours de Yves Bamberger des Ponts ou le DEA de mécanique théorique de Paris 6 dans lequel Georges Duvaut avait eu l'audace d'inscrire un cours sur la théorie de l'homogénéisation alors qu'elle était encore en pleine construction. Il était partout question de variationnel (sauf dans le Truesdell, malheureusement... pour lui), que ce soit pour la construction des modèles, l'établissement de leurs propriétés ou leur traitement numérique. What else? Comme Gilles Francfort suivit quelques années plus tard un itinéraire semblable, il ne fut pas difficile, lorsque l'occasion s'en présenta, que nous convergions vers une approche de la rupture où il serait question de minimisation, de stabilité ou d'homogénéisation. Il est probable que sans cet environnement fertile, l'idée n'aurait jamais germée.

# 1. La genèse : le modèle d'endommagement brutal partiel

L'histoire commence en mars 1989, à Évian, lors du colloque annuel Mécamat. Gilles Francfort y avait été invité à faire une conférence sur l'optimisation de forme. La communauté des mécaniciens des matériaux français découvrit à cette occasion qu'en général un problème d'optimisation de forme n'admet pas de solution classique avec une partie du domaine susceptible d'être occupé vide et la partie complémentaire remplie du matériau choisi, mais qu'il fallait envisager en chaque point matériel du domaine un mélange fin de vide et de matière. On est alors ramené à des problèmes d'homogénéisation et à la recherche de répartition optimale de la microstructure vide-matière pour, par exemple, maximiser la rigidité effective à proportion de matière donnée [4, 25]. Le fait que cette optimisation des propriétés effectives fasse appel aux bornes d'Hashin-Shtrikman ravissait évidemment les micromécaniciens qui consacraient une grande partie de leurs efforts à étudier le comportement effectif des matériaux composites et pour lesquels ces bornes étaient un peu le Graal. Mais de mon point de vue, l'intérêt se situait ailleurs. Suite à son exposé, nous eûmes avec Gilles (Francfort) une discussion portant sur la modélisation de l'endommagement par des approches micro-mécaniques. À l'époque, cela consistait à se donner un type de microstructure (micro-fissures, micro-vides) dépendant d'une famille de paramètres, puis via des techniques d'homogénéisation d'accéder au comportement macrospécique qui en résultait et enfin de postuler une loi d'évolution de ces paramètres microstructuraux en se basant éventuellement sur des formalismes de type matériaux standards généralisés (c.f. par exemple [26, 27]). Le reproche que faisait Gilles Francfort à cette approche était qu'elle était trop contraignante envers la microstructure dont on imposait la forme *a priori*. Il aurait préféré qu'on la laisse libre et que son développement procède d'un principe d'optimisation type minimisation d'énergie. Cette critique rejoignait du reste celle qui était faite en optimisation de forme par les tenants de l'optimisation topologique aux adeptes de l'optimisation paramétrique ou géométrique. Nous nous mîmes alors d'accord pour essayer de faire fonctionner son idée sur le problème modèle de l'endommagement dit brutal partiel dont les principaux ingrédients sont rappelés ci-dessous.

On considère un matériau dont le comportement est linéairement élastique avec une élasticité caractérisée par le tenseur de rigidité $A_0$ tant qu'on ne le déforme pas beaucoup, mais qui s'endommagera, son tenseur de rigidité prenant alors la valeur $A_1$, si les déformations franchissent
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Réponse uniaxiale d'un matériau à endommagement brutal partiel.</figcaption></figure>

un seuil. Cette perte de rigidité est « irréversible », le matériau gardant dès lors la rigidité $A_1$. Pour simplifier l'exposé on supposera que $A_0$ est strictement supérieur à $A_1$ et que $A_1$ est strictement positif (ces inégalités devant être comprises au sens des formes quadratiques), ce qui veut dire qu'il y a effectivement perte de rigidité lors de l'endommagement du matériau, mais que cette perte n'est pas totale : donc endommagement « brutal partiel ». Dans [3] il est postulé que le domaine d'élasticité $\mathcal{R}$, qui définit les déformations que le matériau peut supporter sans s'endommager, est de la forme

$$ \mathcal{R} = \{ \varepsilon : \frac{1}{2} (A_0 - A_1) \varepsilon \cdot \varepsilon \le k \}, \quad (1) $$

où $k$ est une constante positive donnée qui s'interprète comme l'énergie dissipée (par l'élément de volume) lors de son endommagement. Dans un essai uniaxial, cela conduit à une réponse du matériau qui dans un diagramme contrainte-déformation est du type de celle représentée figure 1.

Or, cette forme du critère d'endommagement n'a pas besoin d'être postulée, elle découle naturellement d'un principe physique plus fondamental qui est le postulat de Drucker-Ilyushin [1,2]. Ce postulat dit que lorsqu'un élément de volume subit un cycle en déformation $\mathcal{C}$, le travail de déformation durant ce cycle est nécessairement non négatif. Formellement :

$$ \int_{\mathcal{C}} \sigma \cdot d\varepsilon \ge 0, \quad \forall \mathcal{C} $$

où $\sigma$ désigne le tenseur des contraintes qui vaut donc $A_0\varepsilon$ ou $A_1\varepsilon$ suivant l'état d'endommagement du matériau à l'instant considéré du cycle. Le postulat est indépendant du deuxième principe de la thermodynamique et comme on va le voir il est en général beaucoup plus restrictif sur les lois de comportement. Comme toute la structure variationnelle du problème d'évolution de l'endommagement découle de cette forme du critère d'endommagement, il est important de montrer comment on la déduit du postulat de Drucker-Ilyushin. (De plus, la démonstration qui est courte n'a jamais été publiée. On trouvera des implications du postulat en plasticité, endommagement continu ou rupture dans [28].)

Supposons simplement que le domaine d'élasticité $\mathcal{R}$ est un sous-ensemble connexe fermé à frontière « régulière » de l'ensemble des tenseurs d'ordre 2 symétriques (se reporter à [28] pour une définition plus précise des propriétés minimales que doit avoir $\mathcal{R}$ pour justifier les raisonnements qui suivent). Notons $\mathcal{R}'$ l'ensemble défini dans le membre de droite de (1). Il s'agit de montrer que $\mathcal{R} = \mathcal{R}'$. Pour cela considérons le cycle en déformation suivant, le matériau étant dans son état sain au départ : (i) on fait passer la déformation de $\varepsilon_0$ qui est à l'intérieur de
R à ε₁ qui est sur la frontière de R par un chemin contenu dans R. Donc durant cette phase, la rigidité reste à sa valeur initiale A₀; (ii) on endommage le matériau en faisant un incrément de déformation hε*, h > 0 petit, depuis ε₁ dans une direction ε* qui le fait sortir de R; (iii) le matériau étant désormais endommagé, on revient à la déformation initiale ε₀. Calculons le travail de déformation durant les trois phases, en les notant Wᵢ, 1 ≤ i ≤ 3, et en n'explicitant que les termes d'ordre 0 par rapport à h :

$$
W_1 = \frac{1}{2} A_0 \varepsilon_1 \cdot \varepsilon_1 - \frac{1}{2} A_0 \varepsilon_0 \cdot \varepsilon_0, \quad W_2 = O(h), \quad W_3 = \frac{1}{2} A_1 \varepsilon_0 \cdot \varepsilon_0 - \frac{1}{2} A_1 \varepsilon_1 \cdot \varepsilon_1 + O(h).
$$

En vertu du postulat, leur somme doit être non négative. En passant à la limite en $h$, on en déduit

$$
\frac{1}{2}(A_0 - A_1)\epsilon_0 \cdot \epsilon_0 \le \frac{1}{2}(A_0 - A_1)\epsilon_1 \cdot \epsilon_1. \quad (2)
$$

Cette inégalité vaut pour tout ε₀ intérieur à R et tout ε₁ sur le bord de R. L'inégalité reste vraie
si l'on fait tendre ε₀ vers le bord. Mais alors, en intervertissant les rôles de ε₀ et ε₁, on voit que
½(A₀ - A₁)ε · ε doit être une constante (nécessairement positive) k pour tous les ε sur le bord
de R. De plus en reportant cette propriété dans (2), on en déduit que ½(A₀ - A₁)ε · ε ≤ k pour
tout les ε dans R. On a donc obtenu que R ⊂ R' et que ∂R ⊂ ∂R', ce qui n'est possible que
si R = R', compte-tenu des propriétés supposées de R et des propriétés effectives de R'. Ce
résultat est évidemment remarquable puisqu'il dit que le postulat sur le travail de déformation
impose que l'endommagement n'évolue que quand la restitution d'énergie élastique atteint une
valeur critique k caractéristique du matériau.

Ce faisant, en notant $\chi \in \{0, 1\}$ l'état d'endommagement du matériau, il est facile de vérifier que si l'on part de l'état $(\chi = 0, \varepsilon = 0)$ et si l'on déforme progressivement l'élément de volume pour l'amener à l'état $(\chi, \varepsilon)$, alors le travail de déformation ne dépend que de l'état final (indépendamment du chemin suivi) et s'écrit :

$$
W(\chi, \varepsilon) = \frac{1-\chi}{2} A_0 \varepsilon \cdot \varepsilon + \frac{\chi}{2} A_1 \varepsilon \cdot \varepsilon + k\chi. \quad (3)
$$

De plus, si le matériau part de son état sain et si on lui impose une déformation ε, alors l'état
d'endommagement du matériau χ qui en résulte s'obtient en minimisant W(χ*, ε) par rapport à
χ* :

$$
\chi = \operatorname{argmin}_{\chi_* \in \{0,1\}} W(\chi_*, \varepsilon). \quad (4)
$$

Comme dans le même temps, à état d'endommagement donné, les déplacements d'équilibre d'une structure chargée s'obtiennent par minimisation de l'énergie potentielle, on a donc obtenu une structure variationnelle pour le problème d'évolution de l'endommagement.

La condition d'irréversibilité de l'endommagement qui exige que χ croisse dans le temps demande que l'on procède de façon incrémentale. Formulons le problème d'évolution après discrétisation en temps du processus. En se plaçant en petites déformations, supposons qu'une structure Ω, ouvert de Rⁿ, constituée du matériau à endommagement brutal partiel soit soumise à un processus de chargement caractérisé par la forme linéaire fᵢ dépendant du pas de temps i et par l'ensemble des champs de déplacements cinématiquement admissibles V (supposé pour simplifier sous-espace vectoriel fermé fixe de H¹(Ω, Rⁿ)). On peut associer à la structure, à chaque pas de temps i, l'énergie, qui dépend de son état d'endommagement défini par la fonction caractéristique de la zone endommagée χ et le champ de déplacement u ∈ V, donnée par

$$
\mathcal{E}_i(\chi, \mathbf{u}) = \int_{\Omega} W(\chi, \mathcal{E}(\mathbf{u})) \, dx - f_i(\mathbf{u}) \quad (5)
$$

où $\mathcal{E}(\mathbf{u})$ désigne le gradient symétrique de $\mathbf{u}$. Ce que l'on a obtenu jusqu'ici est que, connaissant l'état d'endommagement $\chi_{i-1}$ au pas de temps $i-1$, l'état d'endommagement $\chi_i$ et de déplacement $\mathbf{u}_i$ au pas $i$ doivent vérifier les problèmes de minimisation partielle de $\mathcal{E}_i$ :

$$
\mathcal{E}_i(\chi_i, \mathbf{u}_i) \le \mathcal{E}_i(\chi, \mathbf{u}_i), \quad \forall \chi \ge \chi_{i-1}, \qquad \mathcal{E}_i(\chi_i, \mathbf{u}_i) \le \mathcal{E}_i(\chi_i, \mathbf{u}), \quad \forall \mathbf{u} \in \mathcal{V}.
$$
Le problème est que ceci ne suffit pas à caractériser $\chi_i$ et $\mathbf{u}_i$. Il faut donc introduire une condition plus forte et c'est là qu'intervient l'idée de Gilles (Francfort) : on va minimiser par rapport au « couple » $(\chi, \mathbf{u})$. Ce principe va constituer ce qui sera appelé ensuite dans les approches variationnelles la « condition de stabilité ». On peut la formuler en temps continu et elle s'écrit alors

$$ \mathcal{E}_t(\chi_t, \mathbf{u}_t) \le \mathcal{E}_t(\chi, \mathbf{u}) \quad \forall \chi \ge \chi_t, \quad \forall \mathbf{u} \in \mathcal{V}, \qquad (6) $$

étant entendu que le champ test $\chi$ ne prend ses valeurs que dans $\{0, 1\}$.

En fait, c'est une simple extension à un système dissipatif *rate independent* du principe de stabilité basé sur l'énergie pour les systèmes conservatifs. Mais cette extension exige de tenir compte de la condition d'irréversibilité. Elle dit que la structure ne pourra se trouver à un instant donné dans un certain état que s'il n'existe pas d'état accessible depuis cet état ayant une énergie moindre. Autrement dit, la stabilité d'un état est testée en envisageant des perturbations possibles de cet état et en comparant les énergies. Dans le cas présent, l'accessibilité est restreinte par la condition d'irréversibilité : la structure ne peut tester la stabilité qu'en envisageant des états d'endommagement plus grands (au sens large) que l'état testé. L'autre point essentiel est de savoir si l'on se limite à des états proches de l'état testé ou si on autorise tous les tests possibles. Dans le premier cas on parle de stabilité « locale », dans le deuxième de stabilité « globale ». Telle qu'elle est formulée dans (6) il s'agit de stabilité globale. C'est beaucoup plus confortable sur le plan mathématique car on dispose de plus d'outils pour la traiter et cela évite d'avoir à choisir une norme pour définir la notion de voisinage. (Ici une norme naturelle serait $H^1(\Omega)$ pour les déplacements et $L^1(\Omega)$ pour les fonctions caractéristiques.) Mais c'est beaucoup plus discutable sur le plan physique, car l'accession à un état de moindre énergie peut nécessiter de franchir des barrières. De plus, l'énergie peut ne pas être bornée inférieurement, ce qui n'est pas le cas ici car l'endommagement n'est que partiel mais ça le devient dans le cas de fissures et de chargement à forces imposées. Toutefois pour ce problème modèle, nous nous en tiendrons à (6).

En revenant à la formulation en temps discret, en supposant que la structure est saine à l'instant initial ($\chi_0 = 0$) et en considérant le premier pas de temps, la condition de stabilité devient simplement un problème de minimisation de $\mathcal{E}_1(\chi, \mathbf{u})$ parmi toutes les fonctions caractéristiques et tous les champs de déplacement dans $\mathcal{V}$. Mais comme cette énergie du fait de sa non convexité n'est pas semi-continue inférieurement ce problème de minimisation n'admet pas de solution classique (comme en optimisation de forme), il faut le relaxer. Cette relaxation consiste à envisager des mélanges fins, en chaque point matériel, du matériau sain et du matériau endommagé. En unidimensionnel ($n = 1$), l'énergie relaxée est la convexifiée de l'énergie d'origine. Le coefficient d'élasticité effectif résultant du mélange en proportion $D$ du matériau endommagé avec le matériau sain est donné par la moyenne harmonique des coefficients d'élasticité $A_0$ et $A_1$, soit

$$ \frac{1}{A_*(D)} = \frac{1-D}{A_0} + \frac{D}{A_1} $$

et l'énergie relaxée qui dépend maintenant du champ d'endommagement $D$ qui peut prendre ses valeurs dans tout l'intervalle $[0, 1]$ s'écrit

$$ \mathcal{E}_1^*(D, u) = \int_{\Omega} \left( \frac{1}{2} A_*(D) (u')^2 + kD \right) dx - f_1(u). $$

En dimension supérieure à 1, le calcul de la relaxée fait appel comme en optimisation de forme aux bornes optimales des tenseurs d'élasticité effectifs provenant d'un mélange du matériau sain et du matériau endommagé. En dimension 2, dans le cas de comportements isotropes, l'expression du potentiel élastique relaxé $\varphi_*(D, \boldsymbol{\epsilon})$ (qui n'est plus quadratique en $\boldsymbol{\epsilon}$) est donnée dans [3], l'expression de l'énergie relaxée s'en déduisant

$$ \mathcal{E}_1^*(D, \mathbf{u}) = \int_{\Omega} (\varphi_*(D, \boldsymbol{\epsilon}) + kD) \, dx - f_1(\mathbf{u}). $$
De plus, on peut montrer que, même s'il n'y a pas unicité, on peut obtenir le mélange optimal par mises en couches (éventuellement multiples) des deux phases. La microstructuration de l'endommagement a donc transformé le modèle de départ qui était un endommagement brutal en un endommagement continu avec comme variable d'endommagement macroscopique la fraction locale de matériau endommagé. Et c'est la condition de stabilité (qui traduit l'idée initiale de Gilles) qui force la microstructuration et qui exige son optimalité. On peut y voir là une avancée conceptuelle majeure dans la modélisation de l'endommagement.

Quand on passe aux autres pas de temps et dès le deuxième, les choses se compliquent à cause de la condition d'irréversibilité et du processus de relaxation. On ne peut se contenter de raisonner sur la fraction D de matériau endommagé et d'exiger seulement que $D_2 \ge D_1$, il faut encore s'assurer que les mises en couches successives qui sont nécessaires pour obtenir l'optimalité sont compatibles avec la condition d'irréversibilité. Il faut enfin passer à la limite quand le pas de temps tend vers 0 pour obtenir la réponse en temps continu. Ce travail délicat qui relève de la bravoure a été réalisé par Gilles Francfort et Adriana Garroni dans [29]. Il montre bien toutes les difficultés supplémentaires qu'introduit l'irréversibilité et que l'on ne saurait esquiver.

On s'est longtemps demandé si ce modèle d'endommagement brutal partiel pouvait servir de modèle de rupture. La réponse n'a été apportée que très récemment par Jean-François Babadjian, Flaviana Iurlano et Filip Rindler, et elle est négative, c.f. [30]. Je vais me contenter ici de commenter ce résultat. Bien sûr, pour répondre à la question, il faut envisager un endommagement total, i.e. $A_1 = 0$, mais on ne peut pas le faire directement car alors l'énergie relaxée est nulle. Il faut faire tendre $A_1$ vers 0, puis passer à la limite. Mais dans le même temps il faut toucher à la densité d'énergie dissipée $k$ et la faire tendre vers l'infini de façon à ce que le résultat ne soit pas trivial. Introduisons donc un petit paramètre $\eta$ destiné à tendre vers 0 et notons $A_1^\eta$ et $k^\eta$ la rigidité du matériau endommagé et la densité d'énergie dissipée qui en dépendent. Une analyse rapide montre que pour obtenir un résultat non trivial (une énergie relaxée non nulle et non purement élastique) il faut que si $A_1^\eta$ est de l'ordre de $\eta$, alors $k^\eta$ soit de l'ordre de $1/\eta$. Posons donc $A_1^\eta = \eta \bar{A}_1$ et $k^\eta = \bar{k}/\eta$.

Dans le cas unidimensionnel ($n = 1$), il est facile de conclure car l'expression de la relaxée est explicite (c'est la convexifiée) et fait intervenir la moyenne harmonique $A_*^\eta$ des modules d'élasticité. Le travail de déformation relaxé s'écrit donc

$$W^\eta(D, \varepsilon) = \frac{1}{2} A_*^\eta(D) \varepsilon^2 + k^\eta D.$$

On en déduit par dérivation le critère d'endommagement. Exprimé en déformation il s'écrit

$$-\frac{dA_*^\eta(D)}{dD} \varepsilon^2 \le 2k^\eta,$$

et donc exprimé en contraintes

$$\frac{d}{dD} \left( \frac{1}{A_*^\eta(D)} \right) \sigma^2 \le 2k^\eta,$$

puisque $\sigma = A_*^\eta(D)\varepsilon$. En explicitant les termes on obtient

$$\left( \frac{1}{\eta \bar{A}_1} - \frac{1}{\bar{A}_0} \right) \sigma^2 \le \frac{2 \bar{k}}{\eta}.$$

On voit donc que la contrainte au seuil est indépendante de D, ce qui est un résultat attendu puisqu'elle correspond à la ligne de Maxwell dans la construction de la convexifiée, mais aussi que cette constante tend vers la valeur finie non nulle $\sigma_0 = \sqrt{2\bar{k}\bar{A}_1}$ quand $\eta$ tend vers 0. On a donc obtenu un modèle limite qui est toujours un modèle d'endommagement, dans lequel la rigidité décroît progressivement vers 0, mais sans que les contraintes tendent vers 0. La réponse contrainte-déformation est celle représentée dans la figure 2.
<figure><img src="image_3.png" /><figcaption>**FIGURE 2.** Modèle unidimensionnel d'endommagement limite obtenu en faisant tendre $A_1$ vers 0 et $k$ vers l'infini.</figcaption></figure>

En plusieurs dimensions, la démonstration est d'une toute autre difficulté, mais la conclusion est la même : il est montré dans [30] que le travail de déformation relaxé limite a un composition homogène de degré 1 à l'infini en terme des déformations, le domaine des contraintes admissibles tendant vers un convexe K qui est explicité et ne se réduit pas à 0.

Donc, ce type de modèles d'endommagement brutal ne peut pas espérer rendre compte de la propagation de fissures « à la Griffith », du moins en l'état. Mais là n'était pas son objectif. Il avait constitué un galop d'essai pour tester la pertinence d'une approche de l'endommagement et de la rupture par minimisation d'énergie. Le vrai défi allait être de faire fonctionner cette idée dans le cadre de la théorie de Griffith. C'est la tâche à laquelle nous nous attelâmes avec Gilles dès son recrutement comme professeur à l'Université Paris-Nord en 1994. Nous fûment rapidement re joints par Blaise Bourdin qui commença les premiers calculs dès son stage de DEA en 1995.

## 2. La théorie de Griffith revisitée

Griffith publia, il y a maintenant un siècle, un article [5] jetant les bases d'une théorie de la rupture conçue pour les matériaux fragiles qui dès lors influencera tous les travaux de modélisation qui seront entrepris sur cette problématique. Elle s'appuie sur deux hypothèses clés, toutes deux de nature énergétique :

(1) l'énergie à apporter pour créer de nouvelles surfaces fissurées est proportionnelle à ces surfaces;

(2) les fissures se propagent dans une structure suivant un principe de minimisation d'énergie.

La première hypothèse peut partiellement se justifier à l'échelle atomique à partir de potentiels interatomiques et correspond en première approximation à l'énergie de surface fournie par deux plans atomiques infiniment éloignés l'un de l'autre. Elle permet de n'introduire qu'une constante matériau supplémentaire dans le cas d'un milieu isotrope : la densité d'énergie de surface $G_c$. Il s'avérera que cette hypothèse simplificatrice s'accompagne d'effets secondaires néfastes en ne permettant pas d'expliquer des phénomènes fondamentaux comme la nucléation de fissures dans une structure saine, la propagation par fatigue de fissures existantes ou en prédisant des effets d'échelle manifestement irréalistes. C'est elle qui va conduire à développer la mécanique de l'endommagement pour pallier les lacunes de la théorie de Griffith. Mais adoptons-la tout de même ici.

C'est la deuxième hypothèse, énoncée de façon vague dans [5], qui constitue le coeur du problème et qui demande à être formalisée. Elle contient en germe l'idée développée dans le
modèle d'endommagement brutal, à savoir qu'il ne faut pas imposer de restriction *a priori* sur le trajet spatio-temporel des fissures, mais les considérer comme surfaces de discontinuité libres dont l'évolution va être régie par des principes énergétiques et la condition d'irréversibilité.

Venons-en donc à une formulation plus précise de la loi d'évolution de la fissuration dans le cas de structures bi ou tridimensionnelles constituées de matériaux linéairement élastiques et isotropes. En suivant l'idée de Griffith et en se limitant au cadre quasi-statique des petites déformations, on peut associer, à chaque instant $t$, à un champ de déplacement $\mathbf{u}$ et à un état de fissuration $\Gamma$ l'énergie $\mathcal{E}_t(\mathbf{u}, \Gamma)$ de la structure $\Omega$ soumise à des forces extérieures (surfaciques ou volumiques) caractérisées par la forme linéaire $\mathbf{f}_t$ :

$$ \mathcal{E}_t(\mathbf{u}, \Gamma) = \mathcal{P}_t(\mathbf{u}, \Gamma) + \mathcal{S}(\Gamma), \qquad (7) $$

$\mathcal{P}_t(\mathbf{u}, \Gamma)$ et $\mathcal{S}(\Gamma)$ étant respectivement l'énergie potentielle et l'énergie de surface,

$$ \mathcal{P}_t(\mathbf{u}, \Gamma) = \int_{\Omega \setminus \Gamma} \frac{1}{2} A (\boldsymbol{\varepsilon}(\mathbf{u}) - \boldsymbol{\varepsilon}_t^0) \cdot (\boldsymbol{\varepsilon}(\mathbf{u}) - \boldsymbol{\varepsilon}_t^0) dx - \mathbf{f}_t(\mathbf{u}), \quad \mathcal{S}(\Gamma) = \int_{\Gamma} G_c dS \qquad (8) $$

où A désigne le tenseur de rigidité, $\boldsymbol{\varepsilon}(\mathbf{u})$ le gradient symétrisé du champ de déplacement $\mathbf{u}$ et $\boldsymbol{\varepsilon}_t^0$ la prédéformation imposée localement à l'instant $t$ (typiquement la déformation thermique via la température). En supposant connu l'état de fissuration $\Gamma$ et à condition bien sûr que les efforts extérieurs et l'état de fissuration soient compatibles avec l'équilibre, on peut déterminer le champ de déplacement réalisant l'équilibre à l'instant $t$, ce champ $\mathbf{u}_t(\Gamma)$ dépendant à la fois de $t$ et de $\Gamma$. Ce faisant, l'énergie devient une fonctionnelle de l'état de fissuration $\Gamma$ :

$$ \hat{\mathcal{E}}_t(\Gamma) = \mathcal{E}_t(\mathbf{u}_t(\Gamma), \Gamma) = \min_{\mathbf{u} \in \mathcal{V}} \mathcal{P}_t(\mathbf{u}, \Gamma) + \mathcal{S}(\Gamma). $$

Il reste à préciser le principe de minimisation. Le choix fort qui est d'abord proposé dans [7], dans l'esprit de [3] mais aussi de [5], est de demander que l'état de fissuration à l'instant $t$ corresponde à celui qui conduit à la plus petite énergie de la structure parmi tous les états possibles qui respectent la condition d'irréversibilité. Cette condition de stabilité « globale » s'écrit alors :

$$ \hat{\mathcal{E}}_t(\Gamma_t) \le \hat{\mathcal{E}}_t(\Gamma), \quad \forall \Gamma \supset \Gamma_t. $$

Il faut lui adjoindre en plus un bilan d'énergie pour compléter la formulation, mais ce point qui n'est pas essentiel pour la discussion qui suit ne sera pas abordé ici. Le lecteur intéressé par la formulation complète peut se référer à Alexander Mielke [8] pour son énoncé dans le cadre général des processus *rate independent* et à [6, 14, 15] pour ce qui concerne l'endommagement et la rupture. Comme il a été remarqué dans le cas du modèle d'endommagement brutal partiel, ce critère de stabilité globale est physiquement contestable. Mais ici, ça l'est encore plus car il ne peut pas fonctionner si $\mathbf{f}_t$ n'est pas nul car alors, en général, l'énergie n'est pas bornée inférieurement. (Dit de façon grossière, il suffit d'envisager des fissures qui séparent en plusieurs morceaux la structure et d'envoyer ces morceaux à l'infini en faisant travailler les forces imposées.) Autrement dit, on ne peut le faire fonctionner que dans le cas où les données non nulles sont des déplacements imposés ou des prédéformations. Ce qui évidemment en limite singulièrement la portée. Mais considérons tout de même ce critère de stabilité globale malgré les défauts qui ont été soulignés. Deux points importants sont à noter :

* Des résultats d'existence de solution à ce problème d'évolution ont été obtenus, le plus récent et le plus complet étant celui d'Antonin Chambolle et de Vito Crismale [9]. (Bien sûr un tel résultat d'existence n'est possible que dans le cas où $\mathbf{f}_t = 0$.) En le comparant avec le modèle d'endommagement brutal partiel où un tel résultat n'existait pas car il avait fallu relaxer l'énergie, cela signifie qu'avec le modèle de Griffith on ne peut pas avoir de mélanges fins entre les fissures et le matériau élastique. On ne peut donc pas parler de matériau micro-fissuré.
* En étudiant le problème d'évolution de fissures de Griffith dans des structures composites élastiques et en se limitant au cas anti-plan, Alessandro Giacomini et Marcello Ponsiglione ont montré dans [10] que les fissures ne perturbaient pas le résultat d'homogénéisation classique pour les matériaux élastiques. Autrement dit, le problème d'évolution limite obtenu lorsqu'on fait tendre la taille de la cellule de périodicité vers 0 est celui de fissures de Griffith ayant en général une densité d'énergie de surface anisotrope se propageant dans un matériau homogène linéairement élastique dont le tenseur d'élasticité effectif est celui de l'homogénéisation classique (sans fissure dans la cellule). Ici en-core, on ne pourra pas parler de cellule fissurée. On ne peut donc espérer construire des modèles d'endommagement à partir de micro-fissures régies par le modèle de Griffith.

Ces deux résultats remarquables tiennent pour l'essentiel à l'hypothèse de Griffith sur l'énergie de surface. Essayons de comprendre pourquoi. Pour ce qui est de l'homogénéisation, considérons un réseau de cellules de taille η dans un milieu tridimensionnel qui contiendra donc un nombre de cellules de l'ordre de 1/η³. Avec l'hypothèse de Griffith, si l'on met dans chaque cellule une fissure de taille de l'ordre de η, cela coûtera une énergie de surface de l'ordre de η² par cellule et donc de l'ordre de 1/η en tout ce qui est beaucoup trop cher. On peut au plus envisager de « contaminer » un nombre de cellules de l'ordre de 1/η² si l'on veut que l'énergie de surface totale soit de l'ordre de 1 comme l'énergie potentielle dans le cas d'un chargement d'ordre 1. Ce faisant, la grande majorité des cellules ne voient pas les fissures et se comportent donc comme s'il n'y en avait pas (à des effets de couche limite près qui ne se voient pas à l'ordre dominant mais qu'aux ordres d'après). D'un autre côté, si l'on met dans chaque cellule des fissures de taille de l'ordre de η³/² de façon à avoir une énergie de surface totale d'ordre 1, ces « microfissures » ne se voient pas du point de vue de l'énergie élastique. Donc, dans tous les cas, le comportement élastique effectif n'est pas perturbé. On peut utiliser pratiquement les mêmes arguments pour ce qui est des mélanges fins : les fissures de la taille de la microstructure coûtent trop cher, les fissures plus petites ne rapportent rien.

Il suffit de remplacer l'énergie de surface de Griffith par une énergie de surface de type Barenblatt [6, 31] où la densité n'est pas une constante mais dépend du saut des déplacements, pour qu'aucun des deux résultats précédents reste vrai. Dans le cas d'une énergie de surface de type Barenblatt, il faut nécessairement relaxer l'énergie car il y aura mélanges fins entre le matériau et des surfaces de discontinuités avec des micro-sauts. Comme ce n'est pas tant la taille des fissures qui comptent que la taille des sauts de déplacement pour une énergie de surface de type Barenblatt, on peut mettre des fissures dans chaque cellule et il peut y avoir interaction entre la fissuration et l'homogénéisation.

Comme on a raisonné seulement en ordre de grandeur des énergies, on sent bien que ces propriétés ne sont pas dues au choix d'un critère de stabilité globale et qu'ils resteront vrais même si on l'affaiblit. L'adoption de ce critère et les outils de Gamma-convergence qui l'accompagnent ont au moins eu le mérite de les mettre en évidence. Cette qualité va se confirmer tout de suite quand on va aborder la question de la mise en oeuvre pratique de ce principe de surfaces de discontinuité libres qui cherchent à minimiser l'énergie. En effet, se pose immédiatement la question : comment peut-on envisager de faire tester par la structure tous les états de fissuration possibles ? La réponse va être fournie encore une fois par les méthodes directes du Calcul des Variations en approximant la fonctionnelle énergie de Griffith par une fonctionnelle énergie qui s'interprètera comme un modèle d'endommagement à gradient. L'idée et la première réalisation de ce type d'approximation sont dues à Luigi Ambrosio et Vincenzo Maria Tortorelli [11, 12]. Comme ils avaient été établis dans le cadre de la segmentation d'image, autrement dit pour des problèmes scalaires sans dépendance en temps, il s'agissait de les transposer à la rupture et donc à des problèmes vectoriels avec nécessairement une dépendance en temps à cause de l'irréversibilité de la fissuration. Formellement cette "régularisation elliptique", pour reprendre la
terminologie de [11], consiste à remplacer la fonctionnelle de Griffith (7)-(8) par la fonctionnelle régularisée suivante, appelée désormais AT2 dans la littérature :

$$ \mathcal{E}_{\ell}(\mathbf{u}, D) = \mathcal{P}_{\ell}(\mathbf{u}, D) + \mathcal{S}_{\ell}(D) \qquad (9) $$

avec

$$ \mathcal{P}_{\ell}(\mathbf{u}, D) = \int_{\Omega} \frac{1}{2} (1 - D)^2 A(\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}^0) \cdot (\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}^0) dx - \mathbf{f}(\mathbf{u}) \qquad (10) $$

et

$$ \mathcal{S}_{\ell}(\Gamma) = \int_{\Omega} \frac{G_c}{2} \left( \frac{D^2}{\ell} + \ell \nabla D \cdot \nabla D \right) dx, \qquad (11) $$

où la dépendance en temps de $\boldsymbol{\epsilon}^0$ et $\mathbf{f}$ est implicite. Dans la fonctionnelle régularisée, il n'y a plus de surfaces de discontinuité, le champ de déplacement est continu et les fissures sont remplacées par le champ scalaire $D$ qui prend ses valeurs dans l'intervalle $[0, 1]$. Ce champ qui affecte le tenseur de rigidité dans l'expression de l'énergie potentielle s'interprète donc comme un champ d'endommagement. L'énergie $\mathcal{S}_{\ell}(D)$ représente l'énergie dissipée par endommagement, elle est donc "volumique" et remplace l'énergie de surface de Griffith. Notons qu'elle contient un terme en gradient d'endommagement qui, comme on va le voir, joue un rôle essentiel.

En l'absence de forces extérieures, i.e. quand $\mathbf{f} = \mathbf{0}$, les données non nulles se limitant donc à des prédéformations ou des déplacements imposés, il est possible de retrouver le résultat d'Ambrosio-Tortorelli et de montrer que la fonctionnelle régularisée Gamma-converge vers la fonctionnelle de Griffith lorsque le paramètre $\ell$, qui a ici la dimension d'une longueur, tend vers 0, c.f. [6, 32]. Ce résultat de Gamma-convergence signifie en particulier que les minima de $\mathcal{E}_{\ell}$ convergent vers le minimum de $\mathcal{E}$ et que les minimiseurs $(u_{\ell}, D_{\ell})$ de la fonctionnelle régularisée convergent (dans un sens approprié) vers les minimiseurs $(u, \Gamma)$ de la fonctionnelle de Griffith quand $\ell$ tend vers 0 (voir [13] pour une introduction aux notions de Gamma-convergence). Sans faire la démonstration, on peut comprendre pourquoi on obtient ce résultat de convergence en analysant le rôle de chaque terme dans l'expression de la fonctionnelle régularisée : (i) pour que le terme en $D^2/\ell$ dans l'énergie dissipée n'explose pas quand $\ell$ tend vers 0 il faut que $D$ soit nul "presque partout"; (ii) d'un autre côté, si l'on veut que l'énergie potentielle, qui varie comme le carré du chargement à $D$ fixé, ne devienne pas trop grande, on a intérêt à ce que $D$ soit voisin de 1 dans des bandes, nécessairement minces en raison de (i), afin que la structure s'assouplisse; (iii) ces bandes endommagées ne peuvent pas être trop minces à cause du terme en gradient de $D$ dans l'énergie dissipée. L'épaisseur optimale est de l'ordre de $\ell$ et elle peut se calculer explicitement (ainsi que le profil d'endommagement optimal) en résolvant un problème unidimensionnel, c.f. [6, 16]. (Sans ce terme en gradient, on pourrait localiser l'endommagement sur des bandes infiniment minces, ce qui permettrait de casser la structure sans que cela coûte de l'énergie.) Donc à la limite, on obtient des surfaces à travers lesquelles $\mathbf{u}$ est discontinu, l'énergie dissipée est proportionnelle à l'aire de ces surfaces de discontinuités et le coefficient $G_c/2$ dans l'énergie dissipée est ajusté de façon à retrouver l'énergie de Griffith. Notons au passage que ce résultat règle du même coup la question du lien entre la mécanique de la rupture et la mécanique de l'endommagement qui a préoccupé la communauté durant les années 70, 80 et 90. Le lecteur intéressé par l'historique de cette problématique en trouvera une présentation dans [33].

Du point de vue numérique cette régularisation ne présente que des avantages. Tout d'abord, on n'a plus à gérer des surfaces de discontinuité libres et on peut utiliser des éléments finis classiques. C'est le champ d'endommagement qui en se localisant va indiquer quand et comment les fissures vont se développer dans la structure au cours du temps, sans que l'on ait besoin de faire des hypothèses *a priori* sur leur trajet spatio-temporel. On a même intérêt à utiliser des maillages non structurés pour ne pas biaiser les calculs et favoriser certains trajets. De plus, du fait du caractère variationnel de la formulation, on peut utiliser un algorithme de minimisation alternée qui consiste, à chaque pas de temps, à minimiser $\mathcal{E}_{\ell}$ par rapport à $\mathbf{u}$ à $D$ fixé, puis
par rapport à *D* à **u** fixé, jusqu'à convergence. Avec le modèle AT2, comme l'énergie est une fonction quadratique et convexe de chaque champ pris séparément (alors que l'énergie n'est pas une fonction convexe du couple (*u*, *D*)), ces problèmes de minimisation admettent une solution unique et sont faciles à implémenter : pour **u** cela revient à résoudre un problème d'élasticité linéaire, alors que pour *D* on doit minimiser la fonctionnelle quadratique sous la contrainte d'irréversibilité, *D* ne pouvant que croître au cours du temps. Il est facile de vérifier que c'est un algorithme de descente, l'énergie décroissant à chaque itération. Toutefois, cela n'assure pas qu'il converge vers un minimum local (et encore moins un minimum global) et un post traitement peut s'avérer nécessaire pour juger de la qualité de l'état convergé. En particulier, un algorithme de *backtracking* a été développé par Blaise Bourdin [18] pour (re)trouver le minimum global (quand il existe).

D'un point de vue historique, ces modèles d'endommagement à gradient vont être utilisés dans un premier temps, i.e. à la fin des années 90 et au début des années 2000, comme modèles régularisés de Griffith : le modèle de Griffith est alors considéré comme le bon modèle physique, le modèle d'endommagement en constituant simplement une approximation via le petit paramètre $\ell$ et le résultat de Gamma-convergence. L'objectif était simplement de s'en servir comme outil capable de rendre compte de trajets de fissuration complexes, mais pas de rendre compte de tout le processus de fissuration allant de la nucléation jusqu'à la rupture finale de la structure. C'est Blaise Bourdin qui a été le maître d'oeuvre de ce travail d'implémentation et de tests. On peut trouver quelques exemples dans [6, 17, 18] ou encore sur sa page web (https://ms.mcmaster.ca/~bourdinb/defectmechanics/). Un de ces exemples est présenté dans la section 3.

Toutefois, en procédant ainsi, la longueur $\ell$ n'ayant pas de signification physique précise puisque c'est un paramètre purement numérique que l'on peut (doit) prendre aussi petit que l'on veut (peut), on perd la possibilité de rendre compte correctement de la nucléation des fissures. En effet, comme la contrainte maximale que le matériau peut supporter tend vers l'infini comme $1/\sqrt{\ell}$ dans le modèle d'endommagement qui sert d'approximation, on retombe sur les défauts de la théorie de Griffith. À la limite, le matériau peut supporter n'importe quelle contrainte et la nucléation devient impossible. Il faut donc changer de point de vue et considérer que le modèle d'endommagement est le modèle physique et que la longueur $\ell$ est une longueur caractéristique du matériau qu'il s'agira d'identifier expérimentalement (en général à partir de mesures indirectes). Ce faisant, le modèle de Griffith n'en devient plus qu'une approximation qui régira l'évolution de la fissuration une fois que les fissures auront été créées par le processus de localisation de l'endommagement et que leur longueur deviendra grande devant la longueur caractéristique $\ell$ du matériau [19]. Mais alors, comme le critère de nucléation va porter sur les contraintes et donc du modèle d'endommagement choisi, on a intérêt à envisager des modèles plus généraux que le modèle AT2. C'est en substance ce qui a été proposé dans [14, 15] où est présentée la procédure générale de construction de modèles d'endommagement à gradient avec le problème d'évolution basé sur une condition d'irréversibilité, un principe de stabilité et un bilan d'énergie. Contentons-nous ici d'en extraire une version simplifiée et d'en commenter les principaux ingrédients.

(1) *Choix de l'énergie.* À l'instant $t$, l'énergie totale de la structure dans l'état ($\mathbf{u}$, $D$) s'écrit

$$ \mathcal{E}_t(\mathbf{u}, D) = \int_{\Omega} W_t(\boldsymbol{\epsilon}(\mathbf{u}), D, \nabla D) \, dx - \mathbf{f}_t(\mathbf{u}), $$

avec

$$ W_t(\boldsymbol{\epsilon}, D, \nabla D) = \frac{1}{2} \alpha(D) A(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_t^0) \cdot (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_t^0) + G_c \left( \frac{\beta(D)}{\ell} + \gamma \ell \nabla D \cdot \nabla D \right), \quad (12) $$

ce qui constitue une famille de modèles d'endommagement à gradient qui se distinguent par le choix des fonctions (scalaires) $\alpha(D)$ et $\beta(D)$ et par le choix des constantes $\ell$ et
$\gamma$, en supposant que le tenseur d'élasticité $A$ et la densité d'énergie de surface $G_c$ sont des données matériaux. Les deux fonctions et les deux constantes doivent satisfaire un certain nombre de conditions qui ne sont pas précisées ici pour que ce modèle puisse être considéré comme un modèle "compatible" avec le modèle de rupture de Griffith qui doit en être une limite aux grandes échelles, i.e. quand $\ell$ devient petit par rapport aux longueurs des fissures.

(2) *Principe de stabilité locale.* Comme on l'a déjà souligné à plusieurs reprises, le critère de stabilité globale doit être remplacé par un critère de stabilité locale, ne serait-ce que parce que l'énergie n'est pas bornée inférieurement dès lors que $f_t$ n'est pas nul. Cette condition de stabilité se formule donc comme suit : à l'instant $t$, les champs de déplacement et d'endommagement $(u_t, D_t)$ doivent être tels que

$$ \mathcal{E}_t (\mathbf{u}_t, D_t) \le \mathcal{E} (\mathbf{u}^*, D^*), \quad \forall (\mathbf{u}^*, D^*) \text{ voisin de } (\mathbf{u}_t, D_t), D_t \le D^* \le 1 $$

avec la norme $H^1(\Omega)$ comme choix naturel de la norme fixant la notion de voisinage. Cette condition consiste à comparer l'énergie de la structure dans son état actuel avec l'énergie qu'elle aurait si on la perturbait "un peu", sachant que l'on n'a droit d'utiliser, pour des raisons d'irréversibilité, que des états d'endommagement plus grands que l'état actuel pour tester sa stabilité.

(3) Enfin il faut inclure un *bilan d'énergie* dont on ne se servira pas ici, c.f. [14, 15] pour sa formulation précise.

On déduit facilement du principe de stabilité le critère d'endommagement local qui correspond à une condition de stabilité d'ordre 1 par rapport au petit paramètre de perturbation et qui, exprimé en contraintes, s'écrit :

$$ -\frac{1}{2} \frac{\alpha'(D_t)}{\alpha(D_t)^2} A^{-1} \boldsymbol{\sigma}_t \cdot \boldsymbol{\sigma}_t + 2\gamma\ell G_c \Delta D_t \le \frac{G_c}{\ell} \beta'(D_t), \quad (13) $$

le terme en laplacien du champ d'endommagement provenant du terme en gradient de l'énergie. Dans les zones encore saines, i.e. là où $D = 0$, on déduit de (13) la condition locale que doivent satisfaire les contraintes. Compte-tenu du fait que $\alpha(D)$ décroît de 1 à 0, elle s'écrit

$$ A^{-1} \boldsymbol{\sigma}_t \cdot \boldsymbol{\sigma}_t \le \frac{2\beta'(0)G_c}{|\alpha'(0)|\ell}. \quad (14) $$

On déduit également de (13) le profil de la zone endommagée qui accompagne la création d'une fissure, l'endommagement valant 1 sur les lèvres de la fissure et décroissant vers 0 dans la direction normale, c.f. [16].

En pratique, les deux modèles les plus utilisés ont été le modèle AT2, déjà évoqué, qui correspond à

$$ \alpha(D) = (1 - D)^2, \quad \beta(D) = \frac{D^2}{2}, \quad \gamma = \frac{1}{2} $$

et le modèle AT1, sa variante, qui correspond à

$$ \alpha(D) = (1 - D)^2, \quad \beta(D) = \frac{3D}{8}, \quad \gamma = \frac{3}{8}. $$

Ils ne diffèrent essentiellement l'un de l'autre que par le fait que $\beta(D)$ est linéaire en $D$ dans AT1, et quadratique en $D$ dans AT2 (d'où leur nom). Les constantes sont ajustées de façon que l'énergie dissipée pour créer une fissure de surface unitaire soit dans les deux cas égale à $G_c$. Le modèle AT2 présente un certain nombre de défauts. Comme $\beta'(0) = 0$, on déduit de (14) qu'il ne possède pas de phase élastique, l'endommagement apparaissant à la moindre contrainte. Ceci induit une pollution sur le plan numérique, toute la structure s'endommageant (un peu) instantanément. De plus, on peut tirer de (13) qu'il présente une phase durcissante avant une phase adoucissante, les contraintes commençant par augmenter avant de décroître. Ce phénomène est visible sur la figure 3. Le modèle AT1 corrige ces défauts : il présente une phase élastique, puis est "purement"
adoucissant. On tire de (13) et (14) que la contrainte uniaxiale maximale que le matériau peut
supporter est

$$
\sigma_c = \sqrt{\frac{3G_c E}{8\ell}}
$$

où E représente le module d'Young du matériau sain. On peut se servir de cette relation pour identifier la longueur caractéristique du matériau à partir des mesures du module d'Young, de la contrainte critique et de la densité d'énergie de surface. La comparaison entre les modèles d'endommagement dans un essai uniaxial est faite dans la figure 3.

<figure><img src="image_4.png" /><figcaption>**FIGURE 3.** Réponse du matériau lors d'un essai uniaxial avec des charges et des décharges pour les modèles AT1 (à gauche) et AT2 (à droite) pour les mêmes valeurs de E, G<sub>c</sub> et l', σ<sub>c</sub> = √(3G<sub>c</sub>E/8l', ε<sub>c</sub> = σ<sub>c</sub>/E. En rouge les phases où le matériau s'endommage, en vert celles où il reste élastique.</figcaption></figure>

En adoptant ce point de vue qui considère que le modèle d'endommagement est le bon modèle physique, il devient possible de rendre compte qualitativement et quantitativement de tout le processus de rupture allant de la nucléation de fissures jusqu'à la rupture des structures. Le problème du choc thermique présenté dans la section 3 en est à ce jour la meilleure illustration. De même, le traitement de la nucléation en fond d'entaille avec les modèles AT1 ou AT2 ont permis de rendre compte du phénomène et la comparaison avec des résultats expérimentaux faits sur différentes éprouvettes avec une famille de valeurs de l'angle d'entaille et pour plusieurs matériaux différents est excellente, c.f. [22].

**3. Trois exemples d'illustration**

3.1. Propagation de fissures dans du verre par choc thermique

Lorsqu'on remplit un ballon de verre d'un liquide froid, apparaissent un réseau de fissures réparties de façon quasi-périodique qui commencent par se propager le long d'un trajet oscillant avant de se propager de façon rectiligne, le front de fissure suivant l'évolution de la surface libre du liquide. On obtient à la fin un état de fissuration comme celui qui est représenté à gauche sur la figure 4 [34].

De même si on trempe des lames de verre dans un bain de liquide froid on observe la nucléation et la propagation d'une ou plusieurs fissures, leur nombre et leur trajet dépendant de la vitesse de plongée dans le bain et du gradient initial de température entre la lame et le bain. On peut voir la propagation droite d'une seule fissure, la propagation oscillante avec
des amplitudes voire des périodes d'oscillations variant progressivement au cours de l'essai ou encore la propagation de plusieurs fissures, c.f. les trois exemples à droite de la figure 4 [35].

<figure><img src="image_5.png" /></figure>

**FIGURE 4.** Résultats expérimentaux de fissuration par choc thermique : à gauche, réseau de fissures oscillantes ou rectilignes dans un ballon chaud rempli progressivement d'un liquide froid [34] ; à droite, état de fissuration final de lames de verre chaudes plongées dans un bain froid, le résultat dépendant de la vitesse de plongée et du gradient de température entre la lame et le liquide [35].

<figure><img src="image_6.png" /></figure>

**FIGURE 5.** Traitement numérique du problème de fissuration par choc thermique de lames de verre : résultats numériques obtenus par Blaise Bourdin (c.f. son site web https://ms.mcmaster.ca/~bourdinb/defectmechanics/oscillatingcracks/) à partir du modèle d'endommagement à gradient AT2 en utilisant un algorithme de minimisation alternée et un maillage non structuré. Les quatre figures correspondent à différentes valeurs des paramètres de l'essai (vitesse de plongée et gradient de température).

En utilisant le modèle AT2 avec les valeurs des paramètres (hormis $\ell$) calées sur les données expérimentales du matériau, Blaise Bourdin est arrivé à reproduire une partie de ces différents trajets de fissuration en fonction des deux paramètres de l'essai. Par exemple, en utilisant un maillage non structuré, après avoir nucléé une fissure par un procédé purement numérique,
il a obtenu à partir de l'algorithm de minimisation alternée des trajets comme ceux qui sont représentés sur les quatre images de la figure 5. On remarquera le plus à droite où l'on voit des amorces de bifurcation de fissures qui finalement avortent. Ce problème de choc thermique, qui fut un des premiers à être traité par l'approche variationnelle et dont les résultats ont été insuffisamment diffusés, mériterait d'être revisité aujourd'hui de façon plus approfondie tant du point de vue théorique que numérique (et même du point de vue expérimental si l'on en croit les spécialistes de la question).

## 3.2. Choc thermique d'une céramique

Le problème du choc thermique traité avec ce type de modèles à gradient d'endommagement en montre toute la richesse et la pertinence puisqu'il permet de traiter tout le processus de développement de la fissuration. Rappelons brièvement les points marquants qui sont développés dans [20, 21]. On se place en contraintes planes et on considère une éprouvette rectangulaire de largeur $L$ dont la température est $T_0$. Elle est soumise à un choc thermique à partir de l'instant $t = 0$, le bord supérieur $x_2 = 0$ étant dès lors maintenu à la température $T_1 < T_0$. Par diffusion thermique et en supposant que l'endommagement ne perturbe pas cette propagation, le chargement est donc donné par le champ de déformation thermique dépendant du temps $\epsilon_t^0 = \alpha(T_1 - T_0)\text{Erfc}(\frac{x_2}{2\sqrt{k_f}})$ où $\alpha$ et $k$ sont les coefficients de dilatation et de diffusion thermiques, Erfc étant la fonction d'erreur complémentaire. Le déplacement normal $u_1$ des bords latéraux $x_1 = \pm L/2$ étant bloqué, ce refroidissement induit une contrainte de traction $\sigma_{11}$ dans le matériau qui, s'il restait élastique, vaudrait $\sigma_T = E_0\alpha(T_0 - T_1)$ en $x_2 = 0$. Par conséquent si le choc thermique est assez grand, i.e. si la contrainte $\sigma_T$ est plus grande que la contrainte critique $\sigma_c$ du matériau, alors le bord supérieur de l'éprouvette va s'endommager dès l'instant 0. En se plaçant dans cette situation, on décrit dans la figure 6 ci-dessous l'évolution de l'endommagement prédit par l'algorithm de minimisation alternée en utilisant le modèle AT1 avec un maillage triangulaire non structuré et en partant d'une éprouvette saine.

<figure><img src="image_7.png" /><figcaption>**FIGURE 6.** Différentes phases de l'évolution de l'endommagement lors du choc thermique.</figcaption></figure>

On peut distinguer trois phases : (i) L'endommagement s'initie dès l'instant 0 sur le bord supérieur de l'échantillon $x_2 = 0$, puis se propage progressivement vers l'intérieur tout en restant
uniforme en $x_1$ jusqu'à un instant critique où apparaît une bifurcation avec un mode de bifurcation périodique; (ii) L'endommagement devient dépendant de $x_1$ et tout en augmentant se localise progressivement au centre de chaque période pour donner finalement naissance à des fissures; (iii) Une fois que le réseau de fissures a été établi, et que leur longueur et leur espacement sont suffisamment grands par rapport à la longueur caractéristique du matériau $\ell$, ces fissures se propagent en suivant pratiquement la loi de Griffith : elles avancent en reproduisant toujours le même profil parabolique d'endommagement qui correspond à une énergie dissipée $G_c$ par unité de longueur. On observe un phénomène de doublement de période à intervalles réguliers : une fissure sur deux s'arrêtant alors que la voisine continue à se propager.

Des expériences ont été faites par [36] sur des éprouvettes de céramique (largeur $L = 50$ mm, hauteur = 9.8 mm, épaisseur = 1 mm) pour plusieurs valeurs de choc thermique. (Dans ces expériences les 4 côtés de l'éprouvette sont mis à la température $T_1$.) Les calculs sont faits avec le modèle AT1, la longueur caractéristique $\ell$ se déduisant des mesures de $E_0$, $\sigma_c$ et $G_c$, ce qui donne ici $\ell = 46 \ \mu m$. La comparaison avec ces résultats expérimentaux est montrée sur les figures 7 et 8.

<figure><img src="image_8.png" /><figcaption>**FIGURE 7.** Comparaison des états finaux de fissuration observés et calculés après un choc thermique de 380 °C sur une éprouvette en céramique.</figcaption></figure>

On notera le très bon accord sur la morphologie du réseau de fissures à la fin de l'essai, que ce soit visuellement sur la figure 7 où sont comparés les états finaux des éprouvettes après un choc thermique de 380 °C ou de façon plus quantitative sur la figure 8 qui montre la valeur moyenne des espacements des fissures en fonction de la profondeur de fissuration. La longueur $\ell_0 = G_c/\alpha^2(T_0 - T_1)^2$ représentative de l'intensité du choc thermique, sert de référence. Dans le cas du choc de 380 °C, elle vaut $14 \ \mu m$.

<figure><img src="image_9.png" /><figcaption>**FIGURE 8.** Comparaison des espacements moyens des fissures observés ou calculés à la fin de l'essai en fonction des profondeurs de fissuration.</figcaption></figure>
### 3.3. Déchirure d'une feuille de papier anisotrope

Au laboratoire Physique et Mécanique des Milieux Hétérogènes, un groupe de chercheurs autour de Benoît Roman s'intéressent au problème de la déchirure de feuilles minces dans des situations diverses et variées. En particulier, ils ont étudié dans [23, 24] la déchirure d'une feuille de papier « anisotrope ». En supposant que ce problème relève de la théorie de Griffith, il s'agit de mesurer expérimentalement la densité d'énergie de surface qui dépend ici de l'orientation de la fissure par rapport aux axes de la feuille du fait de l'anisotropie du papier. Le procédé expérimental est celui schématisé sur la figure 9 : (i) une bande rectangulaire est découpée dans une feuille de papier; (ii) cette bande est préfissurée sur une longueur $l_0$; (iii) le point B étant fixé, on exerce au point A une force orientée suivant $e_1$ dont l'intensité F croît progressivement; (iv) on mesure l'orientation de la fissure qui s'amorce à la pointe initiale $P_0$ lorsque la force atteint une force critique $F_c$. Puis on recommence l'expérience en changeant l'orientation de la bande découpée dans la bande de papier. Si on note $\alpha$ l'angle de la bande par rapport à la feuille de papier et $\theta$ l'angle de branchement que fait la fissure qui s'amorce avec la fissure initiale (la fissure branchée fera donc l'angle $\theta + \alpha$ avec la feuille de papier), il s'agit de déduire de la mesure de cet angle $\theta$ et de la force critique $F_c(\alpha)$ (qui dépend de l'orientation de la bande) la valeur $G_c(\theta + \alpha)$ de la densité d'énergie de surface. D'un point de vue théorique on est donc confronté à un problème de branchement de fissure dans un matériau anisotrope. Il est traité dans [23] dans le cadre de l'approche variationnelle en s'appuyant sur des résultats qui avaient été obtenus par [37]. Sans refaire tous les calculs et tous les raisonnements, on va revoir les grandes lignes de la démarche suivie.

<figure><img src="image_10.png" /><figcaption>**FIGURE 9.** Déchirure d'une feuille de papier pré-fissurée : à gauche, configuration de référence plane de la bande orientée d'un angle $\alpha$ par rapport à la feuille de papier et contenant une fissure ayant branché d'un angle $\theta$ par rapport à la fissure initiale; à droite, configuration d'équilibre sous l'effet de la force $Fe_1$ appliquée en A.</figcaption></figure>

D'un point de vue énergétique, la bande de papier étant mince, on peut négliger en première approximation son énergie de flexion. De plus, comme les forces exercées jusqu'à la déchirure ne sont jamais grandes, on peut donc négliger aussi en première approximation son énergie d'extension. Ce faisant, on a affaire à une membrane inextensible parfaitement flexible et l'énergie potentielle se réduit à l'opposé du travail de la force exercée considérée comme charge morte. Ceci oblige de plus à travailler en grands déplacements avec un matériau anisotrope. On est donc très loin du cadre habituellement adopté en mécanique de la rupture fragile. On ne voit pas comment on peut s'affranchir d'une approche variationnelle. Aussi le problème est-il traité dans [23] dans ce cadre-là et il s'avère que le critère de stabilité locale suffit.

On commence par remarquer qu'à état de fissuration donné (en n'envisageant que des courbes « simples » avec une seule « pointe »), les positions de A, B et de la pointe de la fissure P « dans la configuration d'équilibre » sont alignées suivant l'axe $e_1$. Par conséquent, du fait de
l'inextensibilité, l'énergie potentielle ne dépend que de la position de la pointe de la fissure dans la configuration de référence plane et est donnée par

$$ \mathcal{P} = -F(\ell_1 + \ell_2), $$

où $\ell_1$ et $\ell_2$ désignent les distances de P à A et B dans cette configuration de référence. Il est alors facile par des calculs purement géométriques d'en déduire le taux de restitution d'énergie potentielle dû à un branchement de la fissure d'un angle $\theta$. En utilisant le critère de stabilité locale on obtient le critère de Griffith qui s'écrit

$$ 2F \cos\theta \leq G_c(\theta + \alpha)h \quad (15) $$

où $h$ est l'épaisseur de la feuille.

D'un point de vue théorique, si on suppose connue la fonction $G_c(\cdot)$ on en déduit que la force critique d'amorçage du branchement dans la direction $\theta$ est celle pour laquelle l'égalité a lieu dans (15) (sous la condition évidemment qu'elle soit positive). En envisageant toutes les valeurs de $\theta$ possibles, c'est la plus petite de ces valeurs de la force critique qui correspondra à la force d'amorçage et l'angle $\theta^*(\alpha)$ correspondant sera l'angle de branchement. Donc

$$ F_c(\alpha) = \min_{\theta 1 \in (-\pi/2, +\pi/2)} \frac{G_c(\theta + \alpha)h}{2 \cos\theta}. $$

Le critère qui donne l'angle de branchement n'est donc pas celui du taux de restitution d'énergie potentielle maximal, mais dans le cas présent d'un milieu anisotrope c'est le rapport $G/G_c$ qui compte et c'est donc, à chargement fixé, « le taux de restitution d'énergie potentielle par taux de création d'énergie de surface » qui doit être maximal. Ce résultat se déduit directement du critère de stabilité, alors que les critères invoqués usuellement basés sur les singularités ou les symétries sont ici sans objet.

D'un point de vue expérimental, de la mesure de $F_c(\alpha)$ et de $\theta^*(\alpha)$ on en déduit $G_c(\theta^*(\alpha) + \alpha)$. En faisant varier $\alpha$ on en déduit toute ou partie de la fonction $G_c(\cdot)$ suivant les propriétés de cette dernière. Une discussion est faite dans [23] à ce sujet.

## Remarque finale

En parcourant ce bref historique, on voit clairement la contribution essentielle des mathématiciens appliqués et tout particulièrement des spécialistes des méthodes directes en Calcul des Variations dans l'établissement des grandes propriétés qui donnent toute son assise à l'approche variationnelle de la rupture. La contribution des mécaniciens se voit peut-être moins au premier regard, mais pourtant mon sentiment est que les idées de base et les grandes orientations de ce nouveau paradigme leur sont dues pour une part non négligeable. Toujours est-il que le dialogue n'a été possible et fructueux que parce que ces deux communautés parlaient des langues voisines si ce n'est communes. En particulier, il faut souligner le très gros effort de rapprochement qu'a fait l'école italienne de Calcul des Variations en se familiarisant avec les concepts et en acceptant les contraintes imposées par les problématiques des mécaniciens. Symétriquement, les mécaniciens de ma génération ont pu dialoguer parce qu'ils avaient reçu une formation mathématique suffisante. Si l'on veut que la source ne se tarisse pas, il faut donc que les nouvelles générations aient conscience de tout cela et agissent en conséquence. J'espère que ce texte y contribuera.

## Conflit d'intérêt

L'auteur n'a aucun conflit d'intérêt à déclarer.
Références

[1] D. C. Drucker, «A definition of stable inelastic material», ASME J. Appl. Mech. 26 (1959), p. 101-106.

[2] A. A. Ilyushin, «On the postulate of stability», J. Appl. Math. Mech. 25 (1961), p. 746-752.

[3] G. A. Francfort, J.-J. Marigo, « Stable damage evolution in a brittle continuous medium », Eur. J. Mech. A/Solids 12 (1993), n° 2, p. 149-189.

[4] G. Allaire, *Shape Optimization by the Homogenization Method*, Applied Mathematical Sciences, vol. 146, Springer-Verlag, New York, 2002.

[5] A. Griffith, « The phenomena of rupture and flow in solids », Philos. Trans. R. Soc. Lond. A 221 (1921), p. 163-198.

[6] B. Bourdin, G. A. Francfort, J.-J. Marigo, « The variational approach to fracture », J. Elast. 91 (2008), n° 1-3, p. 5-148.

[7] G. A. Francfort, J. J. Marigo, « Revisiting brittle fracture as an energy minimization problem », J. Mech. Phys. Solids 46 (1998), n° 8, p. 1319-1342.

[8] A. Mielke, « Evolution of rate-independent systems », in *Evolutionary Equations*, Handbook of Differential Equations, vol. II, Elsevier/North-Holland, Amsterdam, 2005, p. 461-559.

[9] A. Chambolle, V. Crismale, « Existence of strong solutions to the Dirichlet problem for the Griffith energy », Calc. Var. Partial Differ. Equ. 58 (2019), n° 4, article no. 136.

[10] A. Giacomini, M. Ponsiglione, « A $\Gamma$-convergence approach to stability of unilateral minimality properties in fracture mechanics and applications », Arch. Ration. Mech. Anal. 180 (2006), n° 3, p. 399-447.

[11] L. Ambrosio, V. Tortorelli, « Approximation of functionals depending on jumps by elliptic functional via gamma-convergence », Commun. Pure Appl. Math. 43 (1990), p. 999-1036.

[12] L. Ambrosio, V. Tortorelli, « On the approximation of free discontinuity problems », Bollettino della unione matematica italiana 6-B (1992), n° 1, p. 105-123.

[13] A. Braides, $\Gamma$-convergence for Beginners, Lecture Series in Mathematics and its Applications, vol. 22, Oxford University Press, New York, 2002.

[14] K. Pham, J.-J. Marigo, « Approche variationnelle de l'endommagement : I. Les concepts fondamentaux », C. R. Méc. 338 (2010), n° 4, p. 191-198.

[15] K. Pham, J.-J. Marigo, « Approche variationnelle de l'endommagement : II. Les modèles à gradient », C. R. Méc. 338 (2010), n° 4, p. 199-206.

[16] K. Pham, H. Amor, J.-J. Marigo, C. Maurini, « Gradient damage models and their use to approximate brittle fracture », Int. J. Damage Mech. 20 (2011), n° 4, p. 618-652, Isiweb.

[17] B. Bourdin, G. A. Francfort, J.-J. Marigo, « Numerical experiments in revisited brittle fracture », J. Mech. Phys. Solids 48 (2000), n° 4, p. 797-826.

[18] B. Bourdin, « Numerical implementation of the variational formulation for quasi-static brittle fracture », Interfaces Free Boundaries 9 (2007), p. 411-430.

[19] P. Sicsic, J.-J. Marigo, « From gradient damage laws to Griffith's theory of crack propagation », J. Elast. 113 (2013), n° 1, p. 55-74.

[20] B. Bourdin, J.-J. Marigo, C. Maurini, P. Sicsic, « Morphogenesis and propagation of complex cracks induced by thermal shocks », Phys. Rev. Lett. 112 (2014), n° 1, article no. 014301.

[21] P. Sicsic, J.-J. Marigo, C. Maurini, « Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling », J. Mech. Phys. Solids 63 (2014), p. 256-284.

[22] E. Tanné, T. Li, B. Bourdin, J.-J. Marigo, C. Maurini, « Crack nucleation in variational phase-field models of brittle fracture », J. Mech. Phys. Solids 110 (2018), p. 80-99.

[23] A. Takei, B. Roman, J. Bico, « Forbidden directions for the fracture of thin anisotropic sheets: an analogy with the Wulff plot », Phys. Rev. Lett. 110 (2013), article no. 144301.

[24] A. Ibarra, J. F. Fuentealba, B. Roman, F. Mela, « Predicting tearing paths in thin sheets », Phys. Rev. E 100 (2019), article no. 023002.

[25] G. Allaire, *Conception Optimale de Structures*, Mathématiques & Applications, Springer-Verlag, Berlin Heidelberg, 2007.

[26] J.-J. Marigo, « Modelling of brittle and fatigue damage for elastic material by growth of microvoids », Eng. Fract. Mech. 21 (1985), n° 4, p. 861-874.

[27] S. Andrieux, Y. Bamberger, J.-J. Marigo, « A model of micro-cracked material for concretes and rocks », J. de Méc. Théor. Appl. 5 (1986), n° 3, p. 471-513.

[28] J.-J. Marigo, « Constitutive relations in plasticity, damage and fracture mechanics based on a work property », Nucl. Eng. Des. 114 (1989), p. 249-272.

[29] G. A. Francfort, A. Garroni, « A variational view of partial brittle damage evolution », Arch. Ration. Mech. Anal. 182 (2006), p. 125-152.

[30] J.-F. Babadjian, F. Iurlano, F. Rindler, « Concentration versus oscillation effects in brittle damage », Commun. Pure Appl. Math. 74 (2021), n° 9, p. 1803-1854.
[31] G. Barenblatt, « The mathematical theory of equilibrium cracks in brittle fracture. », *Adv. Appl. Mech.* **7** (1962), p. 55-129.

[32] G. D. Maso, F. Iurlano, « Fracture models as $\Gamma$-limits of damage models », *Commun. Pure Appl. Anal.* **12** (2013), p. 1657-1686.

[33] J.-J. Marigo, « La mécanique de l'endommagement au secours de la mécanique de la rupture : l'évolution de cette idée en un demi-siècle », *Comptes Rendus Mécanique* **351** (2023), n° S1, à paraitre.

[34] D. Bahat, *Tectonofractography*, Springer-Verlag, Berlin Heidelberg, 1991.

[35] A. Yuse, M. Sano, « Transition between crack patterns in quenched glass plates », *Letters to Nature* **362** (1993), p. 329-331.

[36] Y. Shao, Y. Zhang, X. Xu, Z. Zhou, W. Li, B. Liu, « Effect of crack pattern on the residual strength of ceramics after quenching », *J. Am. Ceram. Soc.* **94** (2011), n° 9, p. 2804-2807.

[37] A. Chambolle, G. A. Francfort, J.-J. Marigo, « When and how do cracks propagate? », *J. Mech. Phys. Solids* **57** (2009), n° 9, p. 1614-1622.