r3.06.13 Eléments finis d’interface mixte pour des modèles de zone cohésive (xxx_INTERFACE et xxx_INTERFACE_S)

Résumé:

Ce document présente un élément fini d’interface mixte pour des modèles de zone cohésive. Les mots clés correspondant pour les modélisations dans Code_Aster sont PLAN_INTERFACE, AXIS_INTERFACE,3D_INTERFACE, PLAN_INTERFACE_S, AXIS_INTERFACE_S et 3D_INTERFACE_S.

Table des matières

Construction de l’élément fini d’interface mixte

Principe de minimisation d’énergie en mécanique de la rupture

Dans leur approche de la mécanique de la rupture, Francfort et Marigo 5 décrivent l’état d’une structure \(O\) par un champ de déplacement \(u\) pouvant admettre des discontinuités \(\delta =\left[\left[u\right]\right]\) à travers les surfaces \(\Gamma (u)\) . Pour un chargement donné le déplacement est caractérisé par une condition d’optimalité : il correspond au minimum local de l’énergie potentielle définie comme la somme de l’énergie de déformation élastique, de l’énergie cohésive et du travail des forces extérieures:

\(\begin{array}{}{E}_{\text{pot}}(u)={E}_{\text{el}}(u)+{E}_{\text{fr}}(\left[\left[u\right]\right])-{W}_{\text{ext}}(u)\\ {E}_{\text{el}}(u)=\underset{\Omega /\Gamma (u)}{\int}\Phi (\epsilon (u))d\Omega ;{E}_{\text{fr}}(\delta )=\underset{\Gamma (u)}{\int}\Pi (\delta )\text{dS}\end{array}\)

où \(\epsilon\) désigne le tenseur des déformations, \(\Phi\) la densité d’énergie de déformation (volumique) et \(\Pi\) la densité d’énergie cohésive (surfacique). Dans cette formulation, les mécanismes de rupture sont considérés comme réversibles. L’irréversibilité sera introduite dans la partie suivante, à partir de variables internes. Les caractéristiques du problème sont déduites de la minimalisation, à savoir:

• la prévision du trajet de fissuration, puisque l’on prend en compte tous les déplacements discontinus possibles;

• un critère en contrainte pour l’amorçage de la fissure (voir 5 ).

• la loi cohésive relie la discontinuité du déplacement \(\delta\) au vecteur de traction \(t\) à partir de la dérivée (généralisée) de \(\Pi\) ;

• les conditions de contact, sont gérées via l’énergie cohésive à laquelle on ajoute une fonction d’indicatrice interdisant l’interpénétration des lèvres.

Cependant, tenir compte de tous les déplacements discontinus possibles entraîne des difficultés numériques liées à la discrétisation de l’espace fonctionnel \(\text{BD}(O)\) . En particulier, autoriser des discontinuités à travers chacun des éléments finis conduit à une dépendance des résultats au maillage. Des travaux alternatifs sont réalisés, basés sur la régularisation des discontinuités 5 ou sur l’adaptation du maillage 5 . Ces derniers sont cependant limités à des énergies de surface de type Griffith (sans forces cohésives) et conduisent à des difficultés numériques importantes. C’est pourquoi nous introduisons une simplification dans notre modèle en considérant que les surfaces discontinuités de déplacement potentielles \(\Gamma\) sont postulées a priori et qu’elles ne dépendent plus de \(u\) .

Recherche de point selle par une méthode de décomposition – coordination

Malgré le postulat sur la direction de fissuration, la minimisation de l’énergie n’est pas évidente compte tenu de la non-dérivabilité de l’énergie de surface. On se base ici sur la technique de décomposition – coordination introduite dans 5 , qui condense la non-dérivabilité à un niveau local (points de Gauss).

Lagrangien augmenté

La relation entre le champ de discontinuité \(\delta\) et le champ de déplacement \(u\) est introduite explicitement dans la minimisation, l’énergie totale E dépend donc explicitement de \(u\) et de \(\delta\) :

\(E(u,\delta )=\underset{\Omega \setminus \Gamma }{\int}\Phi (\epsilon (u))d\Omega -{W}_{\mathrm{ext}}(u)+\underset{\Gamma}{\int}\Pi (\delta )\mathrm{d\Gamma }\)

La minimisation de l’énergie potentielle est alors équivalente au problème de minimisation sous contrainte suivant (le déplacement appartenant à l’ensemble des déplacements cinématiquement admissibles):

\(\underset{\begin{array}{}u,\delta \\ \left[\left[u\right]\right]=\delta \end{array}}{\min}E(u,\delta )\)

