r7.02.19 Eléments cohésifs avec X-FEM

Résumé:

Ce document présente les différents éléments cohésifs disponibles avec la méthode des éléments finis étendus (X-FEM). Trois types de lois cohésives sont disponibles. La loi introduite en premier est une loi régularisée, disponible pour les éléments linéaires, fonction du saut de déplacement uniquement. Suite aux limites de cette dernière, deux lois mixtes ont été introduites: l’une pour les éléments linéaires, l’autre pour les éléments quadratiques. Une procédure de propagation de fissures avec éléments cohésifs est implémentée. Elle se base sur les éléments mixtes linéaires.

Les éléments cohésifs sont implémentés dans Code_Aster en 2D et 3D. Ils sont utilisables avec deux types de discontinuités: les interfaces débouchantes (mot-clé TYPE_DISCONTINUITE=”INTERFACE” dans DEFI_FISS_XFEM) et les fissures de type cohésives (mot-clé TYPE_DISCONTINUITE=”COHESIF”) que l’on utilise pour les études de propagation sur trajet inconnu. La loi cohésive est définie dans STAT_NON_LINE, la résolution s’effectue avec la commande STAT_NON_LINE [U4.51.03].Pour une étude de propagation sur trajet inconnu, l’actualisation de la surface fissurée se fait par la commande PROPA_FISS, et le critère directionnel est calculé par la commande CALC_G.

Table des Matières

Formulation forte du problème cohésif

Formulation des équations du problème général

Notons \(\Omega \subset {ℝ}^{d},d\in \lbrace 1,2\rbrace\) le domaine de calcul non fissuré. Sa frontière \(\partial \Omega\) se décompose en deux parties \({\Gamma}_{u}\) et \({\Gamma}_{g}\) , sur lesquelles sont imposées respectivement des déplacements et une densité d’efforts surfaciques \(t\) (voir fig.). Nous notons \(\Gamma\) la zone potentielle de fissuration sur laquelle est définie la loi cohésive, composée d’une lèvre de fissure supérieure \({\Gamma}^{+}\) et d’une lèvre inférieure \({\Gamma}^{-}\) . Si \({u}^{i}\) est le déplacement sur \({\Gamma}^{i}\) , alors on introduit le saut de déplacement comme \(⟦u⟧\left(x\right)={u}^{+}\left(x\right)-{u}^{-}\left(x\right)\) . Notons \(n\) le vecteur normal sortant de \({\Gamma}^{-}\) , \({n}^{+}\) le vecteur normal sortant de \({\Gamma}^{+}\) , et \(-{t}_{c}^{+}={t}_{c}^{-}={t}_{c}\) la force [1] cohésive qui s’applique sur \({\Gamma}^{-}\) (voir fig.). Sur \(\Gamma\) , nous avons donc \(-\sigma \cdot {n}^{+}=\sigma \cdot n={t}_{c}\) . Avec ces conventions, nous avons \(\mathrm{〚}{u}_{n}\mathrm{〛}=\mathrm{〚}u\mathrm{〛}\cdot n\) positif en ouverture et négatif en interpénétration (voir fig.).

Loi de comportement	\(\sigma =C:\epsilon \text{dans}\Omega\)
Équilibre	\(\nabla \cdot \sigma =f\text{dans}\Omega\)
Efforts surfaciques imposés	\(\sigma \cdot {n}_{\mathit{ext}}=t\text{sur}{\Gamma}_{t}\)
Densité des efforts de cohésion	\(\sigma \cdot n={t}_{c}\text{sur}\Gamma\)
Déplacements imposés	\(u=0\text{sur}{\Gamma}_{u}\)

Tableau 2.1-1 : Équations du problème général.

Figure 2.1-1: Formulation forte du problème cohésif.

Lois cohésives

Lois cohésives régularisées

Le premier type de lois que nous pouvons envisager est une loi dans laquelle l’adhérence initiale n’est pas parfaite: la pente initiale est finie. Deux lois cohésives de cette sorte sont disponibles dans Code_Aster, les lois CZM_EXP_REG etCZM_LIN_REG dont les caractéristiques sont détaillées dans [R7.02.11]. Nous en rappelons ici les principaux points, le lecteur pouvant s’y référer pour plus de détails.Nous détaillons ici l’extension à XFEM pour la loiCZM_EXP_REG,en se basant sur l’article de référence [bib 10 ]. L’extension à la loi CZM_LIN_REG est faite en suivant exactement le même paradigme.

L’ouverture d’une fissure en mode mixte se caractérise par un critère d’endommagement défini au moyen du saut de déplacement équivalent et de la variable interne \(\alpha\) . Le matériau reste dans le domaine élastique tant que l’on vérifie l’inégalité :

\(f({\mathrm{〚}u\mathrm{〛}}_{\mathit{eq}},\alpha )={\mathrm{〚}u\mathrm{〛}}_{\mathit{eq}}-\alpha \le 0\)

• \({⟦u⟧}_{\mathit{èq}}=\sqrt{{\langle ⟦{u}_{n}⟧\rangle }_{+}^{2}+{⟦{u}_{\tau}⟧}^{2}}\) est le saut de déplacement équivalent,

• \(\mathrm{〚}{u}_{\tau}\mathrm{〛}=\mathrm{〚}u\mathrm{〛}-\mathrm{〚}{u}_{n}\mathrm{〛}n\) est le saut de déplacement tangent,

• :math:`alpha (t)=maxlbrace {alpha}_{0},underset{vin [0,t]}{max}{mathrm{〚}u(v)mathrm{〛}}_{mathit{éq}}rbrace ` est la variable interne de la loi cohésive,

• \({\alpha}_{0}\) est la valeur initiale de \(\alpha\) . Cette valeur est donnée par l’utilisateur par l’intermédiaire du paramètre matériau PENA_ADHERENCE de telle sorte que \({\alpha}_{0}=\frac{{G}_{c}}{{\sigma}_{c}}\text{PENA\_ADHERENCE}\) .

La contrainte cohésive s’écrit alors comme somme d’une contrainte élastique, une contrainte dissipative et une contrainte de pénalisation qui rend compte du contact :

\({t}_{c}=H({〚u〛}_{\text{èq}}-\alpha ){\sigma}_{\mathrm{lin}}+(1-H({〚u〛}_{\text{èq}}-\alpha )){\sigma}_{\mathrm{dis}}+{\sigma}_{\mathrm{pen}}\)

où \(H\) est la fonction indicatrice de \({\mathrm{\Re }}^{+}\)

• \({\sigma}_{\mathit{pen}}=C{\langle ⟦{u}_{n}⟧\rangle }_{-}n\) est la contrainte de pénalisation.

où \(C\) est un coefficient de pénalisation explicité dans [R7.02.11] déterminé à partir du paramètre matériau PENA_CONTACT.