La contrainte linéaire \(\left[\left[u\right]\right]=\delta\) est traitée par dualisation: une solution \((u,\delta )\) de correspond à un point selle \((u,\delta ,\lambda )\) du Lagrangien suivant, où \(\lambda\) correspond au champ des multiplicateurs de Lagrange:

\(L(u,\delta ,\lambda )\underset{\mathrm{déf}}{=}E(u,\delta )+\underset{\Gamma}{\int}\lambda \cdot (\left[\left[u\right]\right]-\delta )\mathrm{d\Gamma }\)

Enfin, pour gagner en coercivité par rapport à \(\delta\) , ce qui se révèlera être nécessaire par la suite, on ajoute un terme de pénalisation quadratique n’ayant aucune influence sur la solution puisqu’il est égal à zéro lorsque la contrainte en vérifiée. Le Lagrangien augmenté \({L}_{r}\) , avec le coefficient de pénalisation \(r\) s’écrit alors:

\({L}_{r}(u,\delta ,\lambda )\underset{\mathrm{déf.}}{=}E(u,\delta )+\underset{\Gamma}{\int}\lambda \cdot (\left[\left[u\right]\right]-\delta )\mathrm{d\Gamma }+\underset{\Gamma}{\int}\frac{r}{2}{(\left[\left[u\right]\right]-\delta )}^{2}\mathrm{d\Gamma }\)

Caractérisation du point selle

Les conditions d’optimalité d’ordre 1 pour le Lagrangien augmenté , s’écriventde la façon suivante :

\(\forall \mathrm{\delta \delta }\underset{\Gamma}{\int}\left[t-\lambda -r(\left[\left[u\right]\right]-\delta )\right]\cdot \mathrm{\delta \delta }=0\text{avec}t\in \partial \Pi (\delta )\)

\(\forall \mathrm{\delta u}\underset{\Omega \setminus \Gamma }{\int}\sigma :\epsilon (\mathrm{\delta u})+\underset{\Gamma}{\int}\left[\lambda +r(\left[\left[u\right]\right]-\delta )\right]\cdot \left[\left[\mathrm{\delta u}\right]\right]={W}_{\text{ext}}(\mathrm{\delta u})\text{avec}\sigma =\frac{\partial \Phi }{\partial \epsilon }(\epsilon )\)

\(\forall \mathrm{\delta \lambda }\underset{\Gamma}{\int}(\left[\left[u\right]\right]-\delta )\cdot \mathrm{\delta \lambda }=0\)

La signification du sous-différentiel \(\partial \Pi\) est donnée dans la partie 3 . A ce stade, on se contente de dire que \(t\in \partial \Pi (\delta )\) signifie que \(t\) et \(\delta\) sont reliés par la loi cohésive. L’équation donne donc une interprétation du multiplicateur de Lagrange \(\lambda\) : hormis le terme de pénalisation, il s’agit des forces cohésives. L’équation exprime l’équilibre des efforts dans le volume et le long de la surface de la discontinuité \(\Gamma\) . Enfin, impose la contrainte entre le champ de déplacement et ses discontinuités.

Discrétisation de l’élément fini

Puisque la trajectoire de la fissure \(\Gamma\) est définie a priori et grâce à l’hypothèse des petites perturbations, la discrétisation spatiale du système - peut s’appuyer sur un élément fini d’interface. Supposons que les sous-domaines \({O}_{-}\) et \({O}_{+}\) (parties de \(O\) respectivement en dessous et au dessus de \(\Gamma\) ) sont discrétisés par des tétraèdres ou des hexaèdres de sorte que les nœuds de chaque côté de \(\Gamma\) coïncident [1] . Dans ce cas, les prismes ou hexaèdres dégénérés peuvent être utilisés pour discrétiser \(\Gamma\) et relier les deux lèvres \({\Gamma}_{-}\) et \({\Gamma}_{+}\) de la fissure cohésive potentielle (voir Figure 2.3-a ).

Figure 2.3-a : Discrétisation par élément fini d’interface.

Soit un maillage caractérisé par le paramètre \(h\) (par exemple, la taille maximale des éléments finis). Dans le domaine \(O\) on adopte une interpolation quadratique avec des éléments finis de Lagrange classiques (P2-continu). L’espace discrétisé des champs de déplacements \({U}_{h}\) s’écrit:

\({U}_{h}=\left\lbrace u;\forall x\in \Omega u(x)=\left[N(x)\right]\left\lbrace U\right\rbrace \right\rbrace\)