• \({\sigma}_{\mathit{lin}}={\sigma}_{\mathit{dis}}=\frac{{\sigma}_{c}}{\alpha}\exp(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha )⟦u⟧\) est l’expression commune aux contraintes linéaires et dissipative, avec \(\alpha ={\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\) pour \({\sigma}_{\mathit{dis}}\) .

où \({\sigma}_{c}\) est la contrainte critique à la rupture et \({G}_{c}\) est la ténacité énergétique du matériau. Elle correspond en effet à l’énergie nécessaire à l’ouverture complète de l’interface sur une longueur unitaire. Un calcul rapide avec les expressions précédentes permet de confirmer:

\(\underset{-\infty }{\overset{+\infty }{\int}}{t}_{c,\tau }\cdot d⟦{u}_{\tau}⟧+\underset{0}{\overset{+\infty }{\int}}\left({t}_{c}\cdot n\right)d⟦{u}_{n}⟧={G}_{c}\)

On représente sur la figure la contrainte cohésive pour un chargement en mode \(I\) pur en fonction du saut de déplacement normal.

Figure 2.2.1-1 : Evolution de la force de cohésion normale en fonction du saut de déplacement.

Pour \(\mathrm{〚}{u}_{n}\mathrm{〛}<0\) , il est également usuel de définir une force de cohésion équivalente \({t}_{c,\mathrm{eq}}\) grâce à la condition énergétique d’équivalence:

\({t}_{c,\mathit{eq}}{\dot{⟦u⟧}}_{\mathit{eq}}={t}_{c,n}\dot{⟦{u}_{n}⟧}+{t}_{c,\tau }\cdot \dot{⟦{u}_{\tau}⟧}\)

Pour retrouver sa valeur, on dérive \({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\) par rapport au temps.

\(\dot{{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}}=\frac{\dot{\mathrm{〚}{u}_{n}\mathrm{〛}}\mathrm{〚}{u}_{n}\mathrm{〛}+\dot{\mathrm{〚}{u}_{\tau}\mathrm{〛}}\cdot \mathrm{〚}{u}_{\tau}\mathrm{〛}}{{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}}\)

D’où on identifie : \({t}_{c,\mathit{èq}}={t}_{c,n}\frac{{⟦u⟧}_{\mathit{èq}}}{⟦{u}_{n}⟧}=\parallel {t}_{c,\tau }\parallel \frac{{⟦u⟧}_{\mathit{èq}}}{\parallel ⟦{u}_{\tau}⟧\parallel }=\frac{{\sigma}_{c}}{\alpha}\exp\left(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha \right){⟦u⟧}_{\mathit{èq}}\)

La figure représente l’évolution de la force de cohésion équivalente en fonction du saut de déplacement équivalent selon cette loi de comportement.

Figure 2.2.1-2 : Évolution de la force de cohésion équivalente en fonction du saut de déplacement équivalent.

Remarque:

\({t}_{c,\mathit{èq}}=\frac{{\sigma}_{c}}{\alpha}\exp\left(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha \right){⟦u⟧}_{\mathit{èq}}\) On définit parfois un saut de déplacement équivalent \({⟦u⟧}_{\mathit{èq}}=\sqrt{{\langle ⟦{u}_{n}⟧\rangle }_{+}^{2}+{\beta}^{2}{⟦{u}_{\tau}⟧}^{2}}\) où \(\beta\) est un coefficient expérimental qui représente le rapport d’intensité des forces d’ouverture en mode \(I\) et en mode \(\mathit{II}\) . En reprenant alors le raisonnement précédent avec :On déduit les expressions des composantes :

\({t}_{c,n}=\frac{{t}_{c,\mathit{eq}}⟦{u}_{n}⟧}{{⟦u⟧}_{\mathit{èq}}}n=\frac{{\sigma}_{c}}{\alpha}\exp\left(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha \right)⟦{u}_{n}⟧n\) , \({t}_{c,\tau }={\beta}^{2}\frac{{t}_{c,\mathit{eq}}⟦{u}_{\tau}⟧}{{⟦u⟧}_{\mathit{èq}}}={\beta}^{2}\frac{{\sigma}_{c}}{\alpha}\exp\left(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha \right)⟦{u}_{\tau}⟧\)

Lois cohésives mixtes pour éléments quadratiques

Le second type de lois que nous pouvons envisager est, à l’inverse, une loi dans laquelle l’adhérence initiale est parfaite: la pente initiale est infinie. Deux lois cohésives de cette sorte sont disponibles avec XFEM dans Code_Aster, les lois CZM_OUV_MIX et CZM_TAC_MIX dont les caractéristiques sont détaillées dans [R7.02.11]. Nous rappelons ici la loi CZM_OUV_MIX, représentée sur la figure .

En notant \(\delta\) le saut de déplacement (ceci afin de se conformer aux notations de [R7.02.11]), le matériau reste dans le domaine élastique tant que:

\(f(\delta ,\alpha )={\delta}_{n}-\alpha \le 0\)

En revanche, cette fois, la variable interne a une valeur initiale rigoureusement nulle, de sorte que la contrainte n’est plus explicite en fonction du déplacement, comme le montre la figure .

Figure 2.2.2-1: Composante normale du vecteur contrainte en fonction du saut normal pour la loi CZM_OUV_MIX (seuil \(\alpha\) nul à gauche et positif à droite).

Lois cohésives mixtes pour éléments linéaires

La loi a également une pente initiale infinie, mais la décharge est linéaire (voir fig.). Elle est disponible sous le vocabulaire CZM_LIN_MIX. Pour alléger les notations, notons \(w=⟦u⟧\) . Dans le formalisme du lagrangien augmenté (voir documentations [R5.03.52] et [R5.03.54]), une expression générale de la densité d’énergie de surface est:

\(\Pi (w,\lambda )=\varphi (\lambda +rw)-\frac{\lambda \cdot \lambda }{2r}\)

où \(\varphi\) est une fonction dérivable, et \(r\) un paramètre d’augmentation.

La contrainte d’interface est alors classiquement donnée par une dérivée par rapport à la variable primale \({t}_{c}=\frac{\partial \Pi }{\partial w}\) . Il faut une équation additionnelle pour déterminer si on se trouve dans le régime adhérent, et imposer la relation de Dirichlet correspondante si tel est le cas. Dans le formalisme du Lagrangien augmenté, cette loi d’interface additionnelle est déterminée par dérivation par rapport à la variable duale \(\frac{\partial \Pi }{\partial \lambda }=0\) . L’intérêt de la méthode du Lagrangien augmenté est qu’il n’y a pas d’autres inégalités à vérifier au cours de la résolution: le comportement d’interface est entièrement contenu dans ces deux égalités , ce qui permet l’utilisation d’une méthode de Newton pour la résolution.

La force cohésive est alors:

\({t}_{c}\left(\lambda +rw\right)=\frac{\partial \Pi }{\partial w}=r\frac{\partial \varphi }{\partial (\lambda +rw)}\)

La loi d’interface s’écrit alors simplement comme \(\frac{\partial \Pi }{\partial \lambda }=\frac{1}{r}\left[{t}_{c}\left(\lambda +rw\right)-\lambda \right]=0\) , ce qui peut se réécrire plus simplement \(\lambda ={t}_{c}\left(\lambda +rw\right)\) .

Figure 2.2.3-1: Loi cohésive mixte pour des éléments linéaires.

Une force cohésive équivalente est définie par \({\left(\lambda +rw\right)}_{\mathit{eq}}:=\sqrt{{\langle {\lambda}_{n}+r{w}_{n}\rangle }_{+}^{2}+{\left({\lambda}_{\tau}+r{w}_{\tau}\right)}^{2}}\) . Une fonction seuil \(\phi\) est introduite par \(\phi \left({\left(\lambda +rw\right)}_{\mathit{eq}}\right):=\frac{{\left(\lambda +rw\right)}_{\mathit{eq}}-{\sigma}_{c}}{r{w}_{c}-{\sigma}_{c}}\) . Nous pouvons alors définir une variable interne adimensionnelle \(\alpha\) comme maximum de la fonction seuil au cours du temps, projetée sur l’intervalle \([0,1]\) :

\(\stackrel{̃}{\alpha}(t)=\underset{[0,t]}{\max}\phi \left({\left(\lambda +rw\right)}_{\mathit{eq}}\right)\)

\(\alpha ={P}_{[0,1]}(\stackrel{̃}{\alpha})\)

La valeur \(\alpha =0\) (ou \(\stackrel{̃}{\alpha}\le 0\) ) correspond à un matériau sain (zone adhérente), et la valeur \(\alpha =1\) (ou \(\stackrel{̃}{\alpha}\ge 1\) ) à un matériau entièrement fissuré (zone fissurée). Pour des conditions de charge, c’est-à-dire si \(\alpha =\phi \left({\left(\lambda +rw\right)}_{\mathit{eq}}\right)\) , la fonction \(\varphi\) est définie par:

\(\varphi \left(\lambda +rw\right)=2{G}_{c}\left(1-\frac{{\sigma}_{c}}{r{w}_{c}}\right)\alpha \left(1-\frac{\alpha}{2}\right)+\frac{1}{2r}{\langle {\lambda}_{n}+r{w}_{n}\rangle }_{-}^{2}\)

Une expression de la force cohésive \({t}_{c}=\frac{\partial \varphi }{\partial w}\) qui dérive de cette énergie peut alors être obtenue. En définissant une force cohésive équivalente comme \({t}_{c,\mathit{eq}}=\sqrt{{\langle {t}_{c,n}\rangle }_{+}^{2}+{t}_{c,\tau }^{2}}\) , la loi force/ouverture peut être exprimée en terme de quantités équivalentes comme \({t}_{c,\mathit{eq}}=\left(1-{T}_{d}\right){\left(\lambda +rw\right)}_{\mathit{eq}}\) , où \({T}_{d}\) est la fonction d’endommagement qui correspond à un adoucissement linéaire:

\({T}_{d}=\frac{\alpha}{\left(1-\frac{{\sigma}_{c}}{{w}_{c}r}\right)\alpha +\frac{{\sigma}_{c}}{{w}_{c}r}}\)

La relation vectorielle entre force cohésive et saut de déplacement est alors exprimée en termes de composante normale \({t}_{c,n}=\left(1-{T}_{d}\right){\langle {\lambda}_{n}+r{w}_{n}\rangle }_{+}+{\langle {\lambda}_{n}+r{w}_{n}\rangle }_{-}\) (le dernier terme ayant été ajouté pour rendre compte du contact unilatéral) et de composante tangentielle \({t}_{c,\tau }=\left(1-{T}_{d}\right)\left({\lambda}_{\tau}+r{w}_{\tau}\right)\) .

[1]

Force par unité de surface, homogène à une contrainte

Formulations variationnelles

Transformons la forme forte du problème en une formulation faible, mieux adaptée aux éléments finis. Le champ \(u\) doit appartenir à l’ensemble \({V}_{0}\) des champs de déplacements cinématiquement admissibles:

\({V}_{0}=\left\lbrace v\in {H}^{1},v\text{discontinu}à\text{travers}{\Gamma}_{c},v=0\text{sur}{\Gamma}_{u}\right\rbrace\)

Formulation pour la loi cohésive régularisée

Nous notons \(H={H}^{-1/2}(\Gamma )\) pour les lois cohésives. La formulation faible du problème cohésif s’écrit comme suit:

Trouver \(\left(u,{t}_{c}^{+},{t}_{c}^{-}\right)\in {V}_{0}\times H\times H\) tels que:

\({\int}_{\Omega}\sigma (u):\epsilon ({u}^{\text{*}})d\Omega ={\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega +{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}\mathit{d\Gamma }+{\int}_{{\Gamma}^{+}}{t}_{c}^{+}\cdot {u}^{\text{*}+}{\mathit{d\Gamma }}^{+}+{\int}_{{\Gamma}^{-}}{t}_{c}^{-}\cdot {u}^{\text{*}-}{\mathit{d\Gamma }}^{-}\forall {u}^{\text{*}}\in {V}_{0}\)

Puisque \(⟦u⟧\left(x\right)={u}^{+}\left(x\right)-{u}^{-}\left(x\right)\) et \({t}_{c}={t}_{c}^{-}=-{t}_{c}^{+}\) , la formulation faible s’écrit de manière équivalente comme suit:

Trouver \(u\in {V}_{0}\) tel que:

\({\int}_{\Omega}\sigma (u):\epsilon ({u}^{\text{*}})d\Omega ={\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega +{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}d\Gamma -{\int}_{{\Gamma}_{c}}{t}_{c}\cdot ⟦{u}^{\text{*}}⟧d{\Gamma}_{c}\forall {u}^{\text{*}}\in {V}_{0}\)

L’implémentation actuelle comporte en fait une partie post-traitement, afin de donner une valeur aux multiplicateurs de Lagrange inusités:

Trouver \(\left(u,{\lambda}_{n},{\lambda}_{\tau}\right)\in {V}_{0}\times H\times H\) tels que:

\(\forall \left({u}^{\text{*}},{\lambda}_{n}^{\text{*}},{\lambda}_{\tau}^{\text{*}}\right)\in {V}_{0}\times H\times H\)

Équation d’équilibre	\(\begin{array}{c}{\int}_{\Omega}\sigma (u):\epsilon ({u}^{\text{*}})d\Omega -{\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega -{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}\mathit{d\Gamma }\\ +{\int}_{{\Gamma}_{c}}{t}_{c,n}\cdot {⟦{u}^{\text{*}}⟧}_{n}{\mathit{d\Gamma }}_{c}+{\int}_{{\Gamma}_{c}}{t}_{c,\tau }\cdot {⟦{u}^{\text{*}}⟧}_{\tau}{\mathit{d\Gamma }}_{c}=0\end{array}\)
Post-traitement partie normale	\({\int}_{\mathit{\Gamma c}}\left({\lambda}_{n}-{t}_{c,n}\cdot n\right){\lambda}_{n}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\)
Post-traitement partie tangentielle	\({\int}_{\mathit{\Gamma c}}\left({\lambda}_{\tau}-{t}_{c,\tau }\right){\lambda}_{\tau}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\)

Ainsi, les multiplicateurs \({\lambda}_{n}\) et \({\lambda}_{\tau}\) n’interviennent pas dans la résolution. Ils servent uniquement à stocker les contraintes cohésives de manière explicite.

Inconvénients d’une loi cohésive régularisée

Afin d’évaluer la capacité d’une telle loi cohésive à décrire convenablement la zone adhérente, un test d’ouverture d’une inclusion circulaire est mené. Géométrie et chargement sont donnés en figure (dimensions en millimètres). Il s’agit d’une plaque sous tension en déformations planes, avec \(E=36.56\mathit{GPa}\) et \(\nu =0.2\) . Une inclusion circulaire est insérée dans la plaque, qui est supposée pouvoir s’ouvrir selon une loi cohésive de contrainte critique \({\sigma}_{c}=2.7\mathit{MPa}\) , et ténacité énergétique \({G}_{c}=0.095{\mathit{N.mm}}^{-1}\) .

Figure 3.2-1: Test d’ouverture d’une inclusion circulaire.

Figure 3.2-2: Force cohésive normale le long de l’inclusion.

Figure 3.2-3: Saut de déplacement normal le long de l’inclusion.

Le saut de déplacement normal et la force qui en est déduite sont tracés sur les figures et , pour un chargement \({u}_{0}=0.04\mathit{mm}\) , pour des maillages conformes ou non à la surface fissurée, et pour des paramètres de pénalisation \({\alpha}_{0}={10}^{-2}\) et \({\alpha}_{0}={10}^{-6}\) . Comme nous pouvons le voir sur les figures, les lois cohésives régularisées présentent trois difficultés:

1. si la pénalisation est trop faible, le saut de déplacement le long de l’inclusion est faux: on observe une ouverture qui n’est pas physique (marque 1 sur la fig.);

2. si, au contraire, la pénalisation est trop élevé, des oscillations parasites de la contrainte cohésive appraissent dans la zone adhérente (marque 2 sur la fig.). Ce problème a été abondamment signalé dans la littérature (voir par exemple [bib3]). Originellement remarqué dans le cadre du contact pour Code_Aster, et désigné sous le terme de “condition LBB”, il est traité dans la documentation [R5.03.54], §6. En un mot, il vient du fait que l’espace de discrétisation des forces cohésives devient trop riche en comparaison de celui du déplacement, lorsque la raideur de la loi devient telle que l’on s’approche d’une condition de Dirichlet;

3. une mauvaise évaluation du régime de fonctionnement (adhérent ou dissipatif) peut arriver en conséquence de ces oscillations (marque 3 sur la fig).

Pour conclure, les lois pénalisées échouent à donner une représentation fidèle des zones adhérentes, puisque un paramètre de pénalisation \({\alpha}_{0}^{-1}\) trop petit implique une ouverture non-physique de la zone adhérente, ce qui conduit à une solution fausse, tandis qu’un paramètre trop grand \({\alpha}_{0}^{-1}\) génère des oscillations parasites du fait d’un problème de stabilité, dès lors que le paramètre de pénalisation commence à devenir suffisement important pour décrire convenablement l’adhérence.

En considérant le test d’arrachement de la figure (ce test est numérique, il ne correspond à aucune réalité physique), il est possible de mettre en évidence de façon encore plus évidente le phénomène d’oscillations.

Figure 3.2-4: Test d’arrachement: géométrie et chargement.

Les paramètres pour ce test sont les suivants: \(E=30\mathit{GPa}\) , \(\nu =0.2\) , \({\sigma}_{c}=0.2\mathit{MPa}\) , \({G}_{c}=2\mathit{N.}{\mathit{mm}}^{-1}\) et \({\alpha}_{0}={10}^{-6}\) . Nous avons repésenté en figure la contrainte d’interface pour un maillage conforme, qui sert de référence, et la contrainte pour un maillage non conforme. Le phénomène d’oscillations parasites d’origine numérique y est particulièrement visible sur toute la zone adhérente.

Figure 3.2-5: Test d’arrachement: contrainte cohésive le long de l’interface.

Espace réduit pour la discrétisation de la contrainte d’interface

Pour une description détaillée de la discrétisation des inconnues de contact, le lecteur peut se référer à la documentation [R5.03.54], §5. En bref, les composantes initiales du multiplicateur sont définies sur les noeuds sommets \(K\) des éléments parents intersectés (voir fig.). L’implémentation de tels éléments de contact est détaillée en [R5.03.54], §4. On impose ensuite des rekations d’égalité entre certaines de ces composantes initiales afin d’aboutir à un nombre plus faible \({N}_{\lambda}\) de degrés de liberté effectivement indépendants. Ces relations sont portées par certaines arêtes intersectées \(V\) , dites arêtes vitales: un degré de liberté \(I\) réellement indépendant est partagé par un groupe de noeuds de \(K\) (voir fig.), ce qui produit une fonction de forme de contact étendue \({\psi}_{I}:=\sum_{i\in I}{N}_{i}\) (voir fig.). L’algorithme de sélection de telles arêtes vitales, et donc de contruction de l’espace réduit, est détaillé dans le documentation [R5.03.54], §6. Le champ de multiplicateurs est ensuite obtenu par interpolation sur les éléments parents et le multiplicateur discret est la trace de ce champ sur l’interface:

\({M}_{h}:=\left\lbrace \sum_{I}{\mu}_{I}{\psi}_{I}{\text{|}}_{\Gamma},{\mu}_{I}\in {ℝ}^{d}\right\rbrace\)

Figure 3.3-1: Maillage non conforme à la fissure et espace réduit de multiplicateurs.

Formulation pour une loi cohésive mixte pour éléments quadratiques

Par oppositionà la formulation précédente, le traitement d’une telle loi va nécessiter une vraie formulation à plusieurs champs, dans le sens où unchamp dual vectoriel \(\lambda\) va effectivemententrerdans la formulation, au lieu d’être un artifice de post-traitement comme en § 3.1 . Cette formulation suit un raisonnement énergétique, expliqué en détail dans la documentation [R3.06.13]. Résumons-enen les principaux points:

On écrit que l’ouverture de la fissure coûte une énergie proportionnelle à la surface à ouvrir, soit:

\({E}_{\text{fr}}(\delta )=\underset{\Gamma}{\int}\Pi (\delta )\text{dS}\)

où \(\Pi (\delta )\) est la densité d’énergie cohésive. Pour la loi CZM_OUV_MIX, nous avons par exemple \(\Pi (\delta )={\int}_{0}^{{\delta}_{n}}{t}_{c,n}(\delta ')d\delta '\) .

Le champ de discontinuité \(\delta\) apparaissant dans les expressions précédentes est alors défini comme un champ à part entière, intégré dans la formulation en tant que nouvelle inconnue. L’énergie totale s’écrit alors:

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

La solution du problème consiste alors en la minimisation de cette énergie totale sous la contrainte que \(\delta\) corresponde au saut de déplacement. On cherche:

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

Afin de résoudre celle-ci, nous introduisons le lagrangien associé au problème, auquel nous ajoutons un terme d’augmentation dont l’utilité apparaîtra par la suite:

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

Nous pouvons alors écrire la première condition d’optimalité de ce lagrangien:

\(\forall {\delta}^{\text{*}}\underset{\Gamma}{\int}\left[{t}_{c}-\lambda +r\left(\delta -⟦u⟧\right)\right]\cdot {\delta}^{\text{*}}=0\text{avec}{t}_{c}\in \partial \Pi \left(\delta \right)\)

Cette équationfait intervenir la contrainte cohésive \({t}_{c}\) . Or nous ne disposons comme expression de \({t}_{c}\) que d’une loi de comportement locale. Cettepremière équation doit donc, pour avoir du sens,être discrétisée d’une façon qui permette de se ramener à une expression locale. Ceci est possible si \(\delta\) estdiscrétisé par collocation aux points de Gauss de l’interface, de coordonnées \({x}_{g}\) .

En effet, avec une telle discrétisation, la résolution de la première condition d’optimalitérevient à satisfaire la loi cohésive en chaque point de collocation:

\({t}_{c}\left({\delta}_{g},{\alpha}_{g}\right)={\lambda}_{g}+r\left(⟦{u}_{g}⟧-{\delta}_{g}\right)\)

où nous avons noté \({\lambda}_{g}=\lambda ({x}_{g})\) , par exemple, les valeurs d’un champ aux points de Gauss, et où \({t}_{c}({\delta}_{g},{\alpha}_{g})\) suit la loi . La traduction graphique de cette loi de comportement est la suivante: la solution correspond à l’intersection de la fonction linéaire \(\delta \to {\lambda}_{g}+r⟦{u}_{g}⟧-r\delta\) (avec une pente négative donnée par le coefficient de pénalité \(r\) ) avec le graphe \({t}_{c}(\delta ,\alpha )\) . Nous voyons alors que pour \(r\) assez grand, la solution est unique, d’où l’intérêt d’avoir augmenté le Lagrangien.

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

Le champ \(\delta\) s’écrit donc localement comme une fonction de \(\lambda +r⟦u⟧\) , que nous appellerons multiplicateur augmenté et noterons \(p\) . Par conséquent, il disparaît des champs inconnus du problème. La formulation est alors donnée par les deux conditions d’optimalité du lagrangien restantes:

Trouver \((u,\lambda )\in {V}_{0}\times H\) tels que:

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}})\in {V}_{0}\times H\)

Équation d’équilibre	\(\begin{array}{c}{\int}_{\Omega}\sigma (u):\epsilon ({u}^{\text{*}})d\Omega -{\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega -{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}d\Gamma \\ +{\int}_{{\Gamma}_{c}}\left[\lambda +r\left(⟦u⟧-\delta \left(p\right)\right)\right]\cdot ⟦{u}^{\text{*}}⟧d\Gamma =0\text{avec}p=\lambda +r⟦u⟧\end{array}\)
Loi d’interface	\({\int}_{{\Gamma}_{c}}\left[⟦u⟧-\delta (p)\right]\cdot {\lambda}^{\text{*}}d\Gamma =0\)

En ce qui concerne la discrétisation de cesdeux champs d’inconnues, une discrétisation simple respectant la condition de stabilitéinf-sup, et consistante avec la discrétisation de \(\delta\) par collocation, consiste à discrétiser ledéplacement avec des éléments P2et lemultiplicateur d’une façonP1 sur un espace réduit adapté à X-FEM,comme détaillé en § 3.3 ou en [R5.03.54], §6.

Formulation pour une loi cohésive mixte pour éléments linéaires

Comme démontré auparavant, les quantités à l’interface doivent être définies sur un espace réduit \({M}_{h}\) afin de prévenir des oscillations parasites d’origine numérique dans les phases adhérentes. Ceci s’applique aux contraintes cohésives \(\lambda\) et \({t}_{c}\) . De fait, \({t}_{c}\) dépend de \(\lambda +rw\) où \(w\) est le saut de déplacement. En conséquence de cette augmentation, notamment présente lors des phases adhérentes, \(w\) doit lui aussi être écrit sur l’espace consistant \({M}_{h}\) .

Remarque:

Une autre façon de s’en convaincre consiste à considérer un problème d’adhérence pure, l’adhérence étant décrite par une loi régularisée. Ce problème décrit en effet le comportement du terme d’augmentation lors des phases adhérentes. Pour éviter des oscillations, nous avons vu qu’une loi pénalisée doit être convenablement écrite sur l’espace réduit, de sorte que le système à résoudre est: \(\left[\begin{array}{ccc}{K}_{u}^{1}& 0& 0\\ 0& {K}_{u}^{2}& {A}^{T}\\ 0& A& \frac{1}{r}C\end{array}\right]\left\lbrace \begin{array}{c}\lbrace a\rbrace \\ \lbrace b\rbrace \\ \lbrace \lambda \rbrace \end{array}\right\rbrace =\left\lbrace \begin{array}{c}\lbrace {L}_{\text{méca}}^{1}\rbrace \\ \lbrace {L}_{\text{méca}}^{2}\rbrace \\ 0\end{array}\right\rbrace\)

où:

• nous regroupons dans \(\lbrace {L}_{\text{méca}}\rbrace\) et \([{K}_{u}]\) la discrétisation de tous les termes qui ne sont pas des termes d’interface;

• \({C}_{IK}={\int}_{\Gamma}{\psi}_{I}{\psi}_{K}d\Gamma\) ;

• \({A}_{Ij}={\int}_{\Gamma}{\psi}_{I}{\varphi}_{j}d\Gamma\) ;

• \(\lbrace a\rbrace\) est le vecteur regroupant les DDL de déplacement classiques;

• \(\lbrace b\rbrace\) celui des DDL enrichis;

• \(\lbrace \lambda \rbrace\) regroupe des DDL du multiplicateur à l’interface, discrétisation qui s’effectue sur l’espace réduit \({M}_{h}\) ;

• \(r\) est le paramétre de pénalisation.

Dans une formulation pénalisée consistante, nous pouvons interpréter la contrainte d’interface comme proportionnelle à un saut consistant \(\lbrace w\rbrace\) : \(\lbrace \lambda \rbrace =r\lbrace w\rbrace\) *. Avec l’expression matricielle précédente, nous pouvons alors identifier* \(\lbrace w\rbrace ={[C]}^{-1}[A]\lbrace b\rbrace\) .

Nous pourrions alors condenser le système matriciel précédent pour éliminer \(\lbrace \lambda \rbrace\) :

\([{K}_{u}^{2}]\lbrace b\rbrace +{[A]}^{T}\lbrace \lambda \rbrace =[{K}_{u}^{2}]\lbrace b\rbrace +{[A]}^{T}\left(r{[C]}^{-1}[A]\lbrace b\rbrace \right)\)

Nous voyons alors que la résolution nécéssiterait l’assemblage et le calcul d’une matrice \({[A]}^{T}{[C]}^{-1}[A]\) .

Un tel produit matriciel n’est pas réalisable dans un calcul élémentaire, à plus forte raison compte-tenu du caractère non-local de \({M}_{h}\) . Une telle opération globale va donc à l’encontre de l’architecture informatique de Code_Aster, et doit être proscrite dans la mesure du possible.

Pour parer au problème, \(w\) est introduit comme une nouvelle inconnue du problème, qui n’est pas discrétisé comme \(⟦u⟧\) mais en est une projection sur un espace réduit \({M}_{h}\) . Décomposé à l’aide de ce champ désormais disponible dans la formulation, le terme de stabilisation précédent pourra être déterminé par assemblage d’un calcul élémentaire. L’énergie totale du problème s’écrit:

\(E(u,\lambda ,w)=\frac{1}{2}{\int}_{\Omega}ϵ(u):C:ϵ(u)d\Omega -{\int}_{{\Gamma}_{g}}g\cdot ud{\Gamma}_{g}+{\int}_{\Gamma}\Pi \left(w,\lambda \right)d\Gamma\)

La solution du problème continu consiste en une minimisation sous contraintes d’égalité \(\left(u,w,\lambda \right)=\underset{{w}^{\text{*}}=⟦{u}^{\text{*}}⟧}{\mathit{argmin}}E\left({u}^{\text{*}},{\lambda}^{\text{*}},{w}^{\text{*}}\right)\) . Nous pouvons écrire le Lagrangien associé comme:

\(L(u,w,\lambda ,\mu )=\frac{1}{2}{\int}_{\Omega}ϵ(u):C:ϵ(u)d\Omega -{\int}_{{\Gamma}_{g}}g\cdot ud{\Gamma}_{g}+{\int}_{\Gamma}\Pi (w,\lambda )d\Gamma +{\int}_{\Gamma}\mu \cdot \left(⟦u⟧-w\right)d\Gamma\)

L’écriture des conditions d’optimalité de ce Lagrangien conduit à la formulation variationnelle suivante:

Équation d’équilibre	\(\forall {u}^{\text{*}}\in {V}_{h},{\int}_{\Omega}\sigma (u):ϵ({u}^{\text{*}})d\Omega -{\int}_{{\Gamma}_{g}}g\cdot {u}^{\text{*}}d{\Gamma}_{g}+{\int}_{\Gamma}\mu \cdot ⟦{u}^{\text{*}}⟧d\Gamma =0\)
Projection du saut de déplacement	\(\forall {\mu}^{\text{*}}\in {M}_{h},{\int}_{\Gamma}\left(⟦u⟧-w\right)\cdot {\mu}^{\text{*}}d\Gamma =0\)
Expression de la force cohésive	\(\forall {w}^{\text{*}}\in {M}_{h},-{\int}_{\Gamma}\left[\mu -{t}_{c}(\lambda +rw)\right]\cdot {w}^{\text{*}}d\Gamma =0\)
Loi d’interface	\(\forall {\lambda}^{\text{*}}\in {M}_{h},-{\int}_{\Gamma}\frac{\left[\lambda -{t}_{c}(\lambda +rw)\right]}{r}\cdot {\lambda}^{\text{*}}d\Gamma =0\)

Stratégie de résolution

Méthode de Newton-Raphson

La stratégie de résolution n’est autre qu’une simple méthode de Newton. Contrairement au cas du contact-frottement [R5.03.52] ou [R5.03.54], il n’y a ni boucles de point fixe, ni champs de signes. La seule opération à réaliser en plus des itérations de Newton est donc l’actualisation de la variable interne.

• Pour un pas de temps:

  • Itérations de Newton Calcul de la matrice tangente et du second membre

  • Fin des itérations de Newton

Actualisation de la variable interne \(\alpha\)

On pourrait légitimement se demander pourquoi la variable interne n’est pas actualisée au cours des itérations de Newton. En fait, comme il s’agit d’un paramètre mesurant l’irréversibilité, et déterminé par un maximum au cours du temps, il ne doit être actualisé qu’à chaque pas de temps convergé. En effet, dans le cas contraire, si ce paramètre excède sa valeur d’équilibre lors d’une itération de Newton, l’algorithme de Newton sera alors incapable de le diminuer pour trouver la valeur d’équilibre.

Dans le cas unidimensionnel, la méthode de Newton est un processus itératif permettant d’approcher les zéros d’une fonction continue et dérivable. On se ramène à la résolution de \(F(x)=0\) . On construit une suite de points \({x}^{k}\) en faisant un développent de Taylor de \(F\) au voisinage de \({x}^{k}\) , ce qui donne au premier ordre:

\(F({x}^{k+1})\approx F({x}^{k})+{F}^{'}({x}^{k})({x}^{k+1}-{x}^{k})\)

En notant \({\mathrm{\delta x}}^{k}={x}^{k+1}-{x}^{k}\) l’incrément entre deux itérations successives, l’équation linéarisée à l’itération \(k+1\) est alors la suivante:

\({F}^{'}({x}^{k}){\mathrm{\delta x}}^{k}=-F({x}^{k})\)

Dans le cas de la méthode des éléments finis, \({F}^{'}({x}^{k})\) s’apparente à la matrice tangente, qui peut être calculée à chaque itération si nécessaire, \({\mathrm{\delta x}}^{k}\) est le vecteur des incréments des inconnues nodales, et \(F({x}^{k})\) est le second membre. On note que \({F}^{'}({x}^{k})\) et \(F({x}^{k})\) ne font intervenir que des quantités de l’itération \(k\) , qui sont donc des quantités connues.

Différenciation de la loi cohésive

Pour donner un exemple de différenciation de la loi cohésive couplant les modes, considérons la loi cohésive régularisée CZM_EXP_REG. La différentielle de \({t}_{c}(\mathrm{〚}u\mathrm{〛})\) sera \(\frac{\partial {t}_{c}}{\partial \mathrm{〚}u\mathrm{〛}}\cdot \mathrm{〚}\delta u\mathrm{〛}\) avec:

Or \(\frac{\partial {t}_{c}}{\partial \mathrm{〚}u\mathrm{〛}}=H({\mathrm{〚}u\mathrm{〛}}_{\text{èq}}-\alpha )\frac{\partial {\sigma}_{\mathit{lin}}}{\partial \mathrm{〚}u\mathrm{〛}}+(1-H({\mathrm{〚}u\mathrm{〛}}_{\text{èq}}-\alpha ))\frac{\partial {\sigma}_{\mathit{dis}}}{\partial \mathrm{〚}u\mathrm{〛}}+\frac{\partial {\sigma}_{\mathit{pen}}}{\partial \mathrm{〚}u\mathrm{〛}}\)

Nous réutilisons les expressions de ces trois dérivées partielles qui sont données dans [R7.02.11], avec \(\delta =⟦u⟧\) . Dans l’expression de \({\sigma}_{\mathit{dis}}\) , il faut écrire \(\alpha ={\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\) , qui devient ainsi une variable à prendre en compte dans la dérivation. En pratique, dans le code, on distingue quatre cas pour la clarté de lecture:

• • \({⟦u⟧}_{\mathit{èq}}\ge \alpha \mathrm{et}⟦{u}_{n}⟧\ge 0\) (dissipatif non-contactant). Nous avons alors : \(\frac{\partial {t}_{c}}{\partial ⟦u⟧}={\sigma}_{c}\exp\left(-\frac{{\sigma}_{c}}{{G}_{c}}{⟦u⟧}_{\mathit{èq}}\right)\left(\frac{\text{Id}}{{⟦u⟧}_{\mathit{èq}}}-\frac{⟦u⟧}{{⟦u⟧}_{\mathit{èq}}}\otimes \frac{⟦u⟧}{{⟦u⟧}_{\mathit{èq}}}\left(\frac{{\sigma}_{c}}{{G}_{c}}+\frac{1}{{⟦u⟧}_{\mathit{èq}}}\right)\right)\)

• \({⟦u⟧}_{\mathit{èq}}<\alpha \mathrm{et}⟦{u}_{n}⟧<0\) (élastique contactant). Avec \(({\tau}_{1},{\tau}_{2})\) une base du plan tangent, nous avons: \(\frac{\partial {t}_{c}}{\partial ⟦u⟧}=\frac{\partial {\sigma}_{\mathit{lin}}(⟦{u}_{\tau}⟧)}{\partial ⟦u⟧}+\frac{\partial {\sigma}_{\mathit{pen}}(⟦{u}_{n}⟧)}{\partial ⟦u⟧}=Cn\otimes n+\frac{{\sigma}_{c}}{\alpha}\exp(\frac{-{\sigma}_{c}}{{G}_{c}}\alpha )\left({\tau}_{1}\otimes {\tau}_{1}+{\tau}_{2}\otimes {\tau}_{2}\right)\)

• \({⟦u⟧}_{\mathit{èq}}\ge \alpha \mathrm{et}⟦{u}_{n}⟧<0\) (dissipatif contactant). Par un raisonnement similaire en remplaçant \({\sigma}_{\mathit{lin}}\) par \({\sigma}_{\mathit{dis}}\) , nous obtenons : \(\frac{\partial {t}_{c}}{\partial ⟦u⟧}=Cn\otimes n+\exp\left(\frac{-{\sigma}_{c}}{{G}_{c}}{⟦u⟧}_{\mathit{èq}}\right)\left[\frac{{\sigma}_{c}}{{⟦u⟧}_{\mathit{èq}}}\left({\tau}_{1}\otimes {\tau}_{1}+{\tau}_{2}\otimes {\tau}_{2}\right)-\left(\frac{{\sigma}_{c}^{2}}{{G}_{c}}+\frac{{\sigma}_{c}}{{⟦u⟧}_{\mathit{èq}}}\right)\frac{{⟦u⟧}_{\tau}\otimes {⟦u⟧}_{\tau}}{{{⟦u⟧}_{\mathit{èq}}}^{2}}\right]\)

• \({⟦u⟧}_{\mathit{èq}}<\alpha \mathrm{et}⟦{u}_{n}⟧\ge 0\) (élastique non-contactant). \(\frac{\partial {t}_{c}}{\partial ⟦u⟧}=\frac{{\sigma}_{c}}{\alpha}\exp(\frac{-{\sigma}_{c}}{{G}_{c}}\alpha )\text{Id}\)

Linéarisation du problème

Écriture intégrale pour une formulation avec loi cohésive régularisée

Le système linéaire des trois équations à l’itération de Newton \(k+1\) s’écrit de la manière suivante (pour ne par alourdir l’écriture, les références à l’itération de Newton sont omises, car de manière évidente, les inconnues sont notées avec un \(\delta\) devant, et les champs tests sont désormais notés avec une étoile):

Trouver \(\left(\mathit{\delta u},{\mathit{\delta \lambda }}_{n},{\mathit{\delta \lambda }}_{\tau}\right)\in {V}_{0}\times H\times H\) tel que:

\(\forall \left({u}^{\text{*}},{\lambda}_{n}^{\text{*}},{\lambda}_{\tau}^{\text{*}}\right)\in {V}_{0}\times H\times H\)

Équation d’équilibre	\(\begin{array}{c}{\int}_{\Omega}\sigma (\mathit{\delta u}):\epsilon ({u}^{\text{*}})d\Omega \\ +{\int}_{{\Gamma}_{c}}\left[\frac{\partial {t}_{c,n}}{\partial {⟦u⟧}_{n}}{⟦\delta u⟧}_{n}+\frac{\partial {t}_{c,n}}{\partial {⟦u⟧}_{\tau}}{⟦\delta u⟧}_{\tau}\right]{⟦{u}^{\text{*}}⟧}_{n}{\mathit{d\Gamma }}_{c}\\ +{\int}_{{\Gamma}_{c}}\left[\frac{\partial {t}_{c,\tau }}{\partial {⟦u⟧}_{n}}{⟦\delta u⟧}_{n}+\frac{\partial {t}_{c,\tau }}{\partial {⟦u⟧}_{\tau}}{⟦\delta u⟧}_{\tau}\right]{⟦{u}^{\text{*}}⟧}_{\tau}{\mathit{d\Gamma }}_{c}\\ =-{\int}_{\Omega}\sigma (u):\epsilon ({u}^{\text{*}})d\Omega +{\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega +{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}{\mathit{d\Gamma }}_{t}\\ -{\int}_{{\Gamma}_{c}}\left({t}_{c,n}{⟦{u}^{\text{*}}⟧}_{n}+{t}_{c,\tau }{⟦{u}^{\text{*}}⟧}_{\tau}\right){\mathit{d\Gamma }}_{c}\end{array}\)
Interface: partie normale	\({\int}_{\mathit{\Gamma c}}{\lambda}_{n}^{\text{*}}\left({\lambda}_{n}+{\text{δλ}}_{n}-{t}_{c}\cdot n\right){\mathit{d\Gamma }}_{c}=0\)
Interface: partie tangentielle	\({\int}_{\mathit{\Gamma c}}{\lambda}_{\tau}^{\text{*}}\left({\lambda}_{\tau}+\delta {\lambda}_{\tau}-{t}_{c,\tau }\right){\mathit{d\Gamma }}_{c}=0\)

Écriture intégrale pour une formulation avec loi cohésive mixte pour éléments quadratiques

Le système linéaire des trois équations à l’itération de Newton \(k+1\) s’écrit de la manière suivante (pour ne par alourdir l’écriture, les références à l’itération de Newton sont omises, car de manière évidente, les inconnues sont notées avec un \(\delta\) devant, et les champs tests sont désormais notés avec une étoile):

Trouver \((\mathit{\delta u},\mathit{\delta \lambda })\in {V}_{0}\times H\) tel que:

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}})\in {V}_{0}\times H\)

Équation d’équilibre

\(\begin{array}{c}{\int}_{\Omega}\sigma \left(\delta u\right):ϵ\left({u}^{\text{*}}\right)d\Omega +{\int}_{\Gamma}\left(\mathit{Id}-r\frac{\partial \delta }{\partial p}\right)\cdot \delta \lambda \cdot ⟦{u}^{\text{*}}⟧d\Gamma \\ +{\int}_{\Gamma}r\left(\mathit{Id}-r\frac{\partial \delta }{\partial p}\right)\cdot ⟦\delta u⟧\cdot ⟦{u}^{\text{*}}⟧d\Gamma \\ =-{\int}_{\Omega}\sigma \left(u\right):ϵ\left({u}^{\text{*}}\right)d\Omega +{\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega +{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}d\Gamma \\ -{\int}_{{\Gamma}_{c}}\left[\lambda +r\left(⟦u⟧-\delta (p)\right)\right]\cdot ⟦{u}^{\text{*}}⟧d\Gamma \end{array}\)

Loi d’interface

\(\begin{array}{c}{\int}_{{\Gamma}_{c}}\left(1-r\frac{\partial \delta }{\partial p}\right)\cdot ⟦u⟧\cdot {\lambda}^{\text{*}}d\Gamma -{\int}_{\Gamma}\frac{\partial \delta }{\partial p}\cdot \delta \lambda \cdot \lambda \text{*}d\Gamma \\ =-{\int}_{\Gamma}\left(⟦u⟧-\delta \right)\cdot {\lambda}^{\text{*}}d\Gamma \end{array}\)

Écriture intégrale pour une formulation avec loi cohésive mixte pour éléments linéaires

Le système linéaire à résoudre pour une itération de Newton s’écrit:

Trouver \(\left(\mathit{\delta u},\delta \mu ,\mathit{\delta w},\mathit{\delta \lambda }\right)\in {V}_{0}\times H\times H\times H\) tel que:

\(\forall \left({u}^{\text{*}},{\mu}^{\text{*}},{w}^{\text{*}},{\lambda}^{\text{*}}\right)\in {V}_{0}\times H\times H\times H\)

Équation d’équilibre	\(\begin{array}{c}{\int}_{\Omega}\sigma (\delta u):ϵ({u}^{\text{*}})d\Omega +{\int}_{\Gamma}\delta \mu \cdot ⟦{u}^{\text{*}}⟧d\Gamma \\ =-{\int}_{\Omega}\sigma (u):ϵ({u}^{\text{*}})d\Omega +{\int}_{{\Gamma}_{g}}g\cdot {u}^{\text{*}}d{\Gamma}_{g}-{\int}_{\Gamma}\mu \cdot ⟦{u}^{\text{*}}⟧d\Gamma \end{array}\)
Projection du saut de déplacement	\({\int}_{\Gamma}\left(⟦\delta u⟧-\delta w\right)\cdot {\mu}^{\text{*}}d\Gamma =-{\int}_{\Gamma}\left(⟦u⟧-w\right)\cdot {\mu}^{\text{*}}d\Gamma\)
Contrainte cohésive	\(\begin{array}{c}-{\int}_{\Gamma}\left[\delta \mu -\frac{\partial {t}_{c}}{\partial (\lambda +rw)}\cdot \left(\delta \lambda +r\delta w\right)\right]\cdot {w}^{\text{*}}d\Gamma \\ ={\int}_{\Gamma}\left[\mu -{t}_{c}(\lambda +rw)\right]\cdot {w}^{\text{*}}d\Gamma \end{array}\)
Loi d’interface	\(\begin{array}{c}-{\int}_{\Gamma}\left[\frac{\delta \lambda }{r}-\frac{\partial {t}_{c}}{\partial (\lambda +rw)}\cdot \left(\frac{\delta \lambda }{r}+\delta w\right)\right]\cdot {\lambda}^{\text{*}}d\Gamma \\ ={\int}_{\Gamma}\frac{\left[\lambda -{t}_{c}(\lambda +rw)\right]}{r}\cdot {\lambda}^{\text{*}}d\Gamma \end{array}\)

Expression matricielle du problème

Loi cohésive régularisée

Écriture matricielle du problème linéarisé

En reprenant les notations de [R5.03.54], et en considérant l’écriture unifiée adoptée pour les lois de contact frottement et cohésives régularisées, le système linéarisé tel qu’il est résolu à l’itération \(k+1\) de Newton peut se mettre sous forme matricielle:

Équation d’équilibre	\(\begin{array}{c}\left\lbrace {u}^{\text{*}}\right\rbrace \left[{K}_{\text{méca}}\right]\left(\mathit{\delta u}\right)+\left\lbrace {u}^{\text{*}}\right\rbrace \left[{K}_{\mathit{coh}}^{u}\right]\left(\mathit{\delta u}\right)\\ =\left\lbrace {u}^{\text{*}}\right\rbrace \left({L}_{\text{méca}}^{1}\right)+\left\lbrace {u}^{\text{*}}\right\rbrace \left({L}_{\text{coh}}^{1}\right)\\ \end{array}\)
Post-traitement contrainte cohésive normale	\(\left\lbrace {\lambda}_{n}^{\text{*}}\right\rbrace \left[C\right]\left({\text{δλ}}_{n}\right)=\left\lbrace {\lambda}_{n}^{\text{*}}\right\rbrace \left({L}_{\text{post}}^{1}\right)+\left\lbrace {\lambda}_{n}^{\text{*}}\right\rbrace \left({L}_{\text{coh}}^{2}\right)\)
Post-traitement contrainte cohésive tangentielle	\(\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace \left[{F}_{r}\right]\left({\mathit{\delta \lambda }}_{\tau}\right)=\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace \left({L}_{\text{post}}^{2}\right)+\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace \left({L}_{\text{coh}}^{3}\right)\)

où les vecteurs colonne sont notés \((x)\) et les vecteurs ligne \(\left\lbrace x\right\rbrace ={(x)}^{T}\) .Ce système peut être mis sous la forme matricielle suivante:

\(\left[\begin{array}{c}{{\rm K}}_{\text{méca}}+{K}_{\text{coh}}^{u}\\ 0\\ 0\end{array}\begin{array}{c}0\\ C\\ 0\end{array}\begin{array}{c}0\\ 0\\ F\end{array}\right]\left(\begin{array}{c}\mathit{\delta u}\\ {\text{δλ}}_{n}\\ {\mathit{\delta \lambda }}_{\tau}\end{array}\right)=\left(\begin{array}{c}{L}_{\text{méca}}^{1}+{L}_{\text{coh}}^{1}\\ {L}_{\text{post}}^{1}+{L}_{\text{coh}}^{2}\\ {L}_{\text{post}}^{2}+{L}_{\text{coh}}^{3}\end{array}\right)\)

L’inconnue est l’incrément par rapport à l’itération de Newton précédente. On a volontairement omis la référence au numéro de l’itération de Newton.

\({K}_{\text{méca}}\) est la matrice de rigidité mécanique définie au paragraphe [§3.2] de [R7.02.12].

\({K}_{\mathit{coh}}^{u}\) est la matrice de rigidité due aux forces de cohésion.

\(C\) est la matrice de post-traitement pour la direction normale.

\(F\) est la matrice de post-traitement pour la(les) direction(s) tangentielle(s).

\({L}_{\text{méca}}^{1}\) est le second membre représentant les forces internes et les incréments de chargements.

\({L}_{\text{post}}^{1}\) et \({L}_{\text{post}}^{2}\) sont les seconds membres de post-traitement.

\({L}_{\text{coh}}^{1}\) , \({L}_{\text{coh}}^{2}\) et \({L}_{\text{coh}}^{3}\) sont les seconds membres dus aux forces de cohésion.

Remarque:

On rappelle que le système résolu par Code_Aster n’est pas du type \([K][U]=[F]\) mais du type \([K][U]+[F]=0\) . Il existe donc un signe moins entre les seconds membres donnés dans ce document et ceux codés dans les fichiers fortran.

Expression des matrices et vecteurs élémentaires

La matrice \(C\) nécessaire au post-traitement a pour expression (les indices \(n\) sur \(\lambda\) , pour indiquer la composante normale de la contrainte cohésive, sont omis):

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left[C\right]}_{ij}{(\text{δλ})}_{j}={\int}_{\Gamma}{\psi}_{i}{\lambda}_{i}^{\text{*}}{\psi}_{j}{\text{δλ}}_{j}\mathit{d\Gamma }\)

Le vecteur \({L}_{\text{post}}^{1}\) nécessaire au post-traitement a pour expression (indices \(n\) omis sur sur \(\lambda\) ):

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left({L}_{\text{post}}^{1}\right)}_{i}={\int}_{\Gamma}{\psi}_{i}{\lambda}_{i}^{\text{*}}{\lambda}^{k-1}\mathit{d\Gamma }\)

La matrice \(F\) nécessaire au post-traitement a pour expression:

\({\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace }_{i}{\left[F\right]}_{ij}{\left(\delta {\lambda}_{\tau}\right)}_{j}={\int}_{\Gamma}{\psi}_{i}\left\lbrace \begin{array}{cc}{\lambda}_{\tau ,i}^{1\text{*}}& {\lambda}_{\tau ,i}^{2\text{*}}\end{array}\right\rbrace \cdot \left[\begin{array}{cc}{\tau}_{i}^{1}{\tau}_{j}^{1}& {\tau}_{i}^{1}{\tau}_{j}^{2}\\ {\tau}_{i}^{2}{\tau}_{j}^{1}& {\tau}_{i}^{2}{\tau}_{j}^{2}\end{array}\right]\cdot {\psi}_{j}\left(\begin{array}{c}{\lambda}_{\tau ,j}^{1}\\ {\lambda}_{\tau ,j}^{2}\end{array}\right)\mathit{d\Gamma }\)

Le vecteur \({L}_{\text{post}}^{2}\) nécessaire au post-traitement a pour expression:

\(\begin{array}{c}{\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace }_{i}{\left({L}_{\text{post}}^{2}\right)}_{i}={\int}_{\Gamma}{\psi}_{i}\left\lbrace {\lambda}_{\tau ,i}^{1\text{*}}{\tau}_{i}^{1}+{\lambda}_{\tau ,i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace {\lambda}_{\tau}^{k-1}\mathit{d\Gamma }\end{array}\)

Cherchons à présent à écrire \({K}_{\mathit{coh}}^{u}\) . Soient deux directions de la base fixe \(X\) et \(Y\) , de vecteurs unitaires \({e}_{X}\) et \({e}_{Y}\) . Introduisons la matrice tangente de la loi cohésive dans la base fixe \({K}^{\mathit{gl}}\) de coefficients \({K}_{XY}^{\mathit{gl}}={e}_{X}\cdot \frac{\partial {t}_{c}}{\partial ⟦u⟧}\cdot {e}_{Y}\) . Avec l’expression de \(\frac{\partial {t}_{c}}{\partial \mathrm{〚}u\mathrm{〛}}\) donnée au [§ 4.2 ], nous disposons de la matrice tangente de la loi cohésive \({K}^{\mathit{loc}}\) dans la base locale (voir doc [R7.02.11]). Nous obtenons alors \({K}^{\mathit{gl}}\) par \({K}^{\mathit{gl}}={Q}^{T}\cdot {K}_{\mathit{loc}}\cdot Q\) , où \(Q\) est une matrice de passage orthonormale définie par:

\(Q=\left[\begin{array}{ccc}{n}_{X}& {n}_{Y}& {n}_{Z}\\ {\tau}_{X}^{1}& {\tau}_{Y}^{1}& {\tau}_{Z}^{1}\\ {\tau}_{X}^{2}& {\tau}_{Y}^{2}& {\tau}_{Z}^{2}\end{array}\right]\)

Soient \(i\) et \(j\) deux nœuds enrichis. Soit \(\Gamma\) l’intersection des supports de \(i\) et \(j\) . La matrice \([{K}_{\mathit{coh}}^{u}]\) est alors donnée par:

\({\lbrace {u}^{\text{*}}\rbrace }_{i}{[{K}_{\mathit{coh}}^{u}]}_{ij}{\lbrace \delta u\rbrace }_{j}={\int}_{\Gamma}2{\varphi}_{i}{b}_{i}^{\text{*}}\cdot {K}^{\mathit{gl}}\cdot {b}_{j}2{\varphi}_{j}d\Gamma\)

Les seconds membres de cohésion ont les expressions suivantes:

\({\lbrace {u}_{i}\rbrace }^{\text{*}}{\left({L}_{\mathit{coh}}^{1}\right)}_{i}=-{\int}_{\Gamma}2{\varphi}_{i}{b}_{i}^{\text{*}}\cdot {t}_{c}^{k-1}d\Gamma\)

où l’on peut exprimer \({t}_{c}\) dans la base globale par \({t}_{c}^{\mathit{glo}}={Q}^{T}\cdot {t}_{c}^{\mathit{loc}}\) , et:

\({\left\lbrace {\lambda}_{n}^{\text{*}}\right\rbrace }_{i}{\left({L}_{\text{coh}}^{2}\right)}_{i}=-{\int}_{\Gamma}{\psi}_{i}{\lambda}_{n,i}^{\text{*}}\left({t}_{c,n}^{k-1}\right)\mathit{d\Gamma }\)

\(\begin{array}{c}{\left\lbrace {\lambda}_{\tau}^{\text{*}}\right\rbrace }_{i}{\left({L}_{\text{coh}}^{3}\right)}_{i}=-{\int}_{\Gamma}{\psi}_{i}\left\lbrace {\lambda}_{\tau ,i}^{1\text{*}}{\tau}_{i}^{1}+{\lambda}_{\tau ,i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace \cdot {t}_{c,\tau }^{k-1}\mathit{d\Gamma }\end{array}\)

où \(k-1\) représente l’indice de l’itération de Newton précédente.

Loi cohésive mixte pour éléments quadratiques

Ecriture matricielle du problème avec loi cohésive mixte

Le système matriciel tel qu’il est résolu à l’itération \(k+1\) de Newton peut se mettre sous forme matricielle suivante:

\(\left[\begin{array}{c}{{\rm K}}_{\text{méca}}+{{\rm A}}_{u}\\ {\rm A}\end{array}\begin{array}{c}{{\rm A}}^{T}\\ C\end{array}\right](\begin{array}{c}\mathit{\delta u}\\ \text{δλ}\end{array})=(\begin{array}{c}{L}_{\text{méca}}^{1}+{L}_{\text{coh}}^{1}\\ {L}_{\text{coh}}^{2}\end{array})\)

Expression des matrices élémentaires de cohésion :

Soient deux directions de la base fixe \(X\) et \(Y\) , de vecteurs unitaires \({e}_{X}\) et \({e}_{Y}\) . Nous introduisons comme précédemment la matrice tangente de la loi cohésive dans la base fixe \({K}^{\mathit{gl}}\) de coefficients \({K}_{\mathit{XY}}^{\mathit{gl}}={e}_{X}\cdot \frac{\partial \delta }{\partial p}\cdot {e}_{Y}\) . Nous disposons, par la loi de comportement cohésive, de la matrice tangente \({K}_{\mathit{loc}}\) dans la base locale (voir doc [R7.02.11]). Nous obtenons alors \({K}_{\mathit{gl}}\) par \({K}^{\mathit{gl}}={Q}^{T}\cdot {K}^{\mathit{loc}}\cdot Q\) , où \(Q\) est une matrice de passage orthonormale définie par:

\(Q=\left[\begin{array}{ccc}{n}_{X}& {n}_{Y}& {n}_{Z}\\ {\tau}_{X}^{1}& {\tau}_{Y}^{1}& {\tau}_{Z}^{1}\\ {\tau}_{X}^{2}& {\tau}_{Y}^{2}& {\tau}_{Z}^{2}\end{array}\right]\)

Ayant introduit ces notations, nous avons:

\({\lbrace {u}^{\text{*}}\rbrace }_{i}{[{A}_{u}]}_{ij}{\lbrace \delta u\rbrace }_{j}={\int}_{\Gamma}2{\varphi}_{i}{b}_{i}^{\text{*}}r\left(1-r{K}^{\mathit{gl}}\right)2{\varphi}_{j}{b}_{j}d\Gamma\)

Par ailleurs, nous choisissons de discrétiser \(\lambda\) en base locale. Le coefficient \(\left(1,1\right)\) de la matrice \(\lbrace {\lambda}_{n}^{\text{*}}\rbrace {[{A}_{u}]}_{ij}{\lbrace \delta {u}_{X}\rbrace }_{j}\) est alors \({\int}_{\Gamma}{\varphi}_{i}{\lambda}_{i}^{\text{*}}n\cdot \left(1-r\frac{\partial \delta }{\partial p}\right)\cdot {e}_{X}2{\varphi}_{j}{b}_{j}d\Gamma\) . Pour les coefficients \(\left(2,1\right)\) et \(\left(3,1\right)\) , la formule est la même en remplaçant \(n\) par \({\tau}_{1}\) et \({\tau}_{2}\) , respectivement. En exploitant les notations introduites dans le paragraphe précédent, nous en déduisons:

\({\lbrace {\lambda}^{\text{*}}\rbrace }_{i}{[A]}_{ij}{\lbrace \delta u\rbrace }_{j}={\int}_{\Gamma}{\varphi}_{i}{\lambda}_{i}^{\text{*}}\left(1-r{K}^{\mathit{loc}}\right)\cdot Q2{\varphi}_{j}{b}_{j}d\Gamma\)

Quant à la matrice \([C]\) , elle s’écrit simplement:

\({\lbrace {\lambda}^{\text{*}}\rbrace }_{i}{[C]}_{ij}{\lbrace \delta \lambda \rbrace }_{j}=-{\int}_{\Gamma}{\varphi}_{i}{\lambda}_{i}^{\text{*}}{K}^{\mathit{loc}}{\varphi}_{j}{\lambda}_{j}d\Gamma\)

Expression des vecteurs élémentaires de cohésion :

Le coefficient selon \(X\) de \({\left({L}_{\mathit{coh}}^{1}\right)}_{i}\) a pour expression \({\int}_{\Gamma}{\varphi}_{i}\left(-\lambda \cdot {e}_{X}-r\left(⟦u⟧-\delta \right)\cdot {e}_{X}\right)d\Gamma\) . Avec les notations indroduites dans la partie précédentes, nous déduisons:

\({\lbrace {u}_{i}\rbrace }^{\text{*}}{\left({L}_{\mathit{coh}}^{1}\right)}_{i}={\int}_{\Gamma}2{\varphi}_{i}{b}_{i}^{\text{*}}\cdot \left(-{Q}^{T}\cdot \lambda -r\left(⟦u⟧-{Q}^{T}\cdot \delta \right)\right)d\Gamma\)

où \(\delta\) et \(\lambda\) sont donnés en base locale, et \(⟦u⟧\) en base fixe.

Le coefficient \(1\) de \({\left({L}_{\mathit{coh}}^{2}\right)}_{i}\) s’écrit \({\int}_{\Gamma}\left(-⟦u⟧\cdot n+\delta \cdot n\right)d\Gamma\) . Pour les coefficients \(2\) et \(3\) , la formule est la même en remplaçant \(n\) par \({\tau}_{1}\) et \({\tau}_{2}\) , respectivement. Nous en déduisons:

\({\lbrace \lambda \rbrace }_{i}^{\text{*}}{\left({L}_{\mathit{coh}}^{2}\right)}_{i}={\int}_{\Gamma}{\lambda}_{i}^{\text{*}}\cdot \left(-Q⟦u⟧+\delta \right)d\Gamma\)

où \(\delta\) est donné en base locale, et \(⟦u⟧\) en base fixe.

Loi cohésive mixte pour éléments linéaires

Les composantes des inconnues \(u\) et \(\mu\) sont définies dans une base fixe \(\left({e}_{X},{e}_{Y},{e}_{Z}\right)\) , tandis que les composantes de \(w\) et \(\lambda\) sont définies dans la base locale \(\left(n,{\tau}_{1},{\tau}_{2}\right)\) à la surface fissurée \(\Gamma\) en chaque point \(x\in \Gamma\) , de sorte que :

\(w(x)=\sum_{i=1}^{{N}_{\lambda}}{\psi}_{I}(x)\left({w}_{I,n}n(x)+{w}_{I,\tau 1}{\tau}_{1}(x)+{w}_{I,\tau 2}{\tau}_{2}(x)\right)\)

Une définition analogue vaut pour \(\lambda\) . Pour un degré de liberté \(I\) de l’espace réduit (voir § 3.3 ), il est possible de déterminer les composantes \({t}_{c,n}^{I},{t}_{t,\tau 1}^{I},{t}_{c,\tau 2}^{I}\) de la force cohésive à partir de \(\left({w}_{I,n},{w}_{I,\tau 1},{w}_{I,\tau 2}\right)\) , \(\left({\lambda}_{I,n},{\lambda}_{I,\tau 1},{\lambda}_{I,\tau 2}\right)\) et de la loi cohésive du § 2.2.3 . Ces composantes ne sont pas destinées à être associées à une direction particulière \(I\) du degré de liberté, mais destinées à être reliées au sens faible à la contrainte globale \(\mu\) écrite en base fixe.

Le système matriciel peut se mettre sous la forme suivante:

\(\left[\begin{array}{cccc}{K}^{\mathit{uu}}& {\left({K}^{\mu u}\right)}^{T}& 0& 0\\ {K}^{\mu u}& 0& {\left(-{K}^{w\mu }\right)}^{T}& 0\\ 0& {-K}^{w\mu }& {D}^{\mathit{ww}}& {\left({D}^{\lambda w}\right)}^{T}\\ 0& 0& {D}^{\lambda w}& {D}^{\lambda \lambda }\end{array}\right]\left(\begin{array}{c}\delta u\\ \delta \mu \\ \delta w\\ \delta \lambda \end{array}\right)=\left(\begin{array}{c}{L}_{\text{méca}}+{L}_{\mu}^{1}\\ {L}_{u}-{L}_{w}\\ -{L}_{\mu}^{2}+{L}_{\mathit{coh}}^{1}\\ -{L}_{\lambda}+{L}_{\mathit{coh}}^{2}\end{array}\right)\)

where:

• \({K}^{uu}\) est la matrice de rigidité volumique;

• \({K}^{\mu u}\) et \({K}^{w\mu }\) sont des matrices discrétisant les opérateurs «mortier», la dernière gérant également le changement de base;

• les matrices \(D\) sont toutes diagonales par blocs: pour \(I\) et \(J\) deux DDL de Lagrange distincts, elles vérifient \({D}_{\mathit{IJ}}=0\) .

Donnons les expressions des matrices et seconds membres qui ne dépendent pas de la loi cohésive:

\({\lbrace u\rbrace }_{i}^{\text{*}}{\left[{K}^{\mu u}\right]}_{iJ}{(\delta \mu )}_{J}={b}_{i}^{\text{*}}\cdot (\delta {\mu}_{J}){\int}_{\Gamma}2{\psi}_{J}{\varphi}_{i}d\Gamma\)

En introduisant la matrice de changement de base orthonormale \(Q\) définie comme précédemment, nous avons:

\({\lbrace w\rbrace }_{I}^{\text{*}}{\left[{K}^{w\mu }\right]}_{IJ}{(\delta \mu )}_{J}={w}_{I}^{\text{*}}\cdot \left({\int}_{\Gamma}{\psi}_{I}{\psi}_{J}Qd\Gamma \right)\cdot (\delta {\mu}_{J})\)

\({\lbrace u\rbrace }_{i}^{\text{*}}{\left({L}_{\mu}^{1}\right)}_{i}=-{b}_{i}^{\text{*}}\cdot {\int}_{\Gamma}2{\varphi}_{i}\mu d\Gamma\)

\({\lbrace \mu \rbrace }_{I}^{\text{*}}{\left({L}_{u}\right)}_{I}=-{\mu}_{I}^{\text{*}}\cdot {\int}_{\Gamma}{\psi}_{I}⟦u⟧d\Gamma\)

\({\lbrace \mu \rbrace }_{I}^{\text{*}}{\left({L}_{w}\right)}_{I}=-{\mu}_{I}^{\text{*}}\cdot {\int}_{\Gamma}{\psi}_{I}{Q}^{T}\cdot wd\Gamma\)

\({\lbrace w\rbrace }_{I}^{\text{*}}{\left({L}_{\mu}^{2}\right)}_{I}={w}_{I}^{\text{*}}\cdot {\int}_{\Gamma}{\psi}_{I}Q\cdot \mu d\Gamma\)

Détaillons à présent la discrétisation des quantités d’interface. Les matrices sont diagonales par blocs, et ont pour expressions:

\({\lbrace w\rbrace }_{I}^{\text{*}}{\left[{D}^{ww}\right]}_{II}{(\delta w)}_{I}={w}_{I}^{\text{*}}(\delta {w}_{I})r\frac{\partial {t}_{c}}{\partial (\lambda +rw)}({\lambda}_{I}+r{w}_{I}){\int}_{\Gamma}{\psi}_{I}d\Gamma\) \({\lbrace \lambda \rbrace }_{I}^{\text{*}}{\left[{D}^{\lambda w}\right]}_{II}{(\delta w)}_{I}={\lambda}_{I}^{\text{*}}(\delta {w}_{I})\frac{\partial {t}_{c}}{\partial (\lambda +rw)}({\lambda}_{I}+r{w}_{I}){\int}_{\Gamma}{\psi}_{I}d\Gamma\)

\({\lbrace \lambda \rbrace }_{I}^{\text{*}}{\left[{D}^{\lambda \lambda }\right]}_{II}{(\delta \lambda )}_{I}={\lambda}_{I}^{\text{*}}(\delta {\lambda}_{I})\frac{1}{r}\left(\frac{\partial {t}_{c}}{\partial (\lambda +rw)}({\lambda}_{I}+r{w}_{I})-1\right){\int}_{\Gamma}{\psi}_{I}d\Gamma\)

Les vecteurs élémentaires ont pour expressions:

\({\lbrace w\rbrace }_{I}^{\text{*}}{\left({L}_{\mathit{coh}}^{1}\right)}_{I}={w}_{I}^{\text{*}}{t}_{c}({\lambda}_{I}+r{w}_{I}){\int}_{\Gamma}{\psi}_{I}d\Gamma\)

\({\lbrace w\rbrace }_{I}^{\text{*}}{\left({L}_{\mathit{coh}}^{2}\right)}_{I}=\frac{{w}_{I}^{\text{*}}}{r}{t}_{c}({\lambda}_{I}+r{w}_{I}){\int}_{\Gamma}{\psi}_{I}d\Gamma\)

Propagation d’une fissure cohésive

Dans cette méthode, des modèles de zones cohésives sont introduits sur des surfaces potentielles de fissuration étendues. Ainsi, la loi cohésive séparera naturellement les domaines adhérents et ouverts, ce qui demande une formulation pour la loi cohésive qui permette une représentation fidèle de zones adhérentes de grande taille. Le problème a été discuté au § 3 . Une actualisation implicite du front de propagation peut alors être menée, ce qui fait l’originalité de la méthode proposée ici. Pour les algorithmes de propagation, nous utilisons une loi cohésive mixte avec éléments linéaires CZM_LIN_MIX, qui est introduite dans le modèle par la commande DEFI_CONTACT, avec les mot-clés ZONE/RELATION.

Supposons une structure présentant une surface fissurée, décrite par une loi cohésive, pour laquelle on connaisse le front cohésif séparant la zone cohésive effective de la surface potentielle de fissuration, c’est-à-dire la surface au sein de laquelle la fissure pourra se propager au prochain pas de chargement (voir fig.). La réalisation d’une itération du calcul de propagation comporte alors, dans l’ordre, les étapes suivantes (voir fig.):

1. la réalisation du calcul mécanique au pas de temps suivant (commande STAT_NON_LINE),

2. la détection du nouveau front cohésif (commande PROPA_FISS, opération DETECT_COHESIF),

3. la détermination de la direction de propagation le long de ce nouveau front (commande CALC_G),

4. la construction d’une nouvelle zone de fissuration potentielle au delà de ce front, à partir des directions déterminées en c), selon une extension de longueur fixe (commande PROPA_FISS, opération PROPA_COHESIF),

5. l’actualisation du modèle suivant cette nouvelle zone (commande MODI_MODELE_XFEM), et si nécessaire, la projection de l’état cohésif courant sur ce nouveau modèle (mot-clé ETAT_INIT/COHE dans STAT_NON_LINE). Cet état servira d’état initial pour le calcul mécanique de l’étape a) de l’itération suivante.

Figure 6-a: Propagation d’une fissure sur trajet inconnu: étapes pour la réalisation d’un pas de propagation.

Description d’une fissure cohésive

Afin de pouvoir localiser le front cohésif, on définit une discontinuité de type “COHESIF” dans DEFI_FISS_XFEM. Par rapport à une discontinuité de type “INTERFACE”, elle doit contenir des informations supplémentaires pour localiser le front cohésif: les structures de données relatives à la description d’un front (coordonnées des points, base covariante…) et une level-set supplémentaire. En effet, sachant que la position de la surface de fissuration potentielle est décrite avec une level-set normale par \([{\varphi}_{n}=0]\) , le front cohésif est classiquement localisé par une level-set «tangente» supplémentaire, comme \([{\varphi}_{n}=0]\cap [:ref:`{\varphi}_{t}=0 <{\varphi}_{t}=0>\)]` (voir fig.). Pour une description détaillée de l’utilisation de level-sets avec XFEM, voir la documentation [R7.02.12], §2.

Figure 6.1-a: Description du front cohésif par l’intermédiaire de level-sets.

Contrairement au type “FISSURE”, le type “COHESIF” n’a pas d’enrichissement singulier en point, du fait de la régularité des champs cohésifs. La zone de fissuration potentielle s’arrête donc sur les faces (en 3D) ou les arêtes (en 2D, voir fig.) des éléments.

Figure 6.1-b: Caractéristiques d’une fissure de type COHESIF, en sortie d’une opération DETECT_COHESIF de PROPA_FISS

Remarque:

En fait, pour le type COHESIF , le front de propagation ne coïncide pas nécessairement toujours avec l’iso-zéro de level-set tangente. En effet, lors des phases de calcul mécanique (étape a de la fig.), l’iso-zéro de level-set tangente est utilisé pour localiser la fin de la zone potentielle de fissuration, et donc de l’enrichissement (voir fig.), tandis que le front de propagation localise toujours le front cohésif (à partir duquel l’interface est effectivement ouverte, voir fig.).

Figure 6.1-c: Caractéristiques d’une fissure de type COHESIF, en sortie d’une opération PROPA_COHESIF de PROPA_FISS

Remarque – initialisation de la procédure:

L’initialisation d’une étude de propagation (par l’étape a de la fig.) se fait en donnant une zone potentielle de fissuration et un front de propagation initiaux (voir fig.). Ce dernier peut correspondre à un coin de la structure, une simple ligne sur un bord libre, ou un fond de fissure préexistante (dans ce dernier cas, celle-ci doit être maillée). En pratique, on définit une discontinuité initiale de type COHESIFpar DEFI_FISS_XFEM (voir fig.). Le front initial est renseigné par le mot-clé GROUP_MA_BORD.

Figure 6.1-d: Caractéristiques d’une fissure de type COHESIF, en sortie d’une commande DEFI_FISS_XFEM

Réalisation du calcul mécanique

La loi CZM_LIN_MIX et la formulation utilisées ont été détaillées dans les cinq premières parties de cette documentation. La longueur de la zone cohésive est estimée par \({l}_{c}=\frac{E{G}_{c}}{(1-{\nu}^{2}){\sigma}_{c}^{2}}\) . Il est conseillé de choisir une taille de maille \(h\) telle que \({l}_{c}\sim \text{quelques}h\) pour que la zone cohésive soit bien résolue.

Il est conseillé d’utiliser le pilotage du chargement, ce qui permet de contrôler l’avancée de fissure à chaque pas de propagation. En effet, de façon générale, lorsque l’on utilise le pilotage du chargement, l’incrément de déplacement \(\Delta u\) est lié au pas de temps \(\Delta t\) par la relation (voir documentation [R5.03.80]):

\(f(\Delta u)=\frac{\Delta t}{\text{COEF\_MULT}}\)

où:

• \(f\) est la fonction de pilotage, qui dépend de la quantité de l’on cherche à contrôler,

• COEF_MULT est la valeur renseignée sous le mot-clé du même nom, sous le mot-clé facteur PILOTAGE de la commande STAT_NON_LINE.

Plus particulièrement pour la loi cohésive CZM_LIN_MIX, la fonction de pilotage est:

\(f(\Delta u)=\frac{⟦\Delta u⟧}{{w}_{c}}\) où \({w}_{c}=\frac{2{G}_{c}}{{\sigma}_{c}}\) est le saut de déplacement critique.

La longueur de fissure créée sur un pas de propagation est donc approximativement \(\Delta l=\frac{{l}_{c}}{2}\frac{\Delta t}{\text{COEF\_MULT}}\) . Une fois de plus, il est conseillé de choisir \(\Delta l\sim \text{quelques}h\) . Si possible, \(\Delta l\) doit être très inférieur au rayon de courbure attendu pour la surface fissurée (afin de capter précisément la géométrie de la fissure). A défaut, ce critère doit être vérifié à postériori.

Détection du front

Afin de déterminer le front de fissure à partir du résultat de calcul cohésif, nous allons attribuer aux noeuds du maillage un statut pour savoir s’ils sont sains ou déjà endommagés. Ce statut doit être un scalaire et ne pas être modifié dans des situations de décharge. Par conséquent, il semble naturel d’utiliser la variable interne \(\alpha\) de la loi cohésive. Cependant, comme toute variable interne mesurant un processus irréversible, \(\alpha\) n’évolue pas durant les phases adhérentes (\(\alpha =0\) , voir § 2.2.3 ). En l’état, elle ne permet donc pas de déterminer un champ scalaire dont l’iso-zéro définirait le front de fissure. En revanche, on peut utiliser à cet effet la quantité \(\stackrel{̃}{\alpha}\text{}\) , telle que \(\alpha ={P}_{[0,1]}(\stackrel{̃}{\alpha})\) . Pour rappel (cf § 2.2.3 ), on a \(\stackrel{̃}{\alpha}(t)=\underset{[0,t]}{\max}\phi \left({\left(\lambda +rw\right)}_{\mathit{eq}}\right)\) avec \(\phi \left({\left(\lambda +rw\right)}_{\mathit{eq}}\right):=\frac{{\left(\lambda +rw\right)}_{\mathit{eq}}-{\sigma}_{c}}{r{w}_{c}-{\sigma}_{c}}\) . Dés lors, \(\stackrel{̃}{\alpha}\le 0\) correspond à un matériau sain, et \(\stackrel{̃}{\alpha}\ge 0\) à un matériau fissuré.

La construction d’un front de fissure à partir du champ pertinent précédent se fait en assimilant \(\stackrel{̃}{\alpha}\text{}\) à une level-set tangente \({\varphi}_{t}\) . La procédure consiste alors à:

1. récupérer le champ \(\stackrel{̃}{\alpha}\text{}\) , déjà défini comme un champ aux nœuds pour la formulation mixte avec éléments linéaires;

2. repérer les éléments intersectés (voir fig., étapes a et b) par l’iso-zéro de \(\stackrel{̃}{\alpha}\text{}\) . Pour le détail, voir la documentation [R7.02.12], §3.2.5.1 et §3.2.5.2;

3. déterminer les points d’intersection de l’iso-zéro du champ avec les faces (voir fig.c). Pour un détail de la procédure, voir la documentation Code_Aster [R7.02.12], §2.4.1.

Comme mis en évidence avec les références à la documentation Code_Aster, toutes les routines réalisant ces opérations existent déjà, vu que le champ agit comme une level-set tangente, ce qui rend l’implémentation particulièrement rapide.

Figure 6.3-a : Calcul des levels sets aux noeuds projetés sur le fond de fissure.

Un pré-requis pour la phase de détection est la production d’une level-set tangente \({\varphi}_{t}\) régulière: elle doit pouvoir l’être suffisamment pour être utilisée par les algorithmes de propagation de fissure aux pas de propagation suivants. Le champ \({\varphi}_{t}=\stackrel{̃}{\alpha}\text{}\) et le front qui en découle sont encore trop irréguliers pour cet usage.

Afin de gommer ces irrégularités, il convient d’appliquer au front une procédure de lissage.

Pour ce faire, nous utilisons la connaissance de l’ancien front, supposé plus régulier. L’avancée de fissure qui produit le front nouvellement détecté à partir de l’ancien est déterminée à postériori (fig.a-b). La direction étant connue à partir de la propagation précédente, on applique un algorithme d’actualisation de level-sets pour produire une level-set tangente plus régulière (fig.c-d).

En détail, les étapes de cette procédure sont donnés ci-dessous.

1. Pour chaque point \(P\) de l” ancien front, on détermine une distance de propagation en calculant la distance séparant \(P\) du point d’intersection entre le plan \(\left(P,{n}_{P},{t}_{P}\right)\) et le front «brut» nouvellement détecté (voir fig., étape a). Etant donnée une longueur d’influence \(d\) , cette distance est calculée de la façon suivante:

\({A}_{P}=\frac{\sum_{i}\omega ({d}_{i}){A}_{i}}{\sum_{i}\omega ({d}_{i})}\) avec \(\omega ({d}_{i})=\left\lbrace \begin{array}{c}\exp\left(\frac{-{d}_{i}^{2}}{2{d}^{2}}\right)\text{si}{d}_{i}\le \mathrm{2d}\\ 0\text{sinon}\end{array}\right\rbrace\)

Avertissement. Attention, ceci est à distinguer d’une projection 3D de \(P\) sur le nouveau front, qui ne fournirait pas le même résultat. Dans notre cas, on cherche bien l’avancée dans le plan normal à l’ancien front passant par \(P\) .

2. On lisse le champ des normes de vitesse ainsi obtenu (voir fig., étape b).

Remarques. Si le champ est suffisamment régulier, cette étape pourra être omise. La direction de propagation ayant été déterminée au pas de temps précédent (voir partie 6.4 ), on connaît complètement la vitesse de propagation sur l’ancien front (voir fig., étape c), ce qui permet:

3. La mise en oeuvre d’un algorithme d’actualisation des level-sets (voir partie 6.5 ), pour lequel on effectue uniquement les opérations concernant la level-set tangente (voir fig., étape d).

Conclusion. La level-set tangente actualisée fournit la position d’un front régulier et la base covariante associée.

Figure 6.3-b : Lissage du front de fissure.

Remarque – à propos de l’avancée de fissure aux points extrémités:

Figure 6.3-c : Avancée de fissure aux points extrémités. Lors de la détermination de l’avancée de fissure à postériori, on ne tient pas compte des valeurs pour les points extrémités de l’ancien front. En effet, la base covariante \(\left({n}_{P},{t}_{P}\right)\) y est corrigée pour être tangente à la surface libre (voir § 6.5 ). Suite à cette opération, la valeur de l’avancée peut y être sensiblement différente des valeurs aux points voisins (voir fig.), ce qui a pour effet de fausser l’algorithme de propagation qui suit (fig. , étape d). Au lieu de faire le calcul d’avancée pour ces points, on leur attribue la valeur du point dont le propagé reste à l’intérieur le plus proche.

Remarque – choix de la longueur d’influence \(d\) :

Afin de limiter le nombre deparamètres devant être réglés par l’utilisateur, \(d\) est automatiquement déduit du paramètre NB_POINT_FOND donné en argument à la procédure PROPA_FISS et à la commande CALC_G , par \(d=\frac{\mid T\mid }{\text{NB\_POINT\_FOND}}\) où \(\mid T\mid\) est la longueur du front de fissure.

On peut utiliser un paramètre NB_POINT_FOND différent d’un pas de propagation à un autre (par exemple si le front s’agrandit), mais ce réglage reste à la discrétion de l’utilisateur, il n’est pas automatique. Une fois de plus, il est conseillé de choisir \(d\sim \text{quelques}h\) . Si possible, \(d\) doit être très inférieur aux rayons de courbures attendus à la fois pour la surface fissurée et pour le front (afin de capter précisément leur géométries respectives). A défaut, ce critère doit être vérifié à postériori.

Détermination de la direction de propagation

Dans la littérature, la direction de propagation peut être déduite de deux types d’informations: soit les champs de contraintes au voisinage de la pointe, soit les facteurs d’intensité des contraintes équivalents.

La réutilisation de critéres directionnels basés sur les facteurs d’intensité des contraintes même pour des calculs cohésifs a été suggérée par [bib4, bib5]. Les auteurs y expliquent qu’il s’agit d’une hypothèse raisonnable pour le béton, puisque les expériences montrent que le trajet de fissuration, contrairement à la courbe de réponse force/déplacement, n’est presque pas sensible à la taille de la zone cohésive. On pourrait calculer ces facteurs d’intensité des contraintes à partir d’un calcul élastique de fissure libre supplémentaire – on enlève les forces cohésives du modèle – comme en [bib6, bib7].

On peut cependant se passer de cette étape supplémentaire: par exemple, Moës et Belytschko [bib8] les calculent directement à partir du résultat avec forces cohésives. L’argument justifiant cette simplification est le suivant (Planas et Elices [bib9]): loin de la zone cohésive, les champs mécaniques sont proches de ceux ceux d’une fissure libre «équivalente». Des facteurs d’intensité des contraintes é quivalents peuvent alors être définis par des intégrales de contour qui entourent la totalité de la zone cohésive.

Nous définissons un champ virtuel d’extension de la fissure \(\theta\) , tel que \(\parallel \theta \parallel =1\) , tangent à la surface fissurée en chaque point de celle-ci et dirigé selon \(t\) le long du front (voir figures et ). Soit \(C\) un contour qui entoure la totalité de la zone cohésive, et \({\Gamma}_{C}\) la partie de \(\Gamma\) qui se situe entre le front de fissure et les extrémités de \(C\) (voir fig.), de sorte que \(C\cup {\Gamma}_{C}^{+}\cup {\Gamma}_{C}^{-}\) est un contour fermé (fig.). Comme détaillé en [bib8], une intégrale J peut être définie sur ce contour fermé comme \(J={J}_{\text{coh}}+{J}_{\text{ext}}\) avec \({J}_{\text{ext}}:=-{\int}_{C}\theta \cdot E\cdot nd\Gamma\) et \({J}_{\text{coh}}={\int}_{{\Gamma}_{C}}\theta \cdot ⟦E⟧\cdot nd\Gamma\) , le tenseur d’Eshelby étant noté \(E={\nabla u}^{T}\cdot \sigma -\frac{1}{2}\left(\sigma :ϵ\right)1\) et \(⟦E⟧={E}^{+}-{E}^{-}\) . Du fait que \({\sigma}^{+}\cdot n={\sigma}^{-}\cdot n={t}_{c}\) et \(\theta \cdot n=0\) (pour rappel, \(n\) est orienté de \(-\) vers \(+\) ), l’expression de \({J}_{\text{coh}}\) se réduit à:

\({J}_{\text{coh}}={\int}_{{\Gamma}_{C}}{t}_{c}\cdot \nabla ⟦u⟧\cdot \theta d\Gamma\)

L’intégrande de \({J}_{\text{coh}}\) s’annule en dehors de la zone cohésive, puisque \({t}_{c}=0\) sur la zone ouverte et \(⟦u⟧=0\) sur la zone adhérente . Ainsi, \({J}_{\text{coh}}\) est indépendant du contour \(C\) pourvu que ce dernier entoure la zone cohésive, et s’écrit comme:

\({J}_{\text{coh}}={\int}_{\Gamma}{t}_{c}\cdot \nabla ⟦u⟧\cdot \theta d\Gamma\)

De plus, puisqu’il n’y a pas de singularité au voisinage du front en présence de forces cohésives, nous avons \(J=0\) , ce qui confirme que \({J}_{\text{ext}}\) ne dépend pas de \(C\) tant que ce dernier englobe la zone cohésive.

Par la formule d’Irwin, nous avons par ailleurs \({J}_{\text{ext}}=-{J}_{\text{coh}}=\frac{1-{\nu}^{2}}{E}\left({K}_{I,\mathit{eq}}^{2}+{K}_{\mathit{II},\mathit{eq}}^{2}\right)+\frac{1}{2\mu }{K}_{\mathit{III},\mathit{eq}}^{2}\) . En décomposant la force cohésive par \({t}_{c}={t}_{c,n}n+{t}_{c,t}t+{t}_{c,b}b\) et en adoptant la notation \(⟦\nabla u⟧\cdot \theta =\frac{\partial {⟦u⟧}_{n}}{\partial \theta }n+\frac{\partial {⟦u⟧}_{t}}{\partial \theta }t+\frac{\partial {⟦u⟧}_{b}}{\partial \theta }b\) , les facteurs d’intensité des contraintes équivalents peuvent s’écrire:

\({K}_{I,\mathit{eq}}^{2}=-\frac{E}{1-{\nu}^{2}}{\int}_{\Gamma}\frac{\partial {⟦u⟧}_{n}}{\partial \theta }{t}_{c,n}d\Gamma\)

\({K}_{\mathit{II},\mathit{eq}}^{2}=-\frac{E}{1-{\nu}^{2}}{\int}_{\Gamma}\frac{\partial ⟦{u}_{t}⟧}{\partial \theta }{t}_{c,t}d\Gamma\)

\({K}_{\mathit{III},\mathit{eq}}^{2}=-2\mu {\int}_{\Gamma}\frac{\partial ⟦{u}_{b}⟧}{\partial \theta }{t}_{c,b}d\Gamma\)

Physiquement, ces expressions quantifient l’expression de l’énergie dissipée suivant chacun des trois modes de rupture, si l’on fait l’hypothèse d’une propagation homothétique, dans la direction \(\theta\) , et que la taille de la zone cohésive reste petite par rapport à celle de l’échantillon .

Figure 6.4-a : Définitions d’intégrales invariantes en présence d’une zone cohésive.

Le critère alors adopté est le critère de contrainte circonférentielle maximale de Erdogan et Sih [bib12] (voir également [R7.02.13], §2.1):

\(\beta =2\arctan\left[\frac{1}{4}\left({K}_{I}^{\mathit{èq}}/{K}_{\mathit{II}}^{\mathit{èq}}-\mathit{sign}({K}_{\mathit{II}}^{\mathit{èq}})\sqrt{{\left({K}_{I}^{\mathit{èq}}/{K}_{\mathit{II}}^{\mathit{èq}}\right)}^{2}+8}\right)\right]\)

Remarque – correction du champ \(\theta\) pour le rendre tangent à la surface fissurée:

Afin d’assurer \(\theta \cdot \nabla {\varphi}_{n}=0\) en tout point \(M\) de la surface fissurée (voir figures et ), le champ \(\theta\) est corrigé de la manière suivante (voir fig.):

- Initialisation \(\theta \leftarrow \nabla {\varphi}_{t}\) ;

- Vecteur \(\theta\) normal au plan de propagation: \(b\leftarrow \nabla {\varphi}_{n}\times \nabla {\varphi}_{t}\) , puis \(b\leftarrow b/\parallel b\parallel\) ;

- Correction pour assurer \(\theta\) tangent à la surface de fissure: \(\theta =\nabla {\varphi}_{n}\times b\) , puis \(\theta \leftarrow \theta /\parallel \theta \parallel\) *.*

Figure 6.4-b : Correction du champ θ pour le rendre tangent à la surface fissurée *.*

Remarque – champ \(\theta\) aux extrémité du front:

Dans une méthode \(\left(G,\theta \right)\) classique, une correction est appliquée aux valeurs de \(\theta\) sur les surfaces libres, pour assurer \(\theta \cdot \nu =0\) , où \(\nu\) est la normale à la surface libre (voir fig.), ceci afin d’éviter des valeurs de G négatives. Cette correction n’est plus appliquée dans le cas d’un calcul cohésif (voir fig.), car ce risque est exclu du fait de la nature du calcul (des intégrales de dissipation surfaciques).

Figure 6.4-c : Champ θ sur la surface fissurée pour un calcul (G, θ) classique, et pour un calcul cohésif.

Remarque – choix de la base \(\theta\) pour la décomposition des modes:

Si la zone cohésive a une taille importante, elle se répartit sur une surface incurvée. Le choix de la base pour la décomposition des modes n’est alors pas évident: on peut choisir:

-soit la base au front de fissure (figure , solution a),

-soit la base au point considéré (fig., solution b),

-soit la base liée aux coordonnées polaires au front de fissure (fig., solution c).

Tandis que la première solution sous-estime l’angle de bifurcation et la seconde le sur-estime, la troisième solution c) est celle qui produit les résultats les plus équilibrés, et donc celle qui a été implémentée. Le calcul de la base se déroule donc comme suit:

• Vecteur position dans le repère en fond de fissure: \(\mathit{PM}=x(M)-x(P)\) ;

-Projection dans le plan: \(t\leftarrow \mathit{PM}-(\mathit{PM}\cdot b)b\) , puis \(t\leftarrow t/\parallel t\parallel\) ;

-Construction d’un vecteur \(n\) normal: \(n\leftarrow t\times b\) , puis \(n\leftarrow n/\parallel n\parallel\) .

Figure 6.4-d : Choix de bases possibles pour la décomposition des modes.

Extension de la surface potentielle de fissuration

Pour l’actualisation de la surface fissurée, nous utilisons l’algorithme “GEOMETRIQUE” décrit en [R7.02.13], §6. Dans la phase de détection décrite en § 6.3 , seule la level-set tangente est actualisée (OPERATION=”DETECT_COHESIF” dans PROPA_FISS). Dans la phase d’extension de la surface potentielle de fissuration détaillée ici (OPERATION=”PROPA_COHESIF” dans PROPA_FISS), la level-set normale est aussi actualisée, et l’extension se fait selon une longueur fixe DA_MAX, qui doit être choisie très supérieure à \(\Delta l\) .

Les différentes méthodes d’actualisation de level-sets sont répertoriées dans la documentation [R7.02.13]. Tandis que les méthodes classiques se basent sur une discrétisation par différences centrées d’équations de Hamilton-Jacobi (METHODE_PROPA=”SIMPLEXE”/”UPWIND” dans PROPA_FISS), des méthodes plus simples sont apparues depuis (METHODE_PROPA=”GEOMETRIQUE” dans PROPA_FISS), pour lesquelles la seule étape restante est la mise à jour proprement dite. Cette technique tire parti du fait que la surface de fissure créée à l’instant précédent est figée: le problème de propagation d’une surface se réduit donc à la propagation d’un front à partir duquel des distances peuvent être directement évaluées.

Figure 6.5-a : Calcul direct, plan par plan, des nouvelles level-sets par l’algorithme de propagation géométrique.

A partir des équations d’évolutions, on arrive à montrer que pour chaque nœud \(M\) , pour peu que l’on arrive à projeter \(M\) sur un point \(P\) du front tel que \(M\in (P,{n}_{P},{t}_{P})\) , la mise à jour se ramène à un problème dans ce plan \((P,{t}_{P},{n}_{P})\) . On calcule donc les nouvelles valeurs, plan par plan , par translations et rotations planes, comme suit:

Pour chaque point M du maillage parent:

• On projette ce point sur l’ancien fond de fissure (fig., étape a).

• On en déduit la projection sur le nouveau front (fig., étape b).

• On calcule les level-sets correspondantes. Si on avait \({\varphi}_{t}\left(M\right)\le 0\) , on ne modifie pas \({\varphi}_{n}\left(M\right)\le 0\) , afin de «geler» la surface fissurée existante (fig., étape c).

Remarque – correction de la base \(({n}_{p},{t}_{p})\) aux points extrémités du front:

Une hypothèse importante de la méthode de propagation géométrique est que le point \(M\) appartienne au plan \((P,{n}_{p},{t}_{p})\) (pour rappel, \(M\) est le point de l’espace en lequel on cherche à actualiser les level-sets, et \(P\) est sa projection sur le fond de fissure). Si ce n’est pas le cas, il faut corriger la base \(({n}_{p},{t}_{p})\) pour assurer cette propriété avant d’appliquer l’algorithme de propagation.

Il y a trois cas de figure pour lesquels cette propriété n’est pas vérifiée, qui ont en commun le fait que le point \(M\) se projette sur une extrémité du fond de fissure. Par **ordre chronologiquede traitement:*

1. Le cas d’un front de fissure courbe

Dans ce cas, il est nécessaire de remplacer \({t}_{P}\) par un \({t'}_{P}\) tangent au bord de la structure: \({t'}_{P}\cdot \nu =0\) où \(\nu\) est la normale au bord de la structure.

Figure 6.5-b : Cas d’un front de fissure courbe.

Nous considèrerons que \({t}_{P}\) a été corrigé selon cette procédure par la suite: \({t'}_{P}\) sera simplement noté \({t}_{P}\) .

2. Le cas d’une surface fissurée inclinée par rapport au bord de la structure

Dans ce cas, il est nécessaire de corriger \({n}_{P}\) par un \({n}_{P}^{\text{*}}\) , de sorte que \(M\in (P,{n}_{P}^{\text{*}},{t}_{P})\) *. Il n’est pas nécessaire de corriger* \({t}_{P}\) . Cette correction s’écrit:

\({n}_{P}^{\text{*}}\leftarrow \left(1-{t}_{P}\otimes {t}_{P}\right)\cdot \mathit{PM}\) , et \({n}_{P}^{\text{*}}\leftarrow \frac{{n}_{P}^{\text{*}}}{\parallel {n}_{P}^{\text{*}}\parallel }\) .

Pour une fissure qui déboucherait perpendiculairement, on retrouve bien \({n}_{P}^{\text{*}}={n}_{P}\) .

Figure 6.5-c : Cas d’une surface inclinée par rapport au bord de la structure.

Nous considèrerons que \({n}_{P}\) a été corrigé selon cette procédure par la suite: \({n}_{P}^{\text{*}}\) sera simplement noté \({n}_{P}\) .

3. Le cas d’un bord de structure concave

Dans ce cas, il est nécessaire de corriger \({t}_{P}\) par un \({t}_{P}^{\text{*}}\) , de sorte que \(M\in (P,{n}_{P},{t}_{P}^{\text{*}})\) *. Il n’est pas nécessaire de corriger* \({n}_{P}\) . Cette correction s’écrit:

\({t}_{P}^{\text{*}}\leftarrow \left(1-{n}_{P}\otimes {n}_{P}\right)\cdot \mathit{PM}\) , et \({t}_{P}^{\text{*}}\leftarrow \frac{{t}_{P}^{\text{*}}}{\parallel {t}_{P}^{\text{*}}\parallel }\)

Figure 6.5-d : Cas d’un bord de structure concave.

Extension de l’espace de multiplicateurs et variables internes initiales

Comme illustré par la figure , étapes a-b, l’actualisation de la surface fissurée implique que cette dernière a tourné en amont du front. En conséquence, les arêtes qui sont intersectées par la nouvelle surface fissurée ne sont pas les mêmes que celles qui l’étaient par l’ancienne, en particulier en amont du front (voir fig., étapes a-b).

Sur la fig., étape d, nous notons \(Q\) l’ensemble des arêtes intersectées par la nouvelle surface fissurée, \(V\) le sous-ensemble des arêtes vitales (les arêtes qui portent une relation d’égalité), et \(K\) l’ensemble des nœuds enrichis. Nous notons \({Q}_{0}\) , \({V}_{0}\) et \({K}_{0}\) ces différents ensembles si la surface fissurée précédente est envisagée (fig., étape a). L’ensemble d’arêtes vitales \(V\) – qui définit le nouvel espace réduit de multiplicateurs \({M}_{h}\) – n’est pas un sous-ensemble de \({Q}_{0}\) mais de \(Q\) .

Une reconstruction de \(V\) en partant de zéro conduira probablement à des groupes d’arêtes vitales très différents. Comme les variables internes sont définies sur l’espace \({M}_{h}\) précédent, elles devraient être projetées entre deux espaces assez différents, ce qui amènerait une moins bonne conservation de l’énergie. Il est donc préférable de construire le nouvel espace en respectant les combinaisons existantes en deçà du front, et de l’étendre par de nouveaux groupes uniquement dans la zone non intersectée auparavant, située au delà du front de fissure (voir fig., étapes a-c).

On initialise l’ensemble d’arêtes donné à l’algorithme de restriction à l’ensemble \(Q\) des arêtes intersectées par la nouvelle surface fissurée. Ensuite, on enlève de cet ensemble (voir fig., étape b) les arêtes de \({Q}_{0}\setminus {V}_{0}\) dont les extrémités \(n\) sont toutes deux telles que:

• pour toute arête \(v\in {V}_{0}\) dont \(n\) est une extrémité, \(v\in Q\) .

L’algorithme de sélection des arêtes vitales détaillé dans la documentation [R5.03.54, §6] est ensuite mis en œuvre: par construction, il produira un ensemble \(V\) qui préserve les groupes de nœuds précédents (voir fig., étape c).

Le nouvel espace réduit implique des nœuds dans \(K\setminus {K}_{0}\) qui n’avaient pas de valeurs attribuées pour les variables internes (cercle vides sur la fig., étape c). Pour un tel noeud \(n\) , la valeur initiale est déterminée comme suit (voir fig., étapes c-d):

• Si il existe une arête de \(V\) connectant \(n\) à un noeud \(m\in {K}_{0}\) , la valeur de \(m\) est attribuée à \(n\) .

• Sinon, la variable interne initiale est mise à 0 (matériau sain).

Figure 6.6-a : Extension de l’espace de multiplicateurs et variables internes initiales .

Il y a également quelques nœuds de \({K}_{0}\) dont la valeur doit être changée, car ils cont connectés à un groupe différent entre \(V\) et \({V}_{0}\) (cercles gris sur la fig., étape c). Une fois de plus, si il existe une arête le connectant à un noeud \(m\in {K}_{0}\) , la valeur de \(m\) est attribuée à \(n\) . Notons que la présence de valeurs de variables internes non-nulles bien que faibles en amont du front (telles que \({\alpha}_{5}\) et \({\alpha}_{6}\) sur la figure , étape a) vient du fait que le front de propagation a été lissé lors de sa détection (voir § 6.3 ).

Bibliographie

1. 1. GRIFFITH,“The phenomenon of rupture and flow in solids», philosophical transactions of the royal society, A221, pp. 98-163, 1920

2. MOËS N., DOLBOW J., BELYTSCHKO T., “A finite element method for crack growth without remeshing”, International Journal for Numerical Methods in Engineering, vol. 46, pp. 135-150, 1999

3. Ji H., Dolbow J.E.,“On strategies for enforcing interfacial constraints and evaluating jumps conditions with the extended finite element method”, International Journal for Numerical Methods in Engineering, vol. 61, pp. 2508-2535, 2004

4. BOCCA P., CARPITIERI A., VALENTE S.“Mixed mode fracture of concrete”, International Journal of Solids and Structures, vol. 27, pp. 1139-1153, 1991

5. CENDON D., GALVEZ J., ELICES M., PLANAS J.“Modelling the fracture of concrete under mixed loading”, International Journal of Fracture, vol. 103, pp. 293-310, 2000

6. MESCHKE G., DUMSTORFF P., FLEMING W., JOX S. Computational failure analysis of concrete structures using the extended finite element method. Neue Bauweisen, trends in Statik und Dynamik, 2, 395 – 408, 2006.

7. ZAMANI A., GRACIE R., ESLAMI M. «Cohesive and non-cohesive fracture by higher order enrichment of XFEM». Int. J. Num. Meth. Engng, 90, 452 – 483, 2012.

8. MOËS N., BELYTSCHKO T. «Extended finite element method for cohesive crack growth». Engineering fracture mechanics, 69, 813 – 833, 2002.

9. ELICES M., PLANAS J.«Asymptotic analysis of a cohesive crack: 1.theoretical background». International Journal of Fracture, 55, 153 – 177, 1992.

10. MESCHKE G., DUMSTORFF P. Energy-based modelling of cohesive and cohesionless cracks via X-FEM. Computational Methods in Applied Mechanics Engineering, 196, 2338 – 2557, 2007.

11. FRIES T.P., BAYDOUN M. Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description. International Journal for Numerical Methods in Engineering, 89, 1527 – 1558, 2012.

12. ERDOGAN F., SIH G.C. On the crack entension in plane loading and transverse shear. Journal Basic Engng, 85, 519 – 527, 1963.

Description des versions

Indice document	Version Aster	Auteur(s) Organisme(s)	Description des modifications
A	12	G. FERTE, P.MASSIN EDF/R&D AMA	Texte initial