où \(\left\lbrace U\right\rbrace\) désigne le vecteur des déplacements nodaux et \(\left[N\right]\) la matrice des fonctions de formes quadratiques. La trace du déplacement interpolé sur \({\Gamma}_{-}\) et \({\Gamma}_{+}\) est également quadratique par morceaux et la discontinuité du déplacement s’exprime de la manière suivante:

\(\forall x\in \Gamma {u}_{\mid {\Gamma}_{-}}(x)=\left[{N}_{-}(x)\right]\left\lbrace U\right\rbrace ;{u}_{\mid {\Gamma}_{+}}(x)=\left[{N}_{+}(x)\right]\left\lbrace U\right\rbrace\)

\(\forall x\in \Gamma \left[\left[u(x)\right]\right]=\left[D(x)\right]\left\lbrace U\right\rbrace \text{avec}\left[D(x)\right]\underset{\text{déf.}}{=}\left[{N}_{+}(x)\right]-\left[{N}_{-}(x)\right]\)

où \(\left[{N}_{-}\right]\) et \(\left[{N}_{+}\right]\) correspondent respectivement à la trace de \(\left[N\right]\) sur \({\Gamma}_{-}\) et \({\Gamma}_{+}\) et où \(\left[D\right]\) désigne la matrice des fonctions des formes quadratiques qui interpole le saut de déplacement. Notons qu’il peut être utile d’introduire une rotation dans \(\left[D\right]\) afin d’obtenir les composantes de \(\left[\left[u\right]\right]\) dans un système de coordonnées locales et distinguer ainsi les parties normales et tangentielles.

Le champ des multiplicateurs de Lagrange \(\lambda\) est interpolé sur \(\Gamma\) par des fonctions de formes linéaires par morceaux (P1–continues), conduisant à l’espace discrétisés des multiplicateurs de Lagrange \({L}_{h}\) suivant :

\({L}_{h}=\left\lbrace \lambda ;\forall x\in \Gamma \lambda (x)=\left[L(x)\right]\left\lbrace \Lambda \right\rbrace \right\rbrace\)

où \(\left\lbrace \Lambda \right\rbrace\) correspond aux inconnues nodales du multiplicateur de Lagrange et \(\left[L\right]\) désigne la matrice des fonctions de formes linéaires sur \(\Gamma\) . De cette façon, la contrainte \(\left[\left[u\right]\right]=\delta\) imposée par n’est réalisée qu’au sens faible.

Pour finir, la discrétisation du champ de discontinuité \(\delta \in {D}_{h}\) est basée sur les points de collocation sur \(\Gamma\) , de coordonnées \({x}_{g}\) . Ces points correspondent aux points de Gauss, de poids \({\omega}_{g}\) , utilisés pour calculer les intégrales. On distingue ici une modélisation standard, et une modélisation dite sous-intégrée.

• Dans la modélisation standard, les points de Gauss sont au nombre de 3 (segments), 6 (triangles) ou 9 (quadrilatères) par élément. Ceci correspond à une interpolation P2-discontinu de \(\delta\) . Ce choix présente l’avantage de limiter le risque d’apparition de pivots nuls dans la matrice tangente, mais il ne vérifie pas la condition LBB dans la limite d’un paramètre de pénalisation \(r\) infini.

• Dans la modélisation sous-intégrée, les points de Gauss sont au nombre de 2 (segments), 3 (triangles) ou 4 (quadrilatères) par élément. Ceci correspond à une interpolation P1-discontinu de \(\delta\) . Ce choix vérifie mieux la condition LBB dans la limite \(r\to \infty\) , mais il est susceptible de faire apparaître des pivots nuls dans la matrice tangente. Ces pivots nuls apparaissent lorsque la rigidité des éléments au contact de l’interface ne suffit pas à assurer la coercivité de la formulation, en particulier lorsque l’interface se trouve sur un plan de symétrie, ou entre un volume et un élément de structure (barre, grille ou membrane).

La discrétisation spatiale, les notations ainsi que les schémas des éléments finis pour la modélisation sous-intégrée sont résumés sur la Figure 2.3-b .

Figure 2.3-b : Discrétisation spatiale

Le champ \(\delta\) disparaît de la formulation globale grâce à la condensation statique. En effet, avec la discrétisation adoptée, la résolution de revient à satisfaire la loi cohésive en chaque point de collocation:

\({t}_{g}={\lambda}_{g}+r(\left[\left[{u}_{g}\right]\right]-{\delta}_{g})\in \partial \Pi ({\delta}_{g})\text{avec}\lbrace \begin{array}{cc}\left[\left[{u}_{g}\right]\right]=\left[{D}_{g}\right]\left\lbrace U\right\rbrace & ;\left[{D}_{g}\right]=\left[D({x}_{g})\right]\\ {\mu}_{g}=\left[{L}_{g}\right]\left\lbrace \Lambda \right\rbrace & ;\left[{L}_{g}\right]=\left[L({x}_{g})\right]\end{array}\)

L’intégration de la relation constitutive (i.e. la résolution de , cf. partie suivante), permet de calculer \({\delta}_{g}\) fonction de \(\left\lbrace U\right\rbrace\) et \(\left\lbrace \Lambda \right\rbrace\) que l’on note \(\delta\) :

\({t}_{g}={\lambda}_{g}+r(\left[\left[{u}_{g}\right]\right]-{\delta}_{g})\in \partial \Pi ({\delta}_{g})\iff {\delta}_{g}=\stackrel{ˆ}{\delta}(\left[\left[{u}_{g}\right]\right],{\lambda}_{g})=\delta (U,\Lambda )\)

Le paramètre de pénalisation \(r\) permet d’assurer l’unicité de quels que soient \(\left\lbrace U\right\rbrace\) et \(\left\lbrace \Lambda \right\rbrace\) (cf. partie suivante). Cela constitue une condition nécessaire pour garantir de la robustesse du modèle.

On a alors le système non linéaire dont les inconnues sont les déplacements nodaux \(\left\lbrace U\right\rbrace\) et les multiplicateur de Lagrange nodaux \(\left\lbrace \Lambda \right\rbrace\) :

\(\underset{\Omega \setminus \Gamma }{\int}{\left[\nabla N\right]}^{T}\text{:σ}(U)+\sum_{g}{\omega}_{g}{\left[{D}_{g}\right]}^{T}\cdot (\left[{L}_{g}\right]\left\lbrace \Lambda \right\rbrace +r\left[{D}_{g}\right]\left\lbrace U\right\rbrace -r\delta (U,\Lambda ))=\left\lbrace {F}_{\text{ext}}\right\rbrace\)

\(\sum_{g}{\omega}_{g}{\left[{L}_{g}\right]}^{T}\cdot (\left[{D}_{g}\right]\left\lbrace U\right\rbrace -\delta (U,\Lambda ))=0\)

L’intégrale de volume et le vecteur nodal des efforts extérieurs sont calculés de façon habituelle. Enfin, le système est résolu simultanément par rapport à \(\left\lbrace U\right\rbrace\) et \(\left\lbrace \Lambda \right\rbrace\) au moyen d’une méthode de Newton (généralisée), où l’opérateur tangent est symétrique puisque - correspond à une recherche de point selle.

[1]

Cela ne pose a priori pas de problème aux algorithmes de maillage, au moins pour les fissures ayant une forme simple, par exemple plane.

Intégration de la loi cohésiveCZM_TAC_MIX

Le comportement cohésif est déterminé par la densité de l’énergie cohésive \(\Pi (\delta )\) . Bien que la formulation variationnelle présentée dans la partie précédente soit indépendante du choix de l’énergie cohésive, nous nous intéressons à présent à un modèle cohésif particulier dont on détaille les spécificités. Ces principales caractéristiques sont:

• conditions de contact (non interpénétration des lèvres) ;

• adhérence initiale parfaite (pas de régularisation de l’énergie de surface);

• couplage entre les modes de rupture (traction et cisaillement) ;

• rupture totale, (forces cohésives nulles au-delà d’un certain seuil d’endommagement) ;

• irréversibilité de la rupture.

Cela signifie en particulier qu’il n’y a aucun frottement final, ni de distinction entre les mécanismes de rupture en traction et en cisaillement.

**Remarque:*

Le mot clé Code_Aster correspondant à cette loi cohésive est CZM_TAC_MIX. Pour plus de détails sur cette loi (variables internes, paramètres d’entrée, matrice tangente cohérente) et sur d’autres lois cohésives, on peut se référer à 5 (R7.02.11).

Notations préliminaires

En raison de la contrainte de non interpénétration des lèvres de la fissure, la direction \(n\) normale pour la surface de discontinuité potentielle \(\Gamma\) (ouverture/compression) doit être distinguée des directions dans le plan (glissement, cisaillement). Pour ce faire, les notations suivantes sont introduites, où \(v\) représente toute quantité vectorielle:

\({v}_{n}=v\cdot n;{v}_{\text{//}}=v-{v}_{n}n;{\langle v\rangle }_{+}=\langle {v}_{n}\rangle n+{v}_{\text{//}};{\parallel v\parallel }_{+}={({\langle v\rangle }_{+}\cdot {\langle v\rangle }_{+})}^{1/2}\)

et où \(\langle .\rangle\) désigne la partie positive d’un scalaire.

Définition de la loi cohésive

La loi cohésive proposée par Talon et Curnier 5 répond aux cinq points mentionnés précédemment et rentre dans le cadre de la formulation énergétique . Les réponses du modèle soumis à une traction pure et à un cisaillement pur, ainsi que le critère d’amorçage en contrainte sont détaillés Figure 3.2-a .

Figure 3.2-a : Réponses de la loi Talon Curnier CZM_TAC_MIX en traction-compression pure (a), cisaillement pur (b) et critère d’amorçage en contrainte (c)

L’énergie cohésive correspondante est définie comme suit:

Tout d’abord, les discontinuités de déplacement en ouverture et en cisaillement sont rassemblées dans un scalaire unique \({\delta}_{\text{eq}}\) qui mesure l’amplitude de la discontinuité, répondant ainsi au point 3 :

\({\delta}_{\text{eq}}=N(\delta )=\sqrt{\delta \cdot \delta }\)

Ensuite, on tient compte de l’irréversibilité 3 via une nouvelle variable interne scalaire \(\kappa\) qui mesure le niveau de chargement maximum courant:

\(\kappa (t)=\underset{{t}^{'}<t}{\max}{\delta}_{\text{eq}}({t}^{'})\)

L’énergie cohésive dépend de la variable interne \(\kappa\) et du saut de déplacement équivalent \({\delta}_{\text{eq}}\) . Les conditions de contact 3 sont gérées par une fonction indicatrice qui impose une discontinuité normale positive \({\delta}_{n}\ge 0\) :

\(\Pi (\delta ,\kappa )={I}_{{ℝ}^{+}}({\delta}_{n})+\psi (\max({\delta}_{\text{eq}},\kappa ))\)

La fonction \(\psi\) caractérise notamment la réaction à une sollicitation en mode I pur. D’après 5 , \(\psi\) doit être une fonction croissante et dérivable, où \({\sigma}^{c}={\psi}^{'}(0)\) définit la contrainte critique, paramètre du critère d’amorçage 3 représenté sur la Figure 3.2-a . Par ailleurs, la stabilité du processus de rupture requiert que \(\psi\) soit concave 5 . Enfin, la rupture finale 3 se produit dès que \(\psi\) atteint sa limite supérieure \({G}^{c}\) pour une valeur finie \({\delta}_{\text{eq}}={\delta}^{c}\) , où \({G}^{c}\) est l’énergie de rupture et \({\delta}^{c}\) l’ouverture critique au-delà de laquelle les forces cohésives s’annulent.

Ainsi, on choisit la fonction \(\psi\) suivante qui correspond aux réponses représentées sur la Figure 3.2-a :

\(\psi ({\delta}_{\text{eq}})=\lbrace \begin{array}{cc}{G}^{c}\frac{{\delta}_{\text{eq}}}{{\delta}^{c}}(2-\frac{{\delta}_{\text{eq}}}{{\delta}^{c}})& \text{if}{\delta}_{\text{eq}}\le {\delta}^{c}\\ {G}^{c}& \text{if}{\delta}_{\text{eq}}\ge {\delta}^{c}\end{array}\)

avec la relation entre les paramètres matériaux;

\({G}^{c}=\frac{1}{2}{\sigma}^{c}{\delta}^{c}\)

A ce stade, l’énergie et la loi cohésive qui en dérive sont entièrement définies. Cependant, on peut exprimer explicitement la relation entre la discontinuité de déplacement \(\delta\) et le vecteur contrainte \(t\) , de manière condensée dans l’inclusion différentielle suivante, où \(\partial \Pi\) est le sous-différentiel de \(\Pi\) (voir Clarke 5 ):

\(t\in \partial \Pi (\delta )\)

Pour une valeur donnée de \(\kappa\) , on peut interpréter le sous-différentiel \(\partial \Pi (\delta )\) comme le cône formé par l’ensemble des pentes des dérivées directionnelles de \(\Pi\) en \(\delta\) pour toutes les directions admissibles. Mathématiquement, cela se formule de la façon suivante:

\(\partial \Pi (\delta )=\left\lbrace t\in {R}^{3};\forall \upsilon \in {R}^{3}t\cdot \upsilon \le {\Pi}^{°}(\delta ,\upsilon )\right\rbrace\)

où \({\Pi}^{°}(\delta ,\upsilon )\) est la dérivée directionnelle de \(\Pi\) par rapport à \(\delta\) dans de la direction \(\upsilon\) :

\({\Pi}^{°}(\delta ,\upsilon )=\underset{\begin{array}{}\rho \to {0}^{+}\\ d\to \delta \end{array}}{\text{limsup}}\frac{\Pi (d+\rho \upsilon )-\Pi (d)}{\rho}\)

Cette définition coïncide avec le gradient de \(\Pi\) partout où celle-ci est dérivable. D’après la définition , les points méritant une attention particulière sont \(\delta =0\) , \({\delta}_{n}=0\) et \({\delta}_{\text{eq}}=\kappa\) . On déduit donc de les quatre cas de figure suivants (illustrés en couleur sur la Figure 3.2-a ):

1. Point (\(\kappa =0\) et \(\delta =0\) ): adhérence parfaite, i.e. critère d’amorçage (branche rose)

\(\partial \Pi (\delta )=\left\lbrace t\in {R}^{3};{\parallel t\parallel }_{+}\le {\sigma}^{c}\right\rbrace\)

2. Domaine où \({\delta}_{\text{eq}}<\kappa\) : retour à zéro à contrainte nulle (branche verte)

\(\partial \Pi (\delta )=\left\lbrace {t}_{n}n;{t}_{n}\le 0\text{et}{\delta}_{n}\ge 0\text{et}{t}_{n}{\delta}_{n}=0\right\rbrace\)

3. Hyper cône \({\delta}_{\text{eq}}=\kappa >0\) : décharge verticale de la contrainte (branche orange)

\(\partial \Pi (\delta )=\left\lbrace {t}_{n}n+\rho \delta ;0\le \rho \le \frac{{\psi}^{'}(\kappa )}{\kappa}\text{et}{t}_{n}\le 0\text{et}{\delta}_{n}\ge 0\text{et}{t}_{n}{\delta}_{n}=0\right\rbrace\)

4. Domaine où \({\delta}_{\text{eq}}>\kappa\) : endommagement (branche bleue)

\(\partial \Pi (\delta )=\left\lbrace {t}_{n}n+{\psi}^{'}({\delta}_{\text{eq}})\frac{\delta}{{\delta}_{\text{eq}}};{t}_{n}\le 0\text{et}{\delta}_{n}\ge 0\text{et}{t}_{n}{\delta}_{n}=0\right\rbrace\)

Remarques:

• Il existe un domaine pour la force cohésive, où la discontinuité est nulle : il s’agit du critère d’amorçage représenté Figure 3.2-a lié à la non dérivabilité de \(\Pi\) en \(\delta =0\) *. La forme du domaine dépend de l’expression de* \({\delta}_{\text{eq}}\) *.*

• La condition de Kuhn et Tucker, qui apparaît dans - permet de décrire les conditions de contact (branche rouge Figure 3.2-a ). Dans ce cas, \({t}_{n}\) est négative, et correspond à une compression dont la valeur n’est pas définie par la loi cohésive.

• Il existe un saut entre le retour à zéro à contrainte nulle (branche verte) et l’endommagement (branche bleue), la force cohésive n’est donc pas continue par rapport au saut de déplacement .

• Dans le cas d’un mode I pur ou d’un mode de cisaillement pur, les réponses de la loi cohésive sont celles représentées sur la Figure 3.2-a . Les valeurs maximales en traction et en cisaillement sont égales en raison du choix de la norme .

Intégration numérique

D’après , l’intégration numérique de la loi cohésive revient à calculer \({\delta}_{g}\) pour des valeurs données de \(\left[\left[{u}_{g}\right]\right]\) et de \({\lambda}_{g}\) , et ce pour chaque point de Gauss (on omet à présent l’indice \(g\) ). En fait, est une caractérisation du minimum suivant:

\(\underset{\delta \in {ℝ}^{3}}{\min}\left[\lambda \cdot (\left[\left[u\right]\right]-\delta )+\frac{r}{2}{(\left[\left[u\right]\right]-\delta )}^{2}+\Pi (\delta ,\kappa )\right]\iff \lambda +r\left[\left[u\right]\right]-r\delta \in \partial \Pi (\delta ,\kappa )\)

De plus, l’évolution de la variable interne \(\kappa\) (voir ) doit être prise en compte. Une discrétisation temporelle est alors nécessaire, considérons une série d’instants \({t}^{0}<{t}^{1}<\dots <{t}^{n}\) et les quantités correspondantes \({\lambda}^{n}\) , \(\left[\left[{u}^{n}\right]\right]\) , \({\delta}^{n}\) et \({\kappa}^{n}\) . Lors de la résolution, le processus itératif est divisé en deux parties:

\({\delta}^{n}=\underset{\delta \in {ℝ}^{3}}{\text{argmin}}\left[{\lambda}^{n}\cdot (\left[\left[{u}^{n}\right]\right]-\delta )+\frac{r}{2}{(\left[\left[{u}^{n}\right]\right]-\delta )}^{2}+\Pi (\delta ,{\kappa}^{n-1})\right]\)

\({\kappa}^{n}=\max({\kappa}^{n-1},{\delta}_{\text{eq}}^{n})\)

La seconde partie est triviale, la minimisation revient donc à résoudre la première avec \(\kappa ={\kappa}^{n-1}\) paramètre connu. Une interprétation graphique de l’inclusion différentielle est fournie sur la Figure 3.3-a ( mode I pur sans décharge): la solution correspond à l’intersection de la fonction linéaire \(\delta \to \lambda +r\left[\left[u\right]\right]-r\delta\) (avec une pente négative donnée par le coefficient de pénalité \(r\) ) avec le graphe \(\partial \Pi (\delta ,{\kappa}^{n-1})\) .

Figure 3.3-a : Solution de l’intégration du comportement.

Pour que l’intégration soit robuste, il est nécessaire que la fonction entre crochets dans soit strictement convexe par rapport à \(\delta\) (i.e. que le minimum soit unique). En introduisant \({H}_{c}\) , le module adoucissant de la loi cohésive, on exhibe une condition suffisante pour satisfaire cette condition (voir 5 pour plus de détails) :

\(r>\underset{x\ge 0}{\max}\mid {\psi}^{''}(x)\mid =\frac{{\sigma}^{c}}{{\delta}^{c}}={H}^{c}\)

**Remarque:*

\(r=100\times {H}^{c}\) est une valeur recommandée à l’utilisateur de Code_Aster .

On suppose à présent que la condition est remplie, de plus la fonction entre crochets dans est semi-continue inférieurement, ce qui garantit l’existence et l’unicité du minimum.

Détaillons maintenant l’intégration numérique de la loi, c’est-à-dire la résolution de en fonction de \(\tau =\lambda +r\left[\left[u\right]\right]\) . Compte tenu de l’existence et de l’unicité de la solution, il est intéressant d’utiliser l’inclusion différentielle et l’écriture du sous-différentiel fournie par à . On retrouve ainsi les quatre cas de figure:

• \(\text{si}\kappa =0\text{et}{\parallel \tau \parallel }_{+}\le {\sigma}^{c}={\psi}^{'}(0)\) : adhérence parfaite

\(\delta =0\)

• \(\text{si}\kappa >0\text{et}{\parallel \tau \parallel }_{+}<\mathrm{r\kappa }\) : retour à zéro à contrainte nulle

\(\delta =\frac{{\langle \tau \rangle }_{+}}{r}\)

• \(\text{si}\kappa >0\text{et}\mathrm{r\kappa }\le {\parallel \tau \parallel }_{+}\le \mathrm{r\kappa }+{\psi}^{'}(\kappa )\) : décharge verticale

\(\delta =\kappa \frac{{\langle \tau \rangle }_{+}}{{\parallel \tau \parallel }_{+}}\)

• \(\text{si}\mathrm{r\kappa }+{\psi}^{'}(\kappa )<{\parallel \tau \parallel }_{+}\) : endommagement

\(\delta ={\delta}_{\text{eq}}\frac{{\langle \tau \rangle }_{+}}{{\parallel \tau \parallel }_{+}};{\delta}_{\text{eq}}\text{solution de}{\psi}^{'}({\delta}_{\text{eq}})+{\mathrm{r\delta }}_{\text{eq}}={\parallel \tau \parallel }_{+}\)

**Remarque:*

la distinction \(\kappa =0\) ou pas n’est pas nécessaire puisque apparaît comme un cas particulier de .

Pour conclure, les différents régimes de la loi cohésive correspondent aux solutions \(\delta\) fournies par -. Ces dernières sont fournies en fonction de \(\tau =\lambda +r\left[\left[u\right]\right]\) . Aucune méthode numérique n’est nécessaire, en effet, la seule équation à résoudre (dans ) est linéaire par morceaux.

Conclusion

Un modèle élément fini d’interface, est proposé afin de modéliser l’évolution de fissures cohésives le long d’un trajet prédéfini. Ce dernier est compatible avec l’utilisation d’éléments finis volumiques classiques. Ses inconnues sont les déplacements nodaux sur les lèvres de la fissure cohésive ainsi que les multiplicateurs de Lagrange nodaux, correspondant à la densité surfacique des forces cohésives. Le Lagrangien du problème est augmenté, cela permet d’assurer une condition de convexité (via le paramètre de pénalisation) qui garantit l’unicité de la solution lors de l’intégration locale de la loi.

Cette approche présente toutefois les limites et les inconvénients suivants:

• Les trajectoires potentielles de la fissure doivent être postulées a priori.

• Des degrés de liberté supplémentaires, correspondent aux forces cohésives, sont introduits. Cependant leur nombre reste peu élevé puisqu’ils sont limités aux trajectoires potentielles de la fissure.

• L’introduction de multiplicateurs de Lagrange conduit à une formulation mixte: la résolution du problème revient donc à une recherche de point selle et plus à celle d’un minimum comme c’était le cas dans la formulation énergétique initiale 5 .

• Il est nécessaire d’augmenter le Lagrangien pour avoir une propriété de convexité locale. Cela implique l’introduction d’un paramètre de pénalisation, sans influence sur le problème continu, mais qui pourrait affecter les résultats du problème discrétisé. Cependant, les exemples numériques dans 5 montrent que cette sensibilité reste faible et disparaît avec le raffinement du maillage.

• L’intégration locale de la loi cohésive repose sur le calcul des discontinuités de déplacement à partir des forces cohésives. Une démarche inverse est en général adoptée dans la littérature pour les éléments d’interface.

Cependant un certain nombre de propriétés intéressantes se dégagent :

• Aucune régularisation de la loi cohésive n’est nécessaire, en particulier en ce qui concerne l’adhérence initiale ou la condition de contact.

• Le choix d’une discrétisation quadratique pour les déplacements et linéaire pour les multiplicateurs de Lagrange permet de satisfaire la condition LBB. Cette dernière permet d’assurer la convergence de la solution avec le raffinement du maillage en termes de déplacements et de forces cohésives (voir exemple numérique dans 5 ).

• Le taux de convergence que l’on obtiendrait sans interface n’est pas perturbé par la présence des éléments d’interface (voir 5 ).

• La recherche de point selle conduit à une matrice tangente symétrique.

• Ce modèle d’interface est compatible avec les algorithmes de résolution habituels de Code_Aster tels que la méthode Newton, la recherche linéaire ou le pilotage du chargement. Ceci est illustré par des calculs 2D et 3D dans 5 , démontrant ainsi la robustesse et la fiabilité d’un tel modèle.

Bibliographie

1. Bourdin B., Francfort G.A., Marigo J.-J. (2000) Numerical experiments in revisitted brittle fracture. J. Mech. Phys. Solids 48 (4), 797-826.

2. Charlotte M., Laverne J., Marigo J.-J. (2006) Initiation of cracks with cohesive force models: a variational approach. European Journal of Mechanics A/Solids 25 , 649-669.

3. Clarke F.H. Optimization and non smooth analysis. Classics in applied mathematics 5 , SIAM, 1990.

4. Del Piero G., Truskinovsky L. (2001) Macro- and micro- cracking in one dimensional elasticity. Int. J. Solids Struct. 38 , 1135-1148.

5. Fortin M., Glowinski R. Augmented lagrangian methods : application to the numerical solution of boundary-value problems. Studies in mathematics and its applications 15 , North-Holland, 1983.

6. Francfort G., Marigo J.-J. (1998) Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids 46 (8), 1319-1342.

7. Fraternali F. (2007) Free discontinuity finite element models in two-dimensions for in-plane crack problems. Theoretical and Applied Fracture Mechanics 47 , 274-282.

8. Laverne J. Manuel de référence du Code_Aster R7.02.11 : Lois de comportement cohésives CZM_xxx_xxx et pilotage du chargement.

9. Lorentz E., A mixed interface finite element for cohesive zone models, Comput. Methods Appl. Mech. Engrg. 198 (2008), 302-317.

10. Moes N., Dolbow J., Belytschko T. (1999) A finite element method for crack growth without remeshing. Int. J. Num. Meth. Engng. 46 , 131-150.

11. Talon C., Curnier A. (2003) A model of adhesion coupled to contact and friction. Eur. J. Mech. A/Solids 22 (4), 545-565.

Description des versions du document

Indice doc	Version Aster	Auteur(s) ou contributeur(s), organisme	Description des modifications
A	9.4	J.Laverne EDF/R&D/AMA E.Lorentz EDF/R&D/SINETICS	Texte initial
	11.3	David EDF/R&D/AMA	Modification concernant le nombre de points utilisés pour la discrétisation du saut