r5.03.54 Contact en petits glissements avec X-FEM

Résumé:

Ce document présente une nouvelle approche pour traiter les problèmes de contact en petits glissements avec la eXtended Finite Element Method (X-FEM) [bib1]. On considère la formulation hybride continue de problèmes de contact entre solides [bib2] et la stratégie de résolution est similaire à celle déjà implémentée dans Code_Aster pour le cadre éléments finis classiques [bib3]. Un nouveau type d’élément de contact mixte est introduit, spécifique au cadre X-FEM.

L’approche est implémentée dans Code_Aster en 2D et 3D, et traite à la fois des interfaces complètement coupées par une fissure ainsi que des interfaces avec fond de fissure. Elle est utilisable avec la commande STAT_NON_LINE [U4.51.03]. Le frottement de type Coulomb est pris en compte.

Formulation forte du problème de contact frottant

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

On désignera par \({r}^{1}\) et \({r}^{2}\) les densités des efforts dues aux interactions de contact frottant éventuelles entre les deux surfaces. On désignera par \({t}_{c}^{1}\) et \({t}_{c}^{2}\) les forces [1] dues aux interactions de cohésion éventuelles entre les deux surfaces dans le cas de l’ouverture d’une interface.

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

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

Remarques:

Malgré les apparences, les relations sur les densités d’effort de contact et les densités d’effort de cohésion ne sont pas incompatibles. En effet, la première n’est valide que lorsque les solides sont en contact alors que la seconde n’est applicable que si les solides sont décollés.

Lois du contact

Soit \(P\) un point de \({\Gamma}_{c}\) . On note \({P}^{1}\) et \({P}^{2}\) les points coïncidant sur \({\Gamma}^{1}\) et \({\Gamma}^{2}\) respectivement. La condition de non-interpénétration entre \({P}^{1}\) et \({P}^{2}\) est écrite dans la direction \(n\) , la normale à \({\Gamma}^{1}\) :

\({d}_{n}=(x({P}^{1})-x({P}^{2}))\cdot n\le 0\)

Figure 2.2-1 : Définition du jeu

On décompose la densité d’effort de contact \(r\) en une partie normale \(\lambda\) qui désigne la pression normale de contact et une autre tangentielle \({r}_{\tau}\) . Ainsi, la densité d’effort de contact s’écrit:

\(r=\mathrm{\lambda n}+{r}_{\tau}\)

Figure 2.2-2 : Définition de la densité d’effort de contact

Avec ces notations, les lois de contact (lois de Signorini) s’écrivent sous la forme suivante:

\({d}_{n}\le 0,\lambda \le 0,\lambda {d}_{n}=0\)

Ces lois font intervenir des inéquations, or celles-ci ne se prêtent pas aisément à une formulation faible. Pour cela, on réécrit ces lois sous une autre forme, en les transformant en une seule équation équivalente [bib 1 ]:

\(\lambda -\chi ({g}_{n}){g}_{n}=0\)

Dans cette expression, \(\chi\) est la fonction indicatrice de \({\mathrm{\Re }}^{-}\) définie par

\(\chi (x)=\lbrace \begin{array}{c}1\text{si}x<0\\ 0\text{si}x\ge 0\end{array}\)

Figure 2.2-3 : Définition de la fonction indicatrice de \({\Re }^{-}\)

et \({g}_{n}\) le multiplicateur (dit de contact augmenté [bib 3 ]) défini par:

\({g}_{n}=\lambda -{\rho}_{n}{d}_{n}\)

où \({\rho}_{n}\) est un réel strictement positif.

Le problème de contact posé par les lois de Signorini introduit une relation non-univoque (\(\lambda\) n’est pas une fonction de \({d}_{n}\) ), semi-définie positive et non différentiable en \(\lambda ={d}_{n}=0\) comme l’illustre la [].

Figure 2.2-4 : Graphe de la loi de contact unilatéral de Signorini

Une régularisation par pénalisation de l’interpénétration permet de rendre cette relation univoque, voir [].

Figure 2.2-5 : Graphe de la loi de contact régularisée

Physiquement, on autorise les solides à s’interpénétrer au prix d’une raideur \({\kappa}_{n}\) très élevée. La pression de contact est alors donnée par \(\lambda =-{\kappa}_{n}{d}_{n}\) et devient égale au multiplicateur de contact \({g}_{n}\) . On note que plus on augmente le coefficient de pénalisation \({\kappa}_{n}\) , plus on se rapproche de la loi de contact classique.

Lois du frottement

Pour les phénomènes de frottement, on utilise les lois de Coulomb qui s’écrivent comme suit:

\(\begin{array}{c}\mathrm{\parallel }{r}_{\tau}\mathrm{\parallel }\le \mu \mid \lambda \mid \\ \text{Si}\mathrm{\parallel }{r}_{\tau}\mathrm{\parallel }<\mu \mid \lambda \mid \text{alors}{v}_{\tau}=0\\ \text{Si}\mathrm{\parallel }{r}_{\tau}\mathrm{\parallel }=\mu \mid \lambda \mid \text{alors}\exists \alpha \ge 0;{v}_{\tau}=-\alpha {r}_{\tau}\end{array}\)

où \(\mu\) est le coefficient de frottement de Coulomb et \({v}_{\tau}\) la vitesse relative tangente.

On note par la suite \({x}_{\tau}\) la projection de \(x\) sur le plan tangent à la surface de contact, définie par \({x}_{\tau}=(\text{Id}-n\otimes n)x\) , où le symbole \(\otimes\) désigne le produit tensoriel.

De même que pour les lois de contact, on peut écrire les lois de frottement de manière équivalente comme suit:

\({r}_{\tau}=\mu \lambda \Lambda\)

\(\Lambda -{{\rm P}}_{{\rm B}(0,1)}({g}_{\tau})=0\)

\({g}_{\tau}=\Lambda +{\rho}_{\tau}{v}_{\tau}\)

Dans ces expressions, \(\Lambda\) est un semi-multiplicateur (vectoriel) de frottement, \({g}_{\tau}\) est le semi-multiplicateur (vectoriel) de frottement augmenté, \({P}_{B(0,1)}\) est la projection sur la boule unité et \({\rho}_{\tau}\) un paramètre strictement positif. Le semi-multiplicateur de frottement \(\Lambda\) , dont le module est toujours inférieur ou égal à 1, correspond à la direction de glissement quand son module vaut 1, et correspond à la direction d’adhérence quand son module est strictement inférieur à 1. La [ Figure 2.3-1 ] présente en 3D un cas de glissement et un cas d’adhérence.

Figure 2.3-1 : Projections sur la boule unité en 2D: glissement (à gauche) et adhérence (à droite)

Les lois de frottement sont complétées par l’équation (de type exclusion) suivante:

\({d}_{n}\Lambda =0\text{ou}(1-\chi ({g}_{n}))\Lambda =0\)

Le problème de frottement exprimé via les lois de Coulomb introduit une relation non-univoque (\({R}_{\tau}\) n’est pas une fonction de \({v}_{\tau}\) ), et non différentiable (voir []).

Figure 2.3-2 : Graphe de la loi de frottement de Coulomb

La méthode pénalisée permet de rendre cette relation univoque (voir []).

Figure 2.3-3 : Graphe de la loi de frottement régularisée

Le semi-multiplicateur de frottement vectoriel est alors donné par le relation \(\Lambda ={P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau})\) .

On note que plus on augmente le coefficient de pénalisation \({\kappa}_{\tau}\) , plus on se rapproche de la loi de contact classique.

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 disponible dans Code_Aster, les lois CZM_EXP_REG et CZM_LIN_REG dont les caractéristiques sont détaillées dans [R7.02.11]. Nous détaillons ici l’extension à XFEM pour la loi CZM_EXP_REG, en se basant sur l’article de référence [bib 17 ]. L’extension à la loi CZM_LIN_REG est faite en suivant exactement le même paradigme.

Nous écrivons le saut de déplacement tel que défini pour l’algorithme de contact, avec la notation définie par la figure :

\(\mathrm{〚}u\mathrm{〛}({P}^{1})=u({P}^{1})-u({P}^{2})\)

Nous avons alors \(\mathrm{〚}{u}_{n}\mathrm{〛}=\mathrm{〚}u\mathrm{〛}\cdot n\) négatif en ouverture et positif en interpénétration. Les notations de [R7.02.11] posent un saut de déplacement \(\delta\) tel que \(\delta \cdot n\) soit positif en ouverture et négatif en interpénétration. Pour nous ramener à nos notations, nous réutilisons les résultats de [R7.02.11] en posant \(\delta =-\mathrm{〚}u\mathrm{〛}\) . Nous nous contentons ici d’en rappeler les principaux points. Le lecteur peut se s’y référer pour une compréhension approfondie.

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\)

• \({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}=\sqrt{{\langle \mathrm{〚}{u}_{n}\mathrm{〛}\rangle }_{-}^{2}+{\mathrm{〚}{u}_{\tau}\mathrm{〛}}^{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}_{\mathrm{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 )\mathrm{〚}u\mathrm{〛}\) 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{0}{\int}}{t}_{c,\tau }\cdot d\mathrm{〚}{u}_{\tau}\mathrm{〛}+\underset{-\infty }{\overset{0}{\int}}({t}_{c,n}\cdot n)d\mathrm{〚}{u}_{n}\mathrm{〛}={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.4.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{\mathrm{〚}u\mathrm{〛}}}_{\mathit{eq}}={t}_{c,n}\cdot n\dot{\mathrm{〚}{u}_{n}\mathrm{〛}}+{t}_{c,\tau }\cdot \dot{\mathrm{〚}{u}_{\tau}\mathrm{〛}}\)

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}\cdot n\frac{{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}}{\mathrm{〚}{u}_{n}\mathrm{〛}}=\mathrm{\parallel }{t}_{c,\tau }\mathrm{\parallel }\frac{{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}}{\mathrm{\parallel }\mathrm{〚}{u}_{\tau}\mathrm{〛}\mathrm{\parallel }}=-\frac{{\sigma}_{c}}{\alpha}\exp(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha ){\mathrm{〚}u\mathrm{〛}}_{\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.4.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(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha ){\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\) On définit parfois un saut de déplacement équivalent \({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}=\sqrt{{\langle \mathrm{〚}{u}_{n}\mathrm{〛}\rangle }_{-}^{2}+{\beta}^{2}{\mathrm{〚}{u}_{\tau}\mathrm{〛}}^{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,\mathrm{eq}}〚{u}_{n}〛}{{〚u〛}_{\mathrm{èq}}}n=-\frac{{\sigma}_{c}}{\alpha}\exp(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha )〚{u}_{n}〛n\) , \({t}_{c,\tau }={\beta}^{2}\frac{{t}_{c,\mathit{eq}}\mathrm{〚}{u}_{\tau}\mathrm{〛}}{{\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}}=-{\beta}^{2}\frac{{\sigma}_{c}}{\alpha}\exp(-\frac{{\sigma}_{c}}{{G}_{c}}\alpha )\mathrm{〚}{u}_{\tau}\mathrm{〛}\)

Lois cohésives mixtes

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 .

Comme en partie 2.4.1 , nous introduisons le saut de déplacement \(\mathrm{〚}u\mathrm{〛}\) tel que \(\mathrm{〚}{u}_{n}\mathrm{〛}\) soit négatif en interpénétration et positif en ouverture, et nous nous ramenons aux notations de [R7.02.11] en posant \(\delta =-\mathrm{〚}u\mathrm{〛}\) : on renvoie à cette documentation pour l’expression de la loi et sa dérivée. Notons que de la même façon que pour les lois régularisées, le matériau reste dans le domaine élastique tant que:

\(f({\mathrm{〚}u\mathrm{〛}}_{n},\alpha )=-{\mathrm{〚}u\mathrm{〛}}_{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.4.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).

[1]

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

Formulation variationnelle mixte

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\)

Commençons alors par donner une formulation unifiée commune aux cas de contact-frottement et aux lois cohésives régularisées . Les lois cohésives mixtes obéissent quand à elles à une formulation énergétique suivant une logique différente, expliquée en partie 3.5 .

Formalisme commun au cohésif régularisé et au contact-frottement

Pour ce faire nous noterons \(r={t}_{c}\) dans le cas des lois régularisées. Nous notons \(H={H}^{-1/2}(\Gamma )\) pour les lois cohésives, et nous désignons par \(H\) le sous-espace de \({H}^{-1/2}(\Gamma )\) des champs d’actions de contact dans le cas du contact-frottement. La formulation faible du problème de contact frottant s’écrit comme suit:

Trouver \((u,{r}^{1},{r}^{2})\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}^{1}}{r}^{1}\cdot {u}^{{\text{*}}^{1}}{\mathit{d\Gamma }}^{1}+{\int}_{{\Gamma}^{2}}{r}^{2}\cdot {u}^{{\text{*}}^{2}}{\mathit{d\Gamma }}^{2}\forall {u}^{\text{*}}\in {V}_{0}\)

En écrivant le saut avec la notation suivante

\([[x]]({P}^{1})=x({P}^{1})-x({P}^{2})\) ,

et en notant \(r={r}^{1}\) et par la suite, la formulation faible du problème de contact frottant s’écrit de manière équivalente comme suit:

\({\int}_{\Omega}\sigma (u):\varepsilon ({u}^{\text{*}})d\Omega ={\int}_{\Omega}f\cdot {u}^{\text{*}}d\Omega +{\int}_{{\Gamma}_{t}}t\cdot {u}^{\text{*}}d\Gamma +{\int}_{{\Gamma}_{c}}r\cdot [[{u}^{\text{*}}]]d{\Gamma}_{c}\forall {u}^{\text{*}}\in {V}_{0}\)

Les espaces des inconnues de contact sont les suivants:

\(\begin{array}{c}H=\left\lbrace \lambda \in {H}^{-1/2}({\Gamma}_{c}),\lambda \le 0\text{sur}{\Gamma}_{c}\right\rbrace \\ H=\left\lbrace {r}_{\tau}\in {H}^{-1/2}({\Gamma}_{c}),\mathrm{\parallel }{r}_{\tau}\mathrm{\parallel }\le \mu \lambda \text{sur}{\Gamma}_{c}\right\rbrace \end{array}\)

Méthode du Lagrangien augmenté

La formulation faible à trois champs s’écrit finalement,dans le cas d’une formation Lagrangienne augmentée :

Trouver \((u,\lambda ,\Lambda )\in {V}_{0}\times H\times H\)

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}},{\Lambda}^{\text{*}})\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}}\chi ({g}_{n}){g}_{n}n\cdot [[{u}^{\text{*}}]]{\mathit{d\Gamma }}_{c}-{\int}_{{\Gamma}_{c}}\chi ({g}_{n}){\text{μλ}}_{s}{{\rm P}}_{{\rm B}(0,1)}({g}_{\tau})\cdot [[{u}_{\tau}^{\text{*}}]]{\mathit{d\Gamma }}_{c}=0\end{array}\)
Loi de contact	\({\int}_{\mathit{\Gamma c}}\frac{-1}{{\rho}_{n}}(\lambda -\chi ({g}_{n}){g}_{n}){\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\)
Loi de frottement	\(\begin{array}{c}{\int}_{\mathit{\Gamma c}}\frac{\text{μχ}({g}_{n}){\lambda}_{s}\mathit{\Delta t}}{{\rho}_{\tau}}(\Lambda -{P}_{B(0,1)}({g}_{\tau})){\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ +{\int}_{\mathit{\Gamma c}}(1-\chi ({g}_{n}))\Lambda {\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\end{array}\)

Méthode pénalisée

La formulation faible à trois champs s’écrit finalement, dans le cas d’une formation purement pénalisée :

Trouver \((u,\lambda ,\Lambda )\in {V}_{0}\times H\times H\)

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}},{\Lambda}^{\text{*}})\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}}\chi ({g}_{n})\lambda \cdot [[{u}^{\text{*}}]]\cdot n{\mathit{d\Gamma }}_{c}\\ -{\int}_{{\Gamma}_{c}}\chi ({g}_{n}){\text{μλ}}_{s}{{\rm P}}_{{\rm B}(0,1)}({\kappa}_{\tau}{\nu}_{\tau})\cdot [[{u}_{\tau}^{\text{*}}]]{\mathit{d\Gamma }}_{c}=0\end{array}\)
Loi de contact	\({\int}_{\mathit{\Gamma c}}\frac{-1}{{\kappa}_{n}}(\lambda +\chi ({g}_{n}){\kappa}_{n}{d}_{n}){\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\)
Loi de frottement	\(\begin{array}{c}{\int}_{\mathit{\Gamma c}}\frac{\text{μχ}({g}_{n}){\lambda}_{s}}{{\kappa}_{\tau}}(\Lambda -{P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau})){\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ +{\int}_{\mathit{\Gamma c}}(1-\chi ({g}_{n}))\Lambda {\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}=0\end{array}\)

On remarquera que la valeur de \(\lambda\) obtenue par la loi de contact est une moyenne lorsque l’état est contactant des forces d’interpénétration pour la loi pénalisée. L’utilisation des expressions \(\lambda\) ou \(-{\kappa}_{n}{d}_{n}\) dans l’équation d’équilibre devrait donc conduire au même résultat si ce n’est que la condition LBB ne s’applique que sur les \(\lambda\) . Pour avoir équivalence, il faudrait donc reporter le traitement de la LBB sur les champs de déplacement pour la méthode pénalisée. Comme cela n’est pas fait ici, on utilise le terme \(\lambda\) pour lequel le traitement LBB est fait et on le réinjecte dans l’équation d’équilibre qui prend donc en compte ce traitement. On pourrait tenter le même traitement pour la loi de frottement mais \(\Lambda\) est obtenue comme une moyenne pour des situations glissantes ou adhérentes. Réinjecter cet état moyenné au niveau de l’équation d’équilibre entraîne une non convergence de l’algorithme de Newton. La différence de comportement pour \(\lambda\) et \(\Lambda\) provient du fait que l’on intègre des quantités discontinues suivant l’état contactant ou non contactant pour \(\lambda\) et suivant l’état glissant ou adhérent pour \(\Lambda\) , mais que pour l’état non contactant les contributions sur \(\lambda\) sont nulles dans la loi de contact. On peut en déduire que s’il existe une formulation pénalisée satisfaisant pleinement les conditions LBB pour le contact, il n’existe pour le moment pas de formulation pénalisée satisfaisant pleinement les conditions LBB pour le frottement. Le choix du coefficient \({\kappa}_{\tau}\) en sera du coup d’autant plus important (de fortes valeurs risquant d’exhiber des blocages pour la partie adhérente).

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

Lorsqu’une loi cohésive est utilisée, le contact est géré par le coefficient de pénalisation défini dans la loi cohésive. Connaissant l’expression de la loi cohésive en fonction de \(\mathrm{〚}u\mathrm{〛}\) , la formulation s’écrit:

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

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

Équation d’équilibre	:math:`begin{array}{c}{int}_{Omega}sigma (u):epsilon ({u}^{text{*}})dOmega -{int}_{Omega}fcdot {u}^{text{*}}dOmega -{int}_{{Gamma}_{t}}tcdot {u}^{text{*}}mathit{dGamma }\ -{int}_{{Gamma}_{c}}{t}_{c	n}cdot {mathrm{〚}{u}^{text{*}}mathrm{〛}}_{n}{mathit{dGamma }}_{c}-{int}_{{Gamma}_{c}}{t}_{c	tau }cdot {mathrm{〚}{u}^{text{*}}mathrm{〛}}_{tau}{mathit{dGamma }}_{c}=0end{array}`
Post-traitement partie normale	:math:`{int}_{mathit{Gamma c}}(lambda -{t}_{c	n}cdot n){lambda}^{text{*}}{mathit{dGamma }}_{c}=0`
Post-traitement partie tangentielle	:math:`{int}_{mathit{Gamma c}}(Lambda -{t}_{c	tau }){Lambda}^{text{*}}{mathit{dGamma }}_{c}=0`

On remarque que les multiplicateurs \(\lambda\) et \(\Lambda\) n’interviennent pas dans la résolution. Ils servent uniquement à stocker les contraintes cohésives de manière explicite.

Formulation pour une loi cohésive mixte

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ù un champ dual vectoriel \(\lambda\) va effectivement entrer dans la formulation, au lieu d’être un artifice de post-traitement comme en 3.4 . Cette formulation suit un raisonnement énergétique, expliqué en détail dans la documentation [R3.06.13]. Résumons-en en 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, avec nos conventions de signe, à l’opposé du saut de déplacement. On cherche:

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

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}(u,\delta ,\lambda )\underset{\mathit{déf.}}{=}E(u,\delta )+\underset{\Gamma}{\int}\lambda \cdot (-\mathrm{〚}u\mathrm{〛}-\delta )\mathit{d\Gamma }+\underset{\Gamma}{\int}\frac{r}{2}{(\mathrm{〚}u\mathrm{〛}+\delta )}^{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-\lambda +r(\mathrm{〚}u\mathrm{〛}+\delta )\right]\cdot {\delta}^{\text{*}}=0\text{avec}t\in \partial \Pi (\delta )\)

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. Cette premiè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\) est discré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 conditions d’optimalité revient à satisfaire la loi cohésive en chaque point de collocation:

\({t}_{c}({\delta}_{g},{\alpha}_{g})={\lambda}_{g}-r(\mathrm{〚}{u}_{g}\mathrm{〛}+{\delta}_{g})\)

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\mathrm{〚}{u}_{g}\mathrm{〛}-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.5-a: Solution de l’intégration du comportement.

Le champ \(\delta\) s’écrit donc localement comme une fonction de \(\lambda -r\mathrm{〚}u\mathrm{〛}\) , 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\)

..csv-table:

Équation d’équilibre, :math:`\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(\mathrm{〚}u\mathrm{〛}+\delta (p))\right]\cdot \mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma =0\text{avec}p=\lambda -r\mathrm{〚}u\mathrm{〛}\end{array}`
Loi d'interface     , :math:`-{\int}_{{\Gamma}_{c}}\left[\mathrm{〚}u\mathrm{〛}+\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 le déplacement avec des éléments P2 et le multiplicateur d’une façon P1 adaptée à X-FEM, dont nous donnons le détail dans la partie 4.1 à suivre.

Discrétisations EF

On détaille dans cette partie la discrétisation des inconnues de contact-frottement, qui se fait aux nœuds sommets de l’élément parent.

Multiplicateurs de contact

Les inconnues pour la pression de contact \(\lambda\) et le semi-multiplicateur de frottement \(\Lambda\) sont portées aux nœuds sommets de l’élément parent. L’approximation de la pression de contact fait intervenir les \({\phi }_{i}\) , fonction de forme de l’élément parent linéaire et s’écrit :

\({\lambda}^{h}(x)=\sum_{i=1}^{\mathrm{nno}}{\lambda}_{i}{\phi }_{i}(x)\)

où \(\mathrm{nno}\) est le nombre de nœuds de l’élément parent linéaire. La pression est alors définie comme la trace sur l’interface approximée de ce champ \({\lambda}^{h}\) .

Pour expliciter cette notion d’interface approximée, supposons que l’on découpe, si besoin, l’élément parent en sous-éléments simpliciaux (c’est-à-dire triangles en 2D, tétraèdres en 3D). Par interface approximée, nous entendons en 2D la ligne brisée reliant entre eux les points d’intersection des arêtes de tels sous-éléments avec la ligne de la fissure. En 3D, les points d’intersection des arêtes d’un tel sous-éléments avec la surface de la fissure définissent un polygone, non nécessairement plan au sein de ce sous-élément: il peut y par exemple avoir 4 points d’intersection, pas forcément coplanaires. La méthode retenue a été celle qui consiste à découper ce polygone en facettes triangulaires dont les sommets sont ces points d’intersection ([fig. ]). L’ensemble du procédé décrit dans ce paragraphe est détaillé en partie 4.3 .

Semi-multiplicateurs de frottement

De même que pour les multiplicateurs de contact, les semi-multiplicateurs de frottement sont interpolés avec les \({\phi }_{i}\) fonction de forme de l’élément parent. En 3D, \(\Lambda\) est un vecteur du plan tangent à la surface de la fissure. Les gradients des level sets permettent de définir une base covariante à la surface de la fissure, dans laquelle \(\Lambda\) sera exprimé. On définit les 2 vecteurs de la base covariante par:

\(({n}^{\mathrm{ls}},{\tau}^{1},{\tau}^{2})=(\nabla \mathrm{lsn},\nabla \mathrm{lst},\nabla \mathrm{lsn}\times \nabla \mathrm{lst})\)

où \({n}^{\mathrm{ls}}\) est la normale locale issue du gradient de la level set normale. Les vecteurs \({\tau}^{1},{\tau}^{2}\) étant issus des gradients (nodaux) des level sets, ils peuvent être interpolés au sein des éléments de manière à obtenir des vecteurs aux sommets \(i\) des facettes de contacts, soit \({\tau}_{i}^{1}\) et \({\tau}_{i}^{2}\) , \(i=1,3\) . L’approximation des semi-multiplicateurs de frottement sur une facette de contact s’écrit alors :

\({\Lambda}^{h}(x)=\sum_{i=1}^{\mathrm{nno}}({\Lambda}_{i}^{1}{\tau}_{i}^{1}+{\Lambda}_{i}^{2}{\tau}_{i}^{2}){\phi }_{i}(x)\) ,

où les vecteurs \({\tau}_{i}^{1},{\tau}_{i}^{2}\) pour un nœud i correspondent à ceux du point d’intersection de l’interface avec l’arête de l’élément associée au nœud i. Le point d’intersection associé est soit ce nœud (s’il s’agit d’un nœud par où passe l’interface de contact), soit l’arête intersectée contenant ce nœud (s’il s’agit d’une arête coupée). Dans le cas où le nœud appartient à plusieurs arêtes intersectées, on l’associe au point d’intersection de l’arête vitale.

Élément fini de contact

Cas général

Les degrés de liberté de contact sont portés exclusivement par les nœuds sommets.

A titre d’exemple, la figure montre les nœuds portants les degrés de liberté de contact. On remarque que l’on distingue d’une part les éléments X-FEM coupés par la fissure qui porteront des degrés de liberté de contact et d’autre part les éléments X-FEM non coupés qui n’ont pas besoin de degrés de liberté de contact. On introduit des degrés de liberté de contact uniquement sur les nœuds sommets soit au total 4 degrés de liberté. Deux relations d’égalités relient respectivement les nœuds 1 et 3 ainsi que les nœuds 2 et 4 afin de satisfaire la condition LBB (voir [§ 6 ]). Cela fait un total de 6 variables introduites.

Figure 4.3.1-1 : Nœuds portant les ddl de contact.

FACETTE DE CONTACT

La surface de contact sert seulement à des fins d’intégration mais elle nécessite la mise en place d’algorithme de recherche des points d’intersection et des points milieux si l’élément est quadratique.

L’algorithme de recherche des points d’intersection se présente de la manière suivante :

Sur un élément:

• boucle sur les arêtes de l’élément

soient \({E}^{1}\) et \({E}^{2}\) les deux extrémités de l’arête

si \(\text{lsn}({E}^{1})\text{lsn}({E}^{2})\le 0\) alors

• • boucle sur les deux extrémités

si \(\text{lsn}({E}^{j})=0\) et \(\text{lst}({E}^{j})\le 0\) alors

on ajoute le point \({E}^{j}\) à la liste des \({P}_{i}\) (avec vérification des doublons)

fin si

si élément quadratique 2D : \(\text{lsn}({E}^{3})=0\) et \(\text{lst}({E}^{3})\le 0\) alors

on ajoute le point \({E}^{3}\) à la liste des \({P}_{i}\) (avec vérification des doublons)

fin si

• • fin boucle

si \(\text{lsn}({E}^{k})\ne 0\forall k\in\) nombre de nœuds de l’arête, alors

interpolation des coordonnées de C

si \(\text{lst}(C)\le 0\) alors

on ajoute le point \(C\) à la liste des \({P}_{i}\) (avec vérification des doublons)

fin si

fin si

fin si

• fin boucle

Détails de l’interpolation des coordonnées du point C :

Si l’élément est linéaire :

\(s=\frac{\text{lsn}({E}^{1})}{\text{lsn}({E}^{1})-\text{lsn}({E}^{2})}\)

\(C={E}^{1}-s({E}^{2}-{E}^{1})\)

\(\text{lst}(C)=\text{lst}({E}^{1})-s(\text{lst}({E}^{2})-\text{lst}({E}^{1}))\)

Si l’élément est quadratique 2D :

La position du point sur l’arête renseigne sur la valeur d’une de ses coordonnées de référence \(\xi\) ou \(\eta\) , il suffit ensuite de résoudre l’équation polynomiale \(\sum_{i=1}^{\mathrm{nno}}{\Phi}_{i}(\xi ,\eta ){\mathrm{lsn}}_{i}=0\) pour trouver la valeur de la seconde coordonnées de référence. Par passage dans l’élément réel, on détermine ses coordonnées réelles \((x,y)\) .

fin si

Dans le cas où l’élément est quadratique, il faut, en plus des points d’intersection interpoler les coordonnées des nœuds milieux de la facette de contact. En 2D, la facette de contact s’apparente à un arc de cercle où l’on ne connaît que les coordonnées de ses extrémités suite à l’algorithme précédent. Pour déterminer son point milieu, on utilise la méthode de Newton qui évalue le point d’intersection entre la médiatrice du segment reliant les deux extrémités et l’izo- zéro de la level set normale.

Cas des éléments en fond de fissure

Pour un élément contenant le fond de fissure, il faut un traitement particulier pour déterminer les points de contact. En effet, de tels éléments ne sont pas entièrement coupés par la fissure. Les points de contact sont alors de deux types:

• soit des intersections entre la surface de la fissure et les arêtes de l’élément (cas général évoqué au paragraphe précédent),

• soit des intersections entre le fond de fissure et les faces de l’élément (cas spécifique aux éléments contenant le fond de fissure).

Les points du 1er type sont déterminés par l’algorithme précédent, et les points du 2nd type par l’algorithme de recherche des points du fond de fissure (voir le paragraphe [§2.4] dans [R7.02.12]).

Les points de contact du 2nd type ne sont associés à aucun nœud ni arête, et ne sont donc portés nativement par aucun ddl. Ils n’interviennent donc pas dans l’écriture de l’approximation. Cette situation correspond à une approximation P1 des inconnues de contact sur la facette de contact en fond, ayant une valeur nulle en fond de fissure (voir schéma du haut de la Figure 4.3.2-1 ). Une autre solution est d’envisager un raccord constant. Dans ce cas, la pression en fond de fissure est considérée égale à la pression du point de contact sur une arête le plus proche. C’est cette solution qui a été retenue. Dans l’exemple de la Figure 4.3.2-1 , le schéma du bas illustre l’approximation de la pression de contact sur la facette \(\text{ABC}\) :

\(\begin{array}{ccc}\lambda (x)& =& {\lambda}_{A}{\psi}_{A}+{\lambda}_{B}{\psi}_{B}+{\lambda}_{C}{\psi}_{C}\\ & =& {\lambda}_{D}{\psi}_{A}+{\lambda}_{C}{\psi}_{B}+{\lambda}_{C}{\psi}_{C}\end{array}\)

Puisque \(\text{BC}<\text{BD}\) et \(\text{AD}<\text{AC}\) , l’inconnue de pression en \(B\) est portée par le point \(C\) et l’inconnue de pression en \(A\) est portée par le point \(D\) , il n’y a donc pas de ddl supplémentaire sur les faces. Cette approche pourrait être qualifiée de «P0-P1», car l’approximation est un mélange de P0 et de P1. La pression de contact est P1 le long du fond de fissure, et P0 le long des autres cotés de la facette.

Figure 4.3.2-1 : Différentes approximations du contact sur la facette en fond de fissure

Une autre alternative serait de considérer une vraie approximation P1 sur les facettes de contact en fond de fissure. Pour cela, il faudrait que les degrés de liberté de contact sur les points du fond soient de vrais degrés de liberté indépendants. Ils pourraient être par exemple portés par les nœuds sommets des arêtes opposées. Sur l’exemple de la Figure 4.3.2-2 , la pression de contact en \(B\) serait portée par le nœud \(F\) et la pression de contact en \(A\) par le nœud \(E\) . Ce cas est généralisable à tout type de configuration. L’intérêt est de disposer de degrés de liberté de contact indépendant. Une telle approximation P1 améliorerait la précision par rapport à une approximation constante.

Figure 4.3.2-2 : Approximation P1 sur la facette de contact en fond de fissure

Sous-découpage en facettes triangulaires de contact

Pour chaque élément, à partir de la liste des points de contact \({P}_{i}\) , il faut créer un sous-découpage en facettes triangulaires. Pour cela, on trie les points de contact par élément avec le même procédé que celui servant à orienter les points du fond fissure (voir le paragraphe [§2.4] dans [R7.02.12]). On détermine \(n\) la moyenne des normales aux nœuds (basées sur les gradients des level sets). On détermine \(G\) le barycentre des points de contact sur cet élément. On projette les points de contact dans le plan de normale \(n\) et passant par \(G\) . Pour chacun des projetés, on détermine l’angle \(\theta\) par rapport au 1er point de la liste, puis on trie les points de la liste suivant \(\theta\) croissant.

Pour illustrer le sous-découpage en facettes triangulaires, prenons un hexaèdre occupant la région \({\left[0,1\right]}^{3}\) . Soit le plan d’équation cartésienne \(\mathrm{4x}-11y-\mathrm{9z}+d=0\) . Examinons les intersections entre ce plan et l’hexaèdre, pour différentes valeurs du paramètre \(d\) .

Pour \(d=-1\) , il y 3 points d’intersection entre le plan et l’hexaèdre. La trace du plan dans l’hexaèdre est un triangle, qui correspond à la facette de contact. Pour \(d=4\) , il y 4 points d’intersection entre le plan et l’hexaèdre. La trace du plan dans l’hexaèdre est un quadrangle, est découpé en deux facettes triangulaires de contact. Pour \(d=6\) , il y 5 points d’intersection entre le plan et l’hexaèdre. La trace du plan dans l’hexaèdre est un pentagone, est découpé en trois facettes triangulaires de contact. Pour \(d=8\) , il y 6 points d’intersection entre le plan et l’hexaèdre. La trace du plan dans l’hexaèdre est un hexagone, est découpé en quatre facettes triangulaires de contact. La Figure 4.3.3-1 présente les divers schémas pour les valeurs de \(d\) précédemment évoquées.

De plus, lorsque le fond de fissure est contenu dans un élément, il se peut que cela rajoute un point d’intersection en plus. Par exemple pour \(d=8\) , si le fond de fissure coupe les segments P1P2 et P2P3 alors cela rajoute un point P2a situé sur P1P2 et un autre point P2b situé sur P2 P3 et supprime le point P2 (2 points ajoutés et 1 point supprimé, soit au total 1 point rajouté). Le nombre maximal de points d’intersection est alors 7.

Cet exemple illustre les différents cas de figure qui peuvent se produire. De manière générale, on peut regrouper les différents cas de figures suivant le nombre de points d’intersection rencontrés. Le Tableau 4.3.3-1 rassemble les sous-découpages effectués selon de le nombre de points d’intersection trouvés entre l’élément (que l’élément soit un tétraèdre, un pentaèdre ou un hexaèdre) et la surface de la fissure.

	3 points d’intersection	4 points d’intersection	5 points d’intersection	6 points d’intersection	7 points d’intersection
Facettes triangulaires	P1 P2 P3	P1 P2 P3 P1 P3 P4	P1 P2 P3 P1 P3 P4 P1 P4 P5	P1 P2 P3 P1 P3 P5 P1 P5 P6 P3 P4 P5	P1 P2 P3 P1 P3 P5 P3 P4 P5 P1 P5 P7 P5 P6 P7

Tableau 4.3.3-1 : Découpage en facettes triangulaires

Figure 4.3.3-1 : Intersections et découpages pour d = -1, 4, 6 et 8 (de haut en bas)

Mise à zéro des degrés de liberté inactifs

Cette procédure sert à mettre à zéro la valeur des degrés de liberté de qui n’interviennent pas dans les équations. Plusieurs solutions peuvent être considérées.

Remarque :

La discussion qui suit est fortement influencée par les possibilités (et les restrictions) de Code_Aster; les solutions envisagées ne sauraient alors être exhaustives.

Ne pas introduire

Tout d’abord, ne pourrait-on pas simplement éviter d’introduire ces degrés de liberté qui ne servent à rien et qu’il faudra annuler par la suite? On pourrait par exemple imaginer un élément initialement sans degré de liberté de contact, et une fois les intersections entre les arêtes de l’élément et la fissure déterminées, rajouter les degrés de liberté de contact nécessaires. Cela voudrait dire transformer l’élément fini en élément fini possédant des degrés de liberté de contact en certains nœuds précis (par exemple aux nœuds \(\mathit{N1}\) , \(B\) , \(C\) , \(\mathit{N4}\) de la fig. ). On entrevoit aisément le nombre de possibilités de cela représente, et le nombre d’éléments finis différents que cela impliquerait (21 pour un tétraèdre avec 3 ou 4 points d’intersection, et plus de 150 pour un hexaèdre avec 3, 4, 5 ou 6 points d’intersection). Cette solution est donc écartée. On est donc bien obligé d’avoir un élément de contact avec des degrés de liberté de contact sur tous les nœuds sommet et milieu.

Élimination

Pour éliminer les degrés de liberté de contact qui ne servent à rien, l’idée la plus simple est celle qui consiste à retirer du système global d’équations les lignes et les colonnes correspondantes (sur le papier, cela revient à «rayer» les lignes et les colonnes inutiles). Le système obtenu est donc de taille inférieure au système total, et de même taille que celui que l’on aurait obtenu avec la méthode du paragraphe précédent. Cependant, le fait de retirer du système global des lignes et des colonnes n’est pas possible à l’heure actuelle dans Code_Aster .

Substitution

Si on ne peut éliminer des degrés de liberté, il faut leur affecter la valeur 0. Pour cela on peut décider de modifier la matrice tangente et le second membre, en mettant une valeur réelle quelconque \({k}_{0}\) (par exemple 1) sur la diagonale de la matrice et 0 dans le second membre, à la position correspondante au degré de liberté à annuler. On a donc un système de même taille, mais numériquement mal-conditionné à cause de la valeur quelconque choisie sur la diagonale. En effet, au niveau du calcul de la matrice élémentaire, on ne connaît pas globalement la matrice de raideur, donc la valeur optimale du paramètre \({k}_{0}\) (en terme de conditionnement) n’est pas connue.

Dualisation

Pour pallier ce genre de problème, le mot-clé DDL_IMPO de l’opérateur AFFE_CHAR_MECA permet d’imposer à un degré de liberté d’un nœud une valeur prédéfinie. Pour résoudre le système linéaire sous contraintes ainsi obtenu, la technique des double multiplicateurs de Lagrange est utilisée [bib 5 ], qui permet un meilleur conditionnement qu’avec la technique simpliste du paragraphe précédent, car le choix des paramètres \({k}_{0}\) est réalisé connaissant la matrice assemblée. Le principal inconvénient est que deux équations supplémentaires sont rajoutées au système pour chaque ddl imposé.

Substitution généralisée

La méthode par substitution est généralisée à l’imposition de ddl à une valeur quelconque (autre que 0) dans Code_Aster (opérateur AFFE_CHAR_CINE). Cependant, cet opérateur ne fonctionne pas à l’heure actuelle avec l’opérateur de statique non linéaire (STAT_NON_LINE) utilisé pour résoudre le système (la non-linéarité du problème traité est due au contact-frottement).

Solution choisie

L’élimination est utile pour les mailles bilinéaires et les fonds de fissure. On le fait en utilisant la méthode par substitution: cependant ce choix doit faire l’objet d’un suivi, voire d’une étude plus approfondi de la robustesse, car il y a de possibles impacts sur la stabilité de la convergence. Pour le reste, à savoir les relations d’égalité, la solution qui a été retenue est celle des doubles multiplicateurs de Lagrange. Notons qu’avec l’utilisation de Mumps comme solveur, un seul multiplicateur est pris en compte. Cependant, on s’est rendu-compte que l’argument sur le mauvais conditionnement qui nous a conduit à ne pas choisir la méthode par substitution ne tient pas. Certes, la matrice peut être mal-conditionnée, mais cela n’a pas de conséquence sur la qualité des résultats car la matrice est diagonale par bloc (par exemple, la matrice diagonale \(\mathit{diag}(1,10.e16)\) est mal conditionnée mais ne pose pas de problème). On pense donc utiliser dans le futur la méthode par substitution (soit en mettant à la main des 1 sur la diagonale, soit en utilisant AFFE_CHAR_CINE quand il sera fonctionnel, ce qui est équivalent) à la place de la dualisation.

Remarque :

L’annulation des degrés de liberté Heaviside et crack-tip enrichis en trop est fait en utilisant la méthode de substitution. En effet pour des problèmes où le maillage est raffiné dans la zone de la fissure, le nombre d’équations rajoutées pour les annuler dans le cas du choix de la dualisation engendrerait des temps de calcul supplémentaires non négligeables.

Calcul de la normale à la facette aux points d’intégration

Tant que les champs des level sets sont interpolés par des fonctions de forme linéaires, on peut admettre une normale \(n\) unique sur la facette de contact, issue du produit vectoriel des côtés de cette facette. Lorsque l’on monte en ordre, la fissure est non- plane et il faut considérer une nouvelle normale en chaque point d’intégration. Celle- ci est issue du gradient de la level set normale, qui résulte de l’approximation aux nœuds de la facette des gradients aux nœuds. Les gradients aux nœuds sont eux-mêmes issus d’une moyenne aux nœuds des gradients élémentaires des éléments connectés au nœud.

Conditionnement pour la méthode pénalisée

Un bon conditionnement de l’équation d’équilibre de la formulation pénalisée impose une «plage conseillée» pour la définition des coefficients de pénalisation, qui est laissée à la discrétion de l’utilisateur. Nous avons:

\({K}_{\text{méca}}\mathrm{~}Eh\text{et}{A}_{u}\mathrm{~}{\kappa}_{n}{h}^{2}\) , ce qui impose \({\kappa}_{n}\) raisonnable devant \(\frac{E}{h}\) .

\({K}_{\text{méca}}\mathrm{~}Eh\text{et}{B}_{u}\mathrm{~}\sigma \mu {\kappa}_{\tau}{h}^{2}\text{}\) , ce qui impose \({\kappa}_{\tau}\) raisonnable devant \(\frac{E}{\mu \sigma h}\) .

Dans les tests, on prend \({\kappa}_{n}\mathrm{~}{10}^{5}\frac{E}{h}\) et \({\kappa}_{\tau}\mathrm{~}{10}^{5}\frac{E}{\mu \sigma h}\) .

Stratégie de résolution

La stratégie de résolution est la même que celle utilisée par la méthode continue dans le cadre éléments finis classique [bib 1 ]. La seule différence est qu’avec X-FEM, nul appariement n’est nécessaire.

Algorithme de résolution

Loi de contact-frottement

Avec X-FEM, les points en vis-à-vis sont connus a priori et cet «appariement» intrinsèque n’évolue pas au cours du calcul (hypothèse de petits déplacements). Ainsi, il n’est pas nécessaire de réaliser une phase d’appariement comme dans le cadre classique. La boucle géométrique est aussi supprimée puisqu’il n’y a pas de ré-actualisation. Pour chaque pas de temps, il reste 3 boucles imbriquées. La boucle sur les seuils de frottement permet de résoudre le problème de frottement par une recherche de point fixe sur les seuils de frottement de Tresca. La boucle sur les statuts de contact (apparentée à la méthode des contraintes actives) permet de déterminer l’espace des points de contact effectifs. Au niveau le plus profond, la boucle sur les itérations de Newton permet de résoudre la non-linéarité restante, celle due à la projection sur la boule unité.

Pour un pas de temps :

Initialisation des seuils de frottement \({\lambda}_{s}\) .

• Boucle sur les seuils de frottementInitialisation des statuts de contact \(\chi\)

  • Boucle sur les statuts de contact

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

    • Fin des itérations de Newton

Actualisation des statuts de contact \(\chi\)

• Fin de la boucle des contraintes actives

Actualisation des seuils de frottement \({\lambda}_{s}\)

• Fin de la boucle sur les seuils de frottement

Loi cohésive.

Dans l’implémentation des lois cohésives, tant mixtes que régularisées, nous n’introduisons pas de champs de statuts comme nous avions pu le faire comme pour le contact: les différents régimes sont gérés par la routine de comportement elle-même, directement dans la méthode de Newton. 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.

Critère d’arrêt de la boucle sur les statuts de contact

Pour les points supposés non-contactant \((\chi =0)\) , on vérifie la condition de non interpénétration \(({d}_{n}\le 0)\) . Le test est le suivant: s’il y a interpénétration \(({d}_{n}>0)\) alors le point est supposé contactant \((\chi =1)\) lors de l’itération de contraintes actives suivante. Numériquement, le test s’écrit \({d}_{n}>{10}^{-16}\) .

Pour les points supposés contactant \((\chi =1)\) , on vérifie que la valeur de la réaction de contact suivant la normale est négative \(({\lambda}_{n}\le 0)\) . Le test est le suivant: s’il y a décollement \(({\lambda}_{n}>0)\) alors le point est supposé non-contactant \((\chi =0)\) lors de l’itération de contraintes actives suivante. Numériquement, le test s’écrit \({\lambda}_{n}>-{10}^{-3}\) .

Remarque :

Le statut de contact étant défini indépendamment pour chaque point d’intégration de chaque facette de contact, les tests d’arrêt sont réalisés en chacun de ces points.

Remarque :

Le jeu en ces points d’intégration est calculé grâce à l’interpolation du champ de déplacement sur l’élément parent (3D) alors que la réaction de contact est calculée grâce à l’interpolation des ddls de contact sur la facette de contact (2D). On pourrait aussi calculer le jeu par interpolation du déplacement aux sommets de la facette de contact (ce dernier étant déterminé avec les fonctions de forme de l’élément parent).

Remarque :

Les valeurs numériques des 2 tests d’arrêts sont délicates à déterminer. Dans certaines configurations, des oscillations de statut de contact apparaissent et empêchent la convergence de l’algorithme. Ce phénomène devrait être identifié, et dans le cas où les valeurs dans les 2 cas (contactant et non-contactant) sont voisines, on pourrait considérer que la convergence de la boucle de contraintes actives est atteinte.

Critère d’arrêt de la boucle sur les seuils de frottement

On considère que la boucle sur les seuils de frottement a convergé si la solution ne change pas trop d’une itération à l’autre. Plus précisément, soit \(v\) le vecteur solution mixte (déplacement, multiplicateurs de contact, semi-multiplicateurs de frottement) , on a convergé à l’itération \(i\) si:

\(\frac{\max\mid {v}^{i}-{v}^{i-1}\mid }{\max\mid {v}^{i-1}\mid }<{10}^{-3}\)

où \(\max(u)\) signifie le max sur toutes les composantes du vecteur \(u\) .

Écriture de la formulation lors d’une itération de Newton

Ré-écrivons la formulation faible à trois champs décrite au paragraphe [§ 3 ] lors d’une itération de Newton. Il faut tenir compte de la boucle sur les seuils de frottement, de celle sur les contraintes actives. Ainsi, à ce niveau, les seuils de frottement deviennent des constantes notées \({\lambda}_{s}\) , les statuts de contact \(\chi ({g}_{n})\) deviennent des constantes \(\chi\) . De plus, le problème restant étant non-linéaire (à cause de la projection sur la boule unité) , il est linéarisé par la méthode de Newton-Raphson.

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.

La projection sur la boule unité s’écrit:

\({P}_{B(0,1)}(x)=\lbrace \begin{array}{}x\text{si}x\in B(0,1)\\ \frac{x}{\parallel x\parallel }\text{sinon}\end{array}\)

La différentielle de cette application, en tout point non situé sur le bord de \(B(0,1)\) , est définie par:

\({\partial}_{x}{P}_{B(0,1)}(x)\mathrm{\delta x}=K(x)\mathrm{\delta x}\)

avec

\(K(x)=\lbrace \begin{array}{ccc}{I}_{d}& \text{si}x\in B(0,1)& (\text{adhérence})\\ \frac{1}{\parallel x\parallel }({I}_{d}-\frac{x\cdot {x}^{T}}{{\parallel x\parallel }^{2}})& \text{sinon}& (\text{glissement})\end{array}\)

Ainsi, la différentielle de \({P}_{B(0,1)}({g}_{\tau})\) sera :

\(K({g}_{\tau}){\mathrm{\delta g}}_{\tau}=K(\Lambda +\frac{{\rho}_{\tau}}{\mathrm{\Delta t}}\Delta [[u]{]}_{\tau})(\mathrm{\delta \Lambda }+\frac{{\rho}_{\tau}}{\mathrm{\Delta t}}[[\mathrm{\delta u}]{]}_{\tau})\)

où \({g}_{\tau}\) est le semi-multiplicateur de frottement issu de l’itération de Newton précédente et \({\mathrm{\delta g}}_{\tau}\) l’incrément d’inconnues. La connaissance du semi-multiplicateur de frottement issu de l’itération de Newton précédente permet de savoir facilement si le point est dans l’état adhérent ou glissant.

De la même façon, dans le cas de la loi cohésive régularisée CZM_EXP_REG, par exemple, 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 =-\mathrm{〚}u\mathrm{〛}\) . 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:

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

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

• \({\mathrm{〚}u\mathrm{〛}}_{\mathit{èq}}\ge \alpha \mathrm{et}\mathrm{〚}{u}_{n}\mathrm{〛}\ge 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(\frac{-{\sigma}_{c}}{{G}_{c}}{〚u〛}_{\mathrm{èq}})\left[\frac{{\sigma}_{c}}{{〚u〛}_{\mathrm{èq}}}({\tau}_{1}\otimes {\tau}_{1}+{\tau}_{2}\otimes {\tau}_{2})-(\frac{{\sigma}_{c}^{2}}{{G}_{c}}+\frac{{\sigma}_{c}}{{〚u〛}_{\mathrm{èq}}})\frac{{〚u〛}_{\tau}\otimes {〚u〛}_{\tau}}{{{〚u〛}_{\mathrm{èq}}}^{2}}\right]\)

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

Linéarisation du problème

Écriture intégrale avec la méthode du Lagrangien augmenté

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 \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) tel que:

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}},{\Lambda}^{\text{*}})\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}}\text{χδλ}[[{u}^{\text{*}}]]\cdot {\mathit{nd\Gamma }}_{c}+{\int}_{{\Gamma}_{c}}{\text{χρ}}_{n}[[\mathit{\delta u}]]\cdot n[[{u}^{\text{*}}]]\cdot {\mathit{nd\Gamma }}_{c}\\ -{\int}_{{\Gamma}_{c}}{\text{χμλ}}_{s}K({g}_{\tau}){\mathit{\delta g}}_{\tau}[[{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}}\chi (\lambda -{\rho}_{n}[[u]]\cdot n)[[{u}^{\text{*}}]]\cdot {\mathit{nd\Gamma }}_{c}\\ +{\int}_{{\Gamma}_{c}}{\text{χμλ}}_{s}{P}_{B(0,1)}({g}_{\tau})\cdot [[{u}^{\text{*}}]{]}_{\tau}{\mathit{d\Gamma }}_{c}\end{array}\)
Loi de contact	\(\begin{array}{c}-{\int}_{\mathit{\Gamma c}}\frac{(1-\chi )}{{\rho}_{n}}\text{δλ}{\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}-{\int}_{\mathit{\Gamma c}}\chi [[\mathit{\delta u}]]\cdot {\mathit{n\lambda }}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ ={\int}_{\mathit{\Gamma c}}\frac{(1-\chi )}{{\rho}_{n}}\lambda {\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}+{\int}_{\mathit{\Gamma c}}\chi [[u]]\cdot n{\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}\end{array}\)
Loi de frottement	\(\begin{array}{c}{\int}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}\mathit{\Delta t}}{{\rho}_{\tau}}\left[\mathit{\delta \Lambda }-K({g}_{\tau}){\mathit{\delta g}}_{\tau}\right]{\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}+{\int}_{\mathit{\Gamma c}}(1-\chi ){\text{δΛΛ}}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ =-{\int}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}\mathit{\Delta t}}{{\rho}_{\tau}}\left[\Lambda -{P}_{B(0,1)}({g}_{\tau})\right]{\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}-{\int}_{\mathit{\Gamma c}}(1-\chi )\text{Λ}{\text{Λ}}^{\text{*}}{\mathit{d\Gamma }}_{c}\end{array}\)

Écriture intégrale avec la méthode pénalisé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 \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) tel que:

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}},{\Lambda}^{\text{*}})\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}}\text{χ}\delta \lambda [[{u}^{\text{*}}]]\cdot {\mathit{nd\Gamma }}_{c}\\ -{\int}_{{\Gamma}_{c}}{\text{χμλ}}_{s}{\kappa}_{\tau}K({\kappa}_{\tau}{\nu}_{\tau})\delta {\nu}_{\tau}[[{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}}\chi \lambda [[{u}^{\text{*}}]]\cdot {\mathit{nd\Gamma }}_{c}\\ +{\int}_{{\Gamma}_{c}}{\text{χμλ}}_{s}{P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau})\cdot [[{u}^{\text{*}}]{]}_{\tau}{\mathit{d\Gamma }}_{c}\end{array}\)
Loi de contact	\(\begin{array}{c}-{\int}_{\mathit{\Gamma c}}\frac{(1-\chi )}{{\kappa}_{n}}\text{δλ}{\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}-{\int}_{\mathit{\Gamma c}}\chi (\frac{\text{δλ}}{{\kappa}_{n}}+[[\mathit{\delta u}]]\cdot n){\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ ={\int}_{\mathit{\Gamma c}}\frac{(1-\chi )}{{\kappa}_{n}}\lambda {\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}+{\int}_{\mathit{\Gamma c}}\chi (\frac{\text{λ}}{{\kappa}_{n}}+[[u]]\cdot n){\lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}\end{array}\)
Loi de frottement	\(\begin{array}{c}{\int}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}}{{\kappa}_{\tau}}\left[\mathit{\delta \Lambda }-\textcolor[rgb]{1,0,0}{K({\kappa}_{\tau}{\nu}_{\tau}){\kappa}_{\tau}\delta {\nu}_{\tau}}\right]{\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}+{\int}_{\mathit{\Gamma c}}(1-\chi ){\text{δΛΛ}}^{\text{*}}{\mathit{d\Gamma }}_{c}\\ =-{\int}_{\mathit{\Gamma c}}\frac{{\text{χμλ}}_{s}}{{\kappa}_{\tau}}\left[\Lambda -{P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau})\right]{\Lambda}^{\text{*}}{\mathit{d\Gamma }}_{c}-{\int}_{\mathit{\Gamma c}}(1-\chi )\text{Λ}{\text{Λ}}^{\text{*}}{\mathit{d\Gamma }}_{c}\end{array}\)

On a choisi ici de résoudre les problèmes de contact et de frottement de façon implicite (les semi-multiplicateurs de contact et de frottement s’expriment en fonction du saut de déplacement de l’itération de Newton courante). Ce choix rend la matrice de rigidité non symétrique, le bloc \({B}_{r}\) (voir [§ 5.6.1 ]) étant non nul à cause du terme en rouge dans l’écriture de la loi de frottement alors que le bloc \({{B}_{r}}^{T}\) est nul.

É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 \((\mathrm{\delta u},\mathrm{\delta \lambda },\mathrm{\delta \Lambda })\in {V}_{0}\times H\times H\) tel que:

\(\forall ({u}^{\text{*}},{\lambda}^{\text{*}},{\Lambda}^{\text{*}})\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}}({t}_{c,n}[[{u}^{\text{*}}]{]}_{n}+{t}_{c,\tau }[[{u}^{\text{*}}]{]}_{\tau}){\mathit{d\Gamma }}_{c}\end{array}\)
Interface: partie normale	\({\int}_{\mathit{\Gamma c}}{\lambda}^{\text{*}}(\lambda +\text{δλ}-{t}_{c,n}\cdot n){\mathit{d\Gamma }}_{c}=0\)
Interface: partie tangentielle	\({\int}_{\mathit{\Gamma c}}{\Lambda}^{\text{*}}(\Lambda +\delta \Lambda -{t}_{c,\tau }){\mathit{d\Gamma }}_{c}=0\)

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

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 (\delta u):\epsilon ({u}^{\text{*}})d\Omega -{\int}_{\Gamma}(\mathit{Id}-r\frac{\partial \delta }{\partial p})\cdot \delta \lambda \cdot \mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \\ +{\int}_{\Gamma}r(\mathit{Id}-r\frac{\partial \delta }{\partial p})\cdot \mathrm{〚}\delta u\mathrm{〛}\cdot \mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \\ =-{\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(\mathrm{〚}u\mathrm{〛}+\delta (p))\right]\cdot \mathrm{〚}{u}^{\text{*}}\mathrm{〛}d\Gamma \end{array}\)

Loi d’interface

\(\begin{array}{c}-{\int}_{{\Gamma}_{c}}(1-r\frac{\partial \delta }{\partial p})\cdot \mathrm{〚}u\mathrm{〛}\cdot {\lambda}^{\text{*}}d\Gamma -{\int}_{\Gamma}\frac{\partial \delta }{\partial p}\cdot \delta \lambda \cdot \lambda \text{*}d\Gamma \\ ={\int}_{\Gamma}(\mathrm{〚}u\mathrm{〛}+\delta )\cdot {\lambda}^{\text{*}}d\Gamma \end{array}\)

Termes élémentaires de contact frottant

Écriture matricielle du problème linéarisé

En reprenant les notations de [bib 1 ], 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}{}\left\lbrace {u}^{\text{*}}\right\rbrace \left[{K}_{\text{méca}}\right](\mathrm{\delta u})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace {\left[A\right]}^{T}(\text{δλ})+\left\lbrace {u}^{\text{*}}\right\rbrace \left[{A}_{u}\right](\mathrm{\delta u})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace {\left[{B}_{r}\right]}^{T}(\mathrm{\delta \Lambda })+\left\lbrace {u}^{\text{*}}\right\rbrace \left[{B}_{u}\right](\mathrm{\delta u})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace \left[{D}_{u}\right](\mathrm{\delta u})\\ =\left\lbrace {u}^{\text{*}}\right\rbrace ({L}_{\text{méca}}^{1})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace ({L}_{\text{cont}}^{1})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace ({L}_{\text{frott}}^{1})\\ +\left\lbrace {u}^{\text{*}}\right\rbrace ({L}_{\text{coh}}^{1})\\ \end{array}\)
Loi de contact	\(\begin{array}{}\left\lbrace {\lambda}^{\text{*}}\right\rbrace \left[C\right](\text{δλ})+\left\lbrace {\lambda}^{\text{*}}\right\rbrace \left[A\right](\mathrm{\delta u})\\ =\left\lbrace {\lambda}^{\text{*}}\right\rbrace ({L}_{\text{cont}}^{2})+\left\lbrace {\lambda}^{\text{*}}\right\rbrace ({L}_{\text{coh}}^{2})\\ \end{array}\)
Loi de frottement	\(\begin{array}{}\left\lbrace {\Lambda}^{\text{*}}\right\rbrace \left[{F}_{r}\right](\mathrm{\delta \Lambda })+\left\lbrace {\Lambda}^{\text{*}}\right\rbrace \left[{B}_{r}\right](\mathrm{\delta u})\\ =\left\lbrace {\Lambda}^{\text{*}}\right\rbrace ({L}_{\text{frott}}^{3})+\left\lbrace {\Lambda}^{\text{*}}\right\rbrace ({L}_{\text{coh}}^{3})\end{array}\)

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}}+{{\rm A}}_{u}+{{\rm B}}_{u}+{D}_{u}\\ {\rm A}\\ {{\rm B}}_{r}\end{array}\begin{array}{c}{{\rm A}}^{T}\\ C\\ 0\end{array}\begin{array}{c}{{\rm B}}_{r}^{T}\\ 0\\ {F}_{r}\end{array}\right](\begin{array}{c}\mathrm{\delta u}\\ \text{δλ}\\ \mathrm{\delta \Lambda }\end{array})=(\begin{array}{c}{L}_{\text{méca}}^{1}+{L}_{\text{cont}}^{1}+{L}_{\text{frott}}^{1}+{L}_{\text{coh}}^{1}\\ {L}_{\text{cont}}^{2}+{L}_{\text{coh}}^{2}\\ {L}_{\text{frott}}^{3}+{L}_{\text{coh}}^{3}\end{array})\)

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.

Bien entendu, les termes de cohésion \({D}_{u}\) , \({L}_{\text{coh}}^{1}\) , \({L}_{\text{coh}}^{2}\) et \({L}_{\text{coh}}^{3}\) n’apparaissent que dans la formulation pour loi cohésive régularisée. Si tel est le cas, tous les autres termes sont nuls, à l’exception des termes \({K}_{\text{méca}}\) et \({L}_{\text{méca}}^{1}\) bien entendu, mais aussi à l’exception de \(C\) , \({F}_{r}\) , \({L}_{\text{cont}}^{2}\) et \({L}_{\text{frott}}^{3}\) qui servent au post-traitement, et qui sont déduits des expressions détaillées à suivre en considérant \(\chi =0\) et \({\rho}_{n}=1\) .

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

\({A}_{u}\) est la matrice de rigidité augmentée due au contact.

\({B}_{u}\) est la matrice de rigidité augmentée due au frottement.

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

\(A\) est la matrice liant les termes de déplacement à ceux de contact (matrice de la loi de contact).

\({B}_{r}\) est la matrice liant les termes de déplacement à ceux de frottement (matrice des lois de frottement). Cette matrice est notée \(B\) dans bib 1 , mais pour ne pas la confondre avec la matrice des dérivées des fonctions de forme, nous la noterons \({B}_{r}\) .

\(C\) est la matrice permettant de déterminer les pressions de contact dans le cas non-contactant.

\({F}_{r}\) est la matrice permettant de déterminer les multiplicateurs de frottement dans le cas de contact non-frottant.

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

\({L}_{\text{cont}}^{1}\) et \({L}_{\text{cont}}^{2}\) sont les seconds membres dus au contact.

\({L}_{\text{frott}}^{1}\) et \({L}_{\text{frott}}^{3}\) sont les seconds membres dus au frottement.

\({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 élémentaires de contact

Méthode du Lagrangien augmenté

Compte-tenu des discrétisations des champs évoquées aux paragraphes [§3.2] de [R7.02.12] et [§ 4.1 ] du présent document, le système matriciel continu ci-dessus est remplacé par un système discret. Plus précisément, la matrice \(A\) a la forme suivante:

\(\begin{array}{ccccccccccc}& & \phantom{[}& \delta {a}_{j}& \delta {b}_{j}& \delta {c}_{j}^{1}& \delta {c}_{j}^{2}& \delta {c}_{j}^{3}& \delta {c}_{j}^{4}& \phantom{]}& \\ {[A]}_{ij}& =& [:ref:`& 0& x& x& 0& 0& 0& <& 0& x& x& 0& 0& 0& >\)]& {lambda}_{i}^{text{*}}end{array}`

En effet, ceci est dû au fait que les termes de contact s’annulent pour les ddls dont la fonction de forme associées est continue. En effet,

\(\begin{array}{c}{\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left[A\right]}_{ij}{\left(\mathit{\delta u}\right)}_{j}=-{\int}_{{\Gamma}^{1}}{\text{χψ}}_{i}{\lambda}_{i}^{\text{*}}{\varphi}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}{\mathit{\delta b}}_{j}+{F}^{1}{\mathit{\delta c}}_{j}^{1}+{F}^{2}{\mathit{\delta c}}_{j}^{2}+{F}^{3}{\mathit{\delta c}}_{j}^{3}+{F}^{4}{\mathit{\delta c}}_{j}^{4}\right)\cdot n\mathit{d\Gamma }\\ +{\int}_{{\Gamma}^{2}}{\text{χψ}}_{i}{\lambda}_{i}^{\text{*}}{\varphi}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}{\mathit{\delta b}}_{j}+{F}^{1}{\mathit{\delta c}}_{j}^{1}+{F}^{2}{\mathit{\delta c}}_{j}^{2}+{F}^{3}{\mathit{\delta c}}_{j}^{3}+{F}^{4}{\mathit{\delta c}}_{j}^{4}\right)\cdot n\mathit{d\Gamma }\\ =-{\int}_{{\Gamma}^{1}}{\text{χψ}}_{i}{\lambda}_{i}^{\text{*}}{\varphi}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}({x}^{-}){\mathit{\delta b}}_{j}-\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right)\cdot n\mathit{d\Gamma }+{\int}_{{\Gamma}^{2}}{\text{χψ}}_{i}{\lambda}_{i}^{\text{*}}{\varphi}_{j}\left({\mathit{\delta a}}_{j}+{H}_{j}({x}^{+}){\mathit{\delta b}}_{j}+\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right).n\mathit{d\Gamma }\\ ={\int}_{\Gamma}{\text{χψ}}_{i}{\lambda}_{i}^{\text{*}}{\varphi}_{j}\left(2{\mathit{\delta b}}_{j}+2\sqrt{r}{\mathit{\delta c}}_{j}^{1}\right)\cdot n\mathit{d\Gamma }\end{array}\)

Ici, on voit clairement apparaître le produit des fonctions de forme \({\psi}_{i}\) du triangle avec les fonctions de forme \({\phi }_{j}\) de l’élément volumique parent.

Remarque:

On note suite à ce calcul l’expression du saut de déplacement en fonction des degrés de libertés enrichis X-FEM : \(\begin{array}{c}{[[u]]}_{j}={\left({a}_{j}+{H}_{j}{b}_{j}+{F}^{1}{c}_{j}^{1}+{F}^{2}{c}_{j}^{2}+{F}^{3}{c}_{j}^{3}+{F}^{4}{c}_{j}^{4}\right)}_{1}\\ \phantom{\rule{14em}{0ex}}-{\left({a}_{j}+{H}_{j}{b}_{j}+{F}^{1}{c}_{j}^{1}+{F}^{2}{c}_{j}^{2}+{F}^{3}{c}_{j}^{3}+{F}^{4}{c}_{j}^{4}\right)}_{2}\\ \phantom{\rule{10em}{0ex}}=\left({a}_{j}+{H}_{j}({x}^{-}){b}_{j}-\sqrt{r}{c}_{j}^{1}\right)-\left({a}_{j}+{H}_{j}({x}^{+}){b}_{j}+\sqrt{r}{c}_{j}^{1}\right)\\ \phantom{\rule{10em}{0ex}}=-\left(2{b}_{j}+2\sqrt{r}{c}_{j}^{1}\right)\end{array}\)

En effet par construction (cf. [R7.02.12]), le coefficient de saut pour les fonctions de sélection de domaine vaut:

\({H}_{j}({x}^{+})-{H}_{j}({x}^{-})=2\)

De même, la matrice de rigidité augmentée due au contact a la forme suivante:

\(\begin{array}{cccc}& & \begin{array}{cccccc}\phantom{[}\delta {a}_{j}& \delta {b}_{j}& \delta {c}_{j}^{1}& \delta {c}_{j}^{2}& \delta {c}_{j}^{3}& \delta {c}_{j}^{4}\phantom{]}\end{array}& \\ {[{A}_{u}]}_{ij}& =& \left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& \times & \times & 0& 0& 0\\ 0& \times & \times & 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& \begin{array}{c}{{a}_{i}}^{\text{*}}\\ {{b}_{i}}^{\text{*}}\\ {{c}_{i}^{1}}^{\text{*}}\\ {{c}_{i}^{2}}^{\text{*}}\\ {{c}_{i}^{3}}^{\text{*}}\\ {{c}_{i}^{4}}^{\text{*}}\end{array}\end{array}\)

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{\left[{A}_{u}\right]}_{ij}{(\text{δu})}_{j}={\int}_{\Gamma}{\text{χρ}}_{n}{\phi }_{i}\left\lbrace ({\mathrm{2b}}_{{i}^{\text{*}}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\cdot n\right\rbrace {\phi }_{j}(2{\text{δb}}_{j}+2\sqrt{r}{\text{δc}}_{j}^{1})\cdot n\mathrm{d\Gamma }\)

L’expression de la matrice \(C\) ne change pas par rapport au cas classique sans X-FEM:

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left[C\right]}_{ij}{(\text{δλ})}_{j}=-{\int}_{\Gamma}\frac{1}{{\rho}_{n}}(1-\chi ){\psi}_{i}{\lambda}_{i}^{\text{*}}{\psi}_{j}{\text{δλ}}_{j}\mathrm{d\Gamma }\)

Méthode pénalisée

La matrice \(C\) a pour expression:

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left[C\right]}_{ij}{(\text{δλ})}_{j}=-{\int}_{\Gamma}\frac{1}{{\kappa}_{n}}{\psi}_{i}{\lambda}_{i}^{\text{*}}{\psi}_{j}{\text{δλ}}_{j}\mathrm{d\Gamma }\)

La matrice \({A}_{u}\) est nulle.

La matrice \(A\) a pour expression:

\(\begin{array}{}{\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{\left[A\right]}_{ij}{(\mathrm{\delta u})}_{j}={\int}_{\Gamma}\chi {\text{ψ}}_{i}{\lambda}_{i}^{\text{*}}{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\cdot n\mathrm{d\Gamma }\end{array}\)

Formulation pour loi cohésive régularisée

Une matrice \(C\) nécessaire au post-traitement est utilisée, dont l’expression est celle du lagrangien augmenté où l’on a écrit \(\chi =0\) et \({\rho}_{n}=1\) . Nous obtenons ainsi:

\({\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 }\)

Expression des seconds membres de contact

Méthode du Lagrangien augmenté

Ces expressions font intervenir des grandeurs à l’itération de Newton précédente (itération \(k-1\) ). C’est pour cela que l’on a fait figurer explicitement la référence à l’indice \(k-1\) :

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{cont}}^{1})}_{i}=-{\int}_{\Gamma}\chi {\phi }_{i}\left\lbrace ({\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\cdot n\right\rbrace ({\lambda}^{k-1}-{\rho}_{n}{d}_{n}^{k-1})\mathrm{d\Gamma }\)

L’expression du vecteur \({L}_{\text{cont}}^{2}\) ne change pas par rapport au cas classique sans X-FEM:

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{cont}}^{2})}_{i}={\int}_{\Gamma}{\psi}_{i}{\lambda}_{i}^{\text{*}}(\frac{1-\chi }{{\rho}_{n}}{\lambda}^{k-1}+{\mathrm{\chi d}}_{n}^{k-1})\mathrm{d\Gamma }\)

Méthode pénalisée

Ces expressions font intervenir des grandeurs à l’itération de Newton précédente (itération \(k-1\) ). C’est pour cela que l’on a fait figurer explicitement la référence à l’indice \(k-1\) :

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{cont}}^{1})}_{i}=-{\int}_{\Gamma}\chi {\phi }_{i}\left\lbrace ({\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}})\cdot n\right\rbrace {\lambda}^{k-1}\mathit{d\Gamma }\)

\({\left\lbrace {\lambda}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{cont}}^{2})}_{i}={\int}_{\Gamma}\frac{1}{{\kappa}_{n}}{\psi}_{i}{\lambda}_{i}^{\text{*}}{\lambda}^{k-1}\mathrm{d\Gamma }+{\int}_{\Gamma}\chi {\psi}_{i}{\lambda}_{i}^{\text{*}}{d}_{n}^{k-1}\mathrm{d\Gamma }\)

Formulation pour loi cohésive régularisée

Un vecteur \({L}_{\text{cont}}^{2}\) nécessaire au post-traitement est utilisé, dont l’expression est celle du lagrangien augmenté où l’on a écrit \(\chi =0\) et \({\rho}_{n}=1\) . Nous obtenons ainsi:

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

Expression des matrices de frottement

Méthode du Lagrangien augmenté

Afin d’exprimer les quantités dans le plan tangent, on utilise l’expression du paragraphe [§ 2.3 ], que l’on écrit sous forme matricielle:

\({u}_{\tau}=(\text{Id}-n\otimes n)u=\left[P\right]u\)

Dans cette expression, la matrice \(P\) désigne l’opérateur de projection sur le plan de normale \(n\) . La matrice de cet opérateur symétrique a pour expression:

\(\left[P\right]=\left[\begin{array}{ccc}1-{n}_{x}^{2}& -{n}_{x}{n}_{y}& -{n}_{x}{n}_{z}\\ -{n}_{x}{n}_{y}& 1-{n}_{y}^{2}& -{n}_{y}{n}_{z}\\ -{n}_{x}{n}_{z}& -{n}_{y}{n}_{z}& 1-{n}_{z}^{2}\end{array}\right]\)

où \({n}_{x}\) , \({n}_{y}\) , \({n}_{z}\) sont les coordonnées de la normale \(n\) telle que définie au [§ 4.5 ]. Avec le choix d’une normale constante par facette, cette matrice, ne dépendant que de la normale, a la même valeur en chaque point de Gauss et peut être calculée une seule fois pour chaque facette.

La matrice de rigidité augmentée due au frottement s’écrit de la manière suivante:

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{\left[{B}_{u}\right]}_{ij}{(\mathrm{\delta u})}_{j}=-{\int}_{\Gamma}{\text{χμλ}}_{s}\frac{{\rho}_{\tau}}{\mathrm{\Delta t}}{\phi }_{i}\left\lbrace {\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\rbrace {\left[P\right]}^{T}\left[{K}_{n}\right]{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\left[P\right]\mathrm{d\Gamma }\)

où la matrice \({K}_{n}\) représente la matrice tangente de la projection sur la boule unité du semi-multiplicateur de frottement augmenté à l’itération de Newton précédente: \({K}_{n}=K({g}_{\tau})\) . C’est une matrice connue.

La matrice \({B}_{r}\) a pour expression:

\({\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{\left[{B}_{r}\right]}_{ij}{(\mathrm{\delta u})}_{j}={\int}_{\Gamma}{\text{χμλ}}_{s}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace \left[{K}_{n}\right]\left[P\right]{\phi }_{j}({\mathrm{2b}}_{j}+2\sqrt{r}{c}_{j}^{1})\mathrm{d\Gamma }\)

La matrice \({F}_{r}\) a pour expression:

\(\begin{array}{}{\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{\left[{F}_{r}\right]}_{ij}{(\mathrm{\delta \Lambda })}_{j}={\int}_{\Gamma}\frac{{\text{χμλ}}_{s}\mathrm{\Delta t}}{{\rho}_{\tau}}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace \left[{I}_{d}-{K}_{n}\right]{\psi}_{j}({\Lambda}_{j}^{1}{\tau}_{j}^{1}+{\Lambda}_{j}^{2}{\tau}_{j}^{2})\mathrm{d\Gamma }\\ +{\int}_{\Gamma}(1-\chi ){\psi}_{i}\left\lbrace \begin{array}{cc}{\Lambda}_{i}^{1\text{*}}& {\Lambda}_{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}(\begin{array}{c}{\Lambda}_{j}^{1}\\ {\Lambda}_{j}^{2}\end{array})\mathrm{d\Gamma }\end{array}\)

Méthode pénalisée

La matrice \({F}_{r}\) a pour expression:

\(\begin{array}{}{\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{\left[{F}_{r}\right]}_{ij}{(\mathrm{\delta \Lambda })}_{j}={\int}_{\Gamma}((1-\chi )+\chi \frac{\mu {\lambda}_{s}}{{\kappa}_{\tau}}){\psi}_{i}\left\lbrace \begin{array}{cc}{\Lambda}_{i}^{1\text{*}}& {\Lambda}_{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}(\begin{array}{c}{\Lambda}_{j}^{1}\\ {\Lambda}_{j}^{2}\end{array})\end{array}\mathrm{d\Gamma }\)

La matrice \({B}_{u}\) a pour expression:

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{\left[{B}_{u}\right]}_{ij}{(\mathrm{\delta u})}_{j}=-{\int}_{\Gamma}\chi {\text{μλ}}_{s}\frac{{\kappa}_{\tau}}{\Delta t}{\phi }_{i}\left\lbrace {\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\rbrace {\left[P\right]}^{T}\left[{K}_{n}\right]{\phi }_{j}(2{\mathrm{\delta b}}_{j}+2\sqrt{r}{\mathrm{\delta c}}_{j}^{1})\left[P\right]\mathrm{d\Gamma }\)

où la matrice \({K}_{n}\) représente la matrice tangente de la projection sur la boule unité du semi-multiplicateur de frottement augmenté à l’itération de Newton précédente: \({K}_{n}=K({\kappa}_{\tau}{v}_{\tau})\) . C’est une matrice connue.

La matrice \({B}_{r}\) a pour expression:

\({\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{\left[{B}_{r}\right]}_{ij}{(\mathrm{\delta u})}_{j}={\int}_{\Gamma}{\text{χμλ}}_{s}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace \left[{K}_{n}\right]\left[P\right]{\phi }_{j}({\mathrm{2b}}_{j}+2\sqrt{r}{c}_{j}^{1})\mathrm{d\Gamma }\)

En méthode pénalisée, la matrice de rigidité n’est pas symétrique. On n’a pas \({B}_{r}^{T}={B}_{r}\) mais une matrice nulle à la place de \({B}_{r}^{T}\) .

Formulation pour loi cohésive régularisée

Une matrice \({F}_{r}\) nécessaire au post-traitement est utilisée, dont l’expression est celle du lagrangien augmenté où l’on a écrit \(\chi =0\) . Nous obtenons ainsi:

\({\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{\left[{F}_{r}\right]}_{ij}{(\mathit{\delta \Lambda })}_{j}={\int}_{\Gamma}{\psi}_{i}\left\lbrace \begin{array}{cc}{\Lambda}_{i}^{1\text{*}}& {\Lambda}_{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}(\begin{array}{c}{\Lambda}_{j}^{1}\\ {\Lambda}_{j}^{2}\end{array})\mathit{d\Gamma }\)

Expression des seconds membres de frottement

Méthode du Lagrangien augmenté

Les seconds membres de frottement ont les expressions suivantes:

\(\begin{array}{}{\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{frott}}^{1})}_{i}=-{\int}_{\Gamma}{\text{χμλ}}_{s}{\phi }_{i}\left\lbrace {\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{1\text{*}}\right\rbrace {\left[P\right]}^{T}{P}_{B(0,1)}({g}_{\tau}^{k-1})\mathrm{d\Gamma }\\ {\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{frott}}^{3})}_{i}=-{\int}_{\Gamma}\frac{{\text{χμλ}}_{s}\mathrm{\Delta t}}{{\rho}_{\tau}}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace ({\Lambda}^{k-1}-{P}_{B(0,1)}({g}_{\tau}^{k-1}))\mathrm{d\Gamma }\\ -{\int}_{\Gamma}(1-\chi ){\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace {\Lambda}^{k-1}\mathrm{d\Gamma }\end{array}\)

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

Méthode pénalisée

Les seconds membres de frottement ont les expressions suivantes:

\(\begin{array}{}{\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{frott}}^{3})}_{i}=-{\int}_{\Gamma}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace ((1-\chi ){\Lambda}^{k-1}+\chi \frac{\mu .{\lambda}_{k}}{{\kappa}_{\tau}}({\Lambda}^{k-1}-{P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau}^{k-1})))\mathrm{d\Gamma }\end{array}\)

\({\left\lbrace {u}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{frott}}^{1})}_{i}=-{\int}_{\Gamma}\chi {\text{μλ}}_{s}{\phi }_{i}\left\lbrace {\mathrm{2b}}_{i}^{\text{*}}+2\sqrt{r}{c}_{i}^{{1}^{\text{*}}}\right\rbrace \left[{P}_{T}\right]{P}_{B(0,1)}({\kappa}_{\tau}{\nu}_{\tau}^{k-1})\mathrm{d\Gamma }\)

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

Formulation pour loi cohésive

Un vecteur \({L}_{\text{frott}}^{3}\) nécessaire au post-traitement est utilisé, dont l’expression est celle du lagrangien augmenté où l’on a écrit \(\chi =0\) . Nous obtenons ainsi:

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

Expression des matrices de cohésion

Soient deux directions de la base fixe \(I\) et \(J\) , de vecteurs unitaires \({e}_{I}\) et \({e}_{J}\) . Introduisons la matrice tangente de la loi cohésive dans la base fixe \([{K}_{\mathit{gl}}]\) de coefficients \({[{K}_{\mathit{gl}}]}_{IJ}={e}_{I}\cdot \frac{\partial \delta }{\partial \mathrm{〚}u\mathrm{〛}}\cdot {e}_{J}\) . Avec l’expression de \(\frac{\partial {t}_{c}}{\partial \mathrm{〚}u\mathrm{〛}}\) donnée au [§ 5.4 ], 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}[{K}_{\mathit{loc}}][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 \([{D}_{u}]\) est alors donnée par:

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

Expression des seconds membres de cohésion

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

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

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

\(\begin{array}{c}{\left\lbrace {\Lambda}^{\text{*}}\right\rbrace }_{i}{({L}_{\text{coh}}^{3})}_{i}=-{\int}_{\Gamma}{\psi}_{i}\left\lbrace {\Lambda}_{i}^{1\text{*}}{\tau}_{i}^{1}+{\Lambda}_{i}^{2\text{*}}{\tau}_{i}^{2}\right\rbrace {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.

Pour exprimer \({t}_{c}\) dans la base globale, on peut utiliser \([Q]\) .

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 \(I\) et \(J\) , de vecteurs unitaires \({e}_{I}\) et \({e}_{J}\) . 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{gl}}]}_{\mathit{IJ}}={e}_{I}\cdot \frac{\partial \delta }{\partial p}\cdot {e}_{J}\) . 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}[{K}_{\mathit{loc}}][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{\phi }_{i}{b}_{i}^{\text{*}}r([\mathit{Id}]-r[{K}_{\mathit{gl}}])2{\phi }_{j}{b}_{j}d\Gamma\)

Par ailleurs, nous choisissons de discrétiser \(\lambda\) en base locale. Le coefficient \((1,J)\) de la matrice \(\lbrace {\lambda}_{n}^{\text{*}}\rbrace {[{A}_{u}]}_{ij}{\lbrace \delta u\rbrace }_{j}\) est alors \({\int}_{\Gamma}2{\phi }_{i}{\lambda}_{i}^{\text{*}}n\cdot (\mathit{Id}-r\frac{\partial \delta }{\partial p})\cdot {e}_{j}2{\phi }_{j}{b}_{j}d\Gamma\) . Pour les coefficients \((2,J)\) et \((3,J)\) , 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}_{u}]}_{ij}{\lbrace \delta u\rbrace }_{j}={\int}_{\Gamma}{\varphi}_{i}{\lambda}_{i}^{\text{*}}\left([\mathit{Id}]-r[{K}_{\mathit{loc}}]\right)\cdot [Q]2{\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}{\phi }_{i}{\lambda}_{i}^{\text{*}}[{K}_{\mathit{loc}}]{\phi }_{j}{\lambda}_{j}d\Gamma\)

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

Le coefficient \(I\) du vecteur \({\lbrace u\rbrace }_{i}^{\text{*}}{({L}_{\mathit{coh}}^{1})}_{i}\) a pour expression \({\int}_{\Gamma}2{b}_{i}^{\text{*}}{\phi }_{i}(-\lambda \cdot {e}_{I}+r(\mathrm{〚}u\mathrm{〛}+\delta )\cdot {e}_{I})d\Gamma\) . Avec les notations indroduites dans la partie précédentes, nous déduisons:

\({\lbrace {u}_{i}\rbrace }^{\text{*}}{({L}_{\mathit{coh}}^{1})}_{i}={\int}_{\Gamma}2{b}_{i}^{\text{*}}{\phi }_{I}(-{[Q]}^{T}\cdot \lbrace \lambda \rbrace +r(\lbrace \mathrm{〚}u\mathrm{〛}\rbrace +{[Q]}^{T}\cdot \lbrace \delta \rbrace ))d\Gamma\)

où \(\lbrace \delta \rbrace ` et :math:\)lbrace lambda rbrace ` sont donnés en base locale, et :math:`lbrace mathrm{〚}umathrm{〛}rbrace ` en base fixe.

Le coefficient \(1\) du vecteur \({\lbrace \lambda \rbrace }_{i}^{\text{*}}{({L}_{\mathit{coh}}^{2})}_{i}\) s’écrit \({\int}_{\Gamma}{\lambda}_{i}^{\text{*}}(\mathrm{〚}u\mathrm{〛}\cdot n+\delta \cdot n)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{*}}{({L}_{\mathit{coh}}^{2})}_{i}={\int}_{\Gamma}{\lambda}_{i}^{\text{*}}([Q]\lbrace \mathrm{〚}u\mathrm{〛}\rbrace +\lbrace \delta \rbrace )d\Gamma\)

où \(\lbrace \delta \rbrace ` est donné en base locale, et :math:\)lbrace mathrm{〚}umathrm{〛}rbrace ` en base fixe.

LBB condition

Les approximations choisies d’une part pour le déplacement et d’autre part pour les pressions de contact ne semblent par satisfaire la condition inf-sup dans tous les cas. Le non-respect de la LBB condition engendre des oscillations des pressions de contact, phénomène comparable à celui rencontré en incompressibilité [bib 13 ]. Physiquement, dans le cas du contact lagrangien, cela revient à vouloir imposer le contact en trop de points de l’interface (sur-contrainte), rendant le système hyperstatique. Pour le relâcher, il faut restreindre l’espace des multiplicateurs de Lagrange, comme cela est fait dans [bib 14 ] pour les conditions de Dirichlet avec X-FEM. L’algorithme proposé par Moës [bib 14 ] pour réduire les oscillations est étendu au cas 3D. Son but est d’imposer des relations d’égalité entre multiplicateurs de Lagrange. Cet algorithme a fait l’objet d’une amélioration pour le rendre plus physique et plus efficace.

Dans le cas du contact pénalisé on retrouve les mêmes oscillations. Pour un maillage triangle par exemple, on peut montrer qu’il n’existe pas de combinaison des degrés de liberté de Heaviside qui permettent une rotation de la surface moyenne de la fissure sans générer des oscillations du saut de déplacement. Pour une raideur d’interface élevée, comme c’est le cas en pénalisation, ceci génère des oscillations de pression. Pour y remédier, il est nécessaire de récupérer la pression explicite dans le degré de liberté \(\lambda\) , de lui appliquer la condition de LBB et de la faire remonter dans l’équilibre, ce qui explique la formulation donnée au § 3.3 . L’actualisation des statuts de contact se fait comme dans le cas lagrangien, où \(\lambda\) et non \({d}_{n}\) est testé pour passer d’un état contactant à un état non contactant, ceci afin d’éviter les oscillations sur les statuts de contact.

Description de l’algorithme de Moës pour les éléments linéaires et quadratiques (algo1)

L’algorithme introduit par Moës [bib 14 ] y est présenté dans le but d’imposer des conditions de Dirichlet sur une interface dans le cadre d’X-FEM. Il montre que la technique des multiplicateurs de Lagrange pour imposer des conditions de Dirichlet doit être utilisée avec soin, car la condition inf-sup n’est pas toujours respectée. Le papier se restreint au cas 2D, mais l’algorithme présenté est facilement généralisable au cas 3D. La première phase est une phase de sélection des nœuds, dans laquelle les nœuds sélectionnés sont ceux qui sont «importants» pour l’approximation des multiplicateurs de Lagrange. Les autres nœuds sont surabondants et amènent à des oscillations des multiplicateurs de Lagrange. Une fois les nœuds «importants» sélectionnés, des relations d’égalités sont imposées entre les multiplicateurs de Lagrange, afin de restreindre l’espace des multiplicateurs. Ainsi, les multiplicateurs de Lagrange des arêtes émanant d’un même nœud sélectionné sont égaux.

De manière plus formelle, soient \(E\) et \(N\) les ensembles contenant toutes les arêtes et tous les nœuds du maillage. Les deux extrémités d’une arête \(e\in E\) sont notées \(({v}_{1}(e),{v}_{2}(e))\in {N}^{2}\) . Tout d’abord, on commence par une phase d’initialisation (itération \(k=0\) de l’algorithme). On détermine tout d’abord \({S}_{e}^{0}\) , l’ensemble des arêtes qui sont coupées par l’interface. L’interface étant représentée par la level set normale \(\text{lsn}\) , une arête \(e\in E\) est coupée par l’interface si et seulement si \(\text{lsn}({v}_{1}(e))\cdot \text{lsn}({v}_{2}(e))\le 0\) . Notons que si l’interface coïncide avec le nœud \({v}_{1}(e)\) ou le nœud \({v}_{2}(e)\) , l’arête \(e\) est appartient à \({S}_{e}^{0}\) :

\({S}_{e}^{0}=\left\lbrace e\in E,\text{lsn}({v}_{1}(e))\cdot \text{lsn}({v}_{2}(e))\le 0\right\rbrace\)

Soit \({N}_{e}\) l’ensemble des nœuds connectés par les éléments de \({S}_{e}^{0}\) :

\({N}_{e}=\left\lbrace n\in N,\exists e\in {S}_{e}^{0}n={v}_{1}(e)\text{ou}n={v}_{2}(e)\right\rbrace\)

Soit \({S}_{n}^{0}\) l’ensemble des nœuds sélectionnés à l’itération \(k=0\) (initialisation). Ces nœuds sont ceux qui coïncident avec l’interface (cet ensemble peut être vide):

\({S}_{n}^{0}=\left\lbrace n\in {N}_{e},\text{lsn}(n)=0\right\rbrace\)

Après cette phase d’initialisation, l’algorithme itère de \(k=1,\text{nmax\_iter}\) .

À chaque itération, on effectue les étapes suivantes:

• Mise à jour de l’ensemble de l’ensemble des arêtes: on supprime celles qui sont connectées à un nœud sélectionné à l’itération précédente

\({S}_{e}^{k}={S}_{e}^{k-1}\setminus \left\lbrace e\in {S}_{e}^{k-1},{v}_{1}(e)\in {S}_{n}^{k-1}\text{ou}{v}_{2}(e)\in {S}_{n}^{k-1}\right\rbrace\)

• Calcul du score des nœuds: pour chaque nœud dans \({N}_{e}\) , on calcule un score composé de 2 chiffres: le premier chiffre correspond au nombres d’arêtes dans \({S}_{e}^{k}\) connectées, et le deuxième correspond à la valeur absolue de la level set normale en ce nœud. Ce score \(\text{sc\_no}\) est donc une matrice à deux colonnes dont les lignes représentent le nœud.

\(\forall n\in {N}_{e}\lbrace \begin{array}{}{\text{sc\_no}}^{k}(n,1)=\text{nombre}\text{d'arêtes}\text{connectées}\text{au}\text{noeud}n\\ {\text{sc\_no}}^{k}(n,2)=\mid \text{lsn}(n)\mid \end{array}\)

• Calcul du score des arêtes: pour chaque arête dans \({S}_{e}^{k}\) , on calcule un score composé de 2 chiffres: le premier chiffre correspond à la valeur absolue de la différence du 1er chiffre du score des 2 nœuds extrémités, et le deuxième correspond à un rapport entre les valeurs du 2ème chiffre des 2 nœuds extrémités (c’est-à-dire un rapport de valeurs de \(\text{lsn}\) . Ce score \(\text{sc\_ar}\) est donc une matrice à deux colonnes dont les lignes représentent l’arête.

\(\begin{array}{}\forall e\in {S}_{e}^{k},\forall j\in \left\lbrace 1,2\right\rbrace ,{s}_{j}={\text{sc\_no}}^{k}({v}_{j}(e),1),{l}_{j}={\text{sc\_no}}^{k}({v}_{j}(e),2)\\ \lbrace \begin{array}{}{\text{sc\_ar}}^{k}(e,1)=\mid {s}_{1}-{s}_{2}\mid \\ {\text{sc\_ar}}^{k}(e,2)=\lbrace \begin{array}{cc}\frac{{l}_{1}}{{l}_{1}+{l}_{2}}& \text{si}{s}_{1}<{s}_{2}\\ \frac{{l}_{1}}{{l}_{1}+{l}_{2}}& \text{si}{s}_{1}>{s}_{2}\\ \frac{\min({l}_{1,}{l}_{2})}{{l}_{1}+{l}_{2}}& \text{si}{s}_{1}={s}_{2}\end{array}\end{array}\end{array}\)

• Recherche de la «meilleure arête» \({b}_{e}\) : c’est l’arête dont le 1er chiffre de son score est le plus grand. En cas d’égalité entre 2 arêtes, c’est celle dont le 2ème chiffre de son score est le plus grand.

• Recherche du «meilleur nœud» \({b}_{n}\) : c’est le nœud extrémité de \({b}_{e}\) dont le 1er chiffre de son score est le plus grand. En cas d’égalité, c’est le nœud dont le 2ème chiffre du score est le plus petit (le nœud le plus proche de l’interface). Le nœud \({b}_{n}\) est le seul nœud sélectionné à cette itération:

\({S}_{n}^{k}=\left\lbrace {b}_{n}\right\rbrace\)

L’algorithme s’arrête si lors d’une itération, l’ensemble \({S}_{e}^{k}\) devient l’ensemble vide. L’ensemble final des nœuds sélectionnés sera alors:

\(W=\underset{k}{\cup}{S}_{n}^{k}\)

Après cette phase de sélection des nœuds, l’algorithme construit l’espace des multiplicateurs de Lagrange, dont la taille est égale à celle de \(W\) . Ainsi, l’espace des multiplicateurs est:

\({S}_{\lambda}=\left\lbrace {\lambda}^{i},i\in \left\lbrace 1,\text{card}(W)\right\rbrace \right\rbrace\)

Avec cet algorithme tout nœud de \({N}_{e}\text{/}W\) est connecté par une arête de \({S}_{e}^{0}\) à un nœud de \(W\) . On impose alors une relation d’égalité entre ces deux nœuds. En cas de conflit (connections à plusieurs nœuds de \(W\) par autant d’arêtes), on discrimine les nœuds de \(W\) par le critère 2 des nœuds (plus faible level-set normale).

A noter que cet algorithme peut aussi être utilisé en grands glissements.

Description de l’algorithme modifié pour les éléments linéaires et quadratiques (algo2)

Basé sur des idées similaires, un nouveau algorithme a été proposé. Ainsi, dans la nouvelle version, on part de l’ensemble des arêtes sur les lesquelles la level set normale s’annule au moins en un point. Ces arêtes relient des points de part et d’autre de l’interface (ou éventuellement des points sur l’interface). L’algorithme cherche le sous-ensemble minimal d’arêtes permettant de relier tous les points extrémités des arêtes. Ensuite, des groupes d’arêtes connectées en sont extraits. Les relations imposées sont alors les suivantes:

• les multiplicateurs sur les nœuds sommets de chaque groupe sont imposés égaux,

De manière plus formelle, soient \(E\) et \(N\) les ensembles contenant toutes les arêtes et tous les nœuds du maillage. Les deux extrémités d’une arête \(e\in E\) sont notées \(({v}_{1}(e),{v}_{2}(e))\in {N}^{2}\) . On détermine tout d’abord \({S}_{e}\) , l’ensemble des arêtes qui sont strictement coupées par l’interface. L’interface étant représentée par la level set normale \(\text{lsn}\) , une arête \(e\in E\) est strictement coupée par l’interface si et seulement si \(\text{lsn}({v}_{1}(e))\cdot \text{lsn}({v}_{2}(e))<0\) . Notons que si l’interface coïncide avec le nœud \({v}_{1}(e)\) ou le nœud \({v}_{2}(e)\) , l’arête \(e\) n’est pas dans \({S}_{e}\) :

\({S}_{e}=\left\lbrace e\in E,\text{lsn}({v}_{1}(e))\cdot \text{lsn}({v}_{2}(e))<0\right\rbrace\)

Soit \({N}_{e}\) l’ensemble des nœuds connectés par les éléments de \({S}_{e}\) . On sépare \({N}_{e}\) en deux parties: les nœuds «dessous» et «dessus» la fissure, selon le signe de \(\text{lsn}\) :

\(\begin{array}{}{N}_{e}=\left\lbrace n\in N,\exists e\in {S}_{e}n={v}_{1}(e)\text{ou}n={v}_{2}(e)\right\rbrace \\ {N}_{e}^{+}=\left\lbrace n\in {N}_{e},\text{lsn}(n)>0\right\rbrace \text{et}{N}_{e}^{-}=\left\lbrace n\in {N}_{e},\text{lsn}(n)<0\right\rbrace \end{array}\)

L’algorithme recherche \({S}_{\text{ve}}\) , le sous-ensemble minimal de \({S}_{e}\) qui permet de relier les nœuds dans \({N}_{e}^{+}\) aux nœuds dans \({N}_{e}^{-}\) . Chaque nœud dans \({N}_{e}^{+}\) doit être relié à au moins un nœud dan \({N}_{e}^{-}\) s , et chaque nœud dans \({N}_{e}^{-}\) doit être relié à au moins un nœud dans \({N}_{e}^{+}\) . Les arêtes dans \({S}_{\text{ve}}\) sont appelées «arêtes vitales», car si une de ces arêtes disparaît, au moins un nœud dans \({N}_{e}\) sera orphelin. Cet ensemble d’arêtes vitales n’est pas nécessairement unique. En présence de choix, l’arête vitale la plus courte est privilégiée. Comme on le verra par la suite, cela revient à minimiser la région d’approximation P0. Pour la recherche de l’ensemble \({S}_{\text{ve}}\) , nous avons choisi un algorithme basé sur la notions de scores de nœuds et d’arête, notion que l’on trouve dans l’algo1. L’algorithme va enlever une à une toutes les arêtes non-vitales, jusqu’à ce qu’il ne reste plus que des arêtes vitales. Plus précisément, on associe un score à chaque nœud, qui correspond au nombre d’arêtes connectées à ce nœud. À chaque arête, on associe un score, qui correspond au minimum des scores des deux nœuds extrémités. Soit \(e\) l’arête ayant le score le plus élevé (à score identique, l’arête la plus longue est privilégiée). Si le score de \(e\) est égal à 1, alors toutes les arêtes qui restent sont des arêtes vitales. \({S}_{\text{ve}}\) est déterminé et l’algorithme s’arrête. Si le score de \(e\) est strictement supérieur à 1, l’arête \(e\) est une arête non-vitale, et elle symboliquement enlevée de la liste des arêtes \({S}_{e}\) . L’algorithme recommence, avec un nouveau calcul du score des nœuds, et ainsi de suite jusqu’à ce qu’il ne reste plus que des arêtes vitales.

Il est important de noter que \({S}_{\text{ve}}\) est composé de certaines arêtes déconnectées, et de certaines arêtes connectées entre elles. Ces groupes d’arêtes vitales connectées sont extraits de \({S}_{\text{ve}}\) . Notons que dans un tel groupe, toutes les arêtes sont connectées par un unique nœud (voir Figure 6.2-1 ). Soit \({G}_{\text{cve}}^{i}\) le groupe d’arêtes vitales connectées par le nœud \(i\) . Alors \({G}_{\text{cve}}^{i}\) est défini par:

\({G}_{\text{cve}}^{i}=\left\lbrace e\in {S}_{\text{ve}},i={v}_{1}(e)\text{ou}i={v}_{2}(e)\right\rbrace\)

Maintenant, on impose des relations entre les multiplicateurs de Lagrange. Tous les multiplicateurs portés par les nœuds sommets d’un même groupe sont imposés égaux.

Nous illustrons ces algorithmes sur le cas 2D de la Figure 6.2-1 . Les groupes de nœuds liés par des relations d’égalité avec la version 1 sont marqués par les cercles bleus. Les groupes de nœuds liés par des relations d’égalité avec la version 2 sont marquées par les arêtes pleines reliées entre elles.

A noter que cet algorithme peut aussi être utilisé en grands glissements.

Figure 6.2-1 : Exemple d’arêtes coupées par une interface et approximation résultante

Description de l’algorithme pour les éléments quadratiques (algo 3) :

Cet algorithme permet d’imposer le nombre de relations d’égalité minimal afin d’atteindre le taux de convergence optimal pour les éléments quadratiques. La détermination des arêtes portant les relations d’égalité entre les «Lagrange» se déroule en deux temps:

• Étude des composantes connexes:

Une composante connexe est un ensemble d’arêtes coupées liées entre elles de façon continue (par l’un des deux sommets de chaque arête). L’algorithme analyse l’ensemble des arêtes coupées afin de déterminer les composantes connexes, ensuite il met une relation d’égalité sur la première arête de chaque composante.

Figure 6.3-1: Relations d’égalité sur composantes connexes

• Étude du critère de proximité:

Il s’agit d’analyser, pour chaque arête coupée, la proximité des nœuds par rapport à l’interface. L’idée est de lier les «Lagrange» des nœuds à «faible» influence. Concrètement on commence par éliminer les arêtes non vitales afin d’aboutir à un sous-ensemble d’arêtes coupées liant tous les nœuds de part et d’autre de l’interface. Ensuite les relations d’égalité sont mises sur chaque arête où le critère de proximité est vérifié. Finalement on réétudie les arêtes non vitales: si l’une parmi elles connecte deux arêtes dont une seule vérifie le critère de proximité, l’arête non vitale en question porte une égalité; dans les autres cas on ne fait rien.

Figure 6.3-2: Selection du groupe d’arêtes vitales

Figure 6.3-3: Relations d’égalité sur composantes connexes et critère de proximité

De manière plus formelle, et en reprenant les notations utilisées au §6.2, on considère \({S}_{\text{e}}\) l’ensemble des arêtes strictement coupées par l’interface.

La première étape est de partitionner Se en une union de composantes connexes disjointes: \({S}_{\text{e}}={U}_{i}{S}_{\text{ei}}\)

Une relation d’égalité est imposée sur la première arête de chaque sous-ensemble \({S}_{\text{ei}}\) .

Ensuite on construit \({S}_{\text{ve}}\) , le groupe des arêtes «vitales» (au sens de l’algo 2), pour étudier le critère de proximité. \({S}_{\text{ve}}\) est composé de deux types d’arêtes: des arêtes indépendantes, et des arêtes connectées entre elles via un sommet commun. On note \({S}_{\text{ind}}\) et \({S}_{\text{cte}}\) les sous-ensembles de \({S}_{\text{ve}}\) correspondant respectivement à ces deux types, tels que: \({S}_{\text{e}}={S}_{\text{ind}}U{S}_{\text{cte}}\)

Le critère de proximité est étudié sur les nœuds suivants:

• pour les arêtes de \({S}_{\text{ind}}\) on cherche le nœud le plus proche de l’interface et si la distance (relative) est inférieure à la distance seuil fixée par le critère, l’arête porte une relation d’égalité;

• pour les arêtes de \({S}_{\text{cte}}\) on calcule la distance minimale par rapport au sommet commun, et si cette distance minimale est inférieure au seuilon impose des relations d’égalité sur toutes les arêtes du paquet.

Ainsi on parvient à lier les «Lagrange» à «faible» influence.

A noter que cet algorithme ne peut pas être utilisé dans le cadre de grands glissements, les différents appariements possibles changeant au fur et à mesure du glissement.

Relations imposées entre les semi-multiplicateurs de frottement

Lorsque le frottement est pris en compte, on observe le même phénomène d’oscillations sur les semi-multiplicateurs de frottement que dans le cas précédent sur les multiplicateurs de contact. Pour que la réaction de contact n’oscille plus, il est nécessaire de supprimer les oscillations de la réaction normale (pression de contact) et de réaction tangentielle (donc des semi-multiplicateurs de frottement). Pour cela, il faut donc aussi activer l’algorithme de restriction des espaces des multiplicateurs pour les semi-multiplicateurs de frottement.

Les relations sont à imposer sur les réactions tangentielles, et font intervenir les inconnues de frottement \({\Lambda}_{1}\) et \({\Lambda}_{2}\) ainsi que les vecteurs de la base covariante \({\tau}_{1}\) et \({\tau}_{2}\) .

Pour le cas d’imposition d’une relation d’égalité entre les nœuds \(A\) et \(B\) , la relation s’écrit:

\({\Lambda}_{1}^{A}{\tau}_{1}^{A}+{\Lambda}_{2}^{A}{\tau}_{2}^{A}={\Lambda}_{1}^{B}{\tau}_{1}^{B}+{\Lambda}_{2}^{B}{\tau}_{2}^{B}\)

Les deux inconnues à déterminer (\({\Lambda}_{1}^{A}\) et \({\Lambda}_{2}^{A}\) par exemple) étant scalaires, il faut transformer la relation vectorielle précédente en deux relations scalaires. Ceci est fait par projection sur la base \(({\tau}_{1}^{A},{\tau}_{2}^{A})\) . Les deux relations à imposer sont donc finalement:

\(\lbrace \begin{array}{}{\Lambda}_{1}^{A}({\tau}_{1}^{A}\cdot {\tau}_{1}^{A})+{\Lambda}_{2}^{A}({\tau}_{2}^{A}\cdot {\tau}_{1}^{A})={\Lambda}_{1}^{B}({\tau}_{1}^{B}\cdot {\tau}_{1}^{A})+{\Lambda}_{2}^{B}({\tau}_{2}^{B}\cdot {\tau}_{1}^{A})\\ {\Lambda}_{1}^{A}({\tau}_{1}^{A}\cdot {\tau}_{2}^{A})+{\Lambda}_{2}^{A}({\tau}_{2}^{A}\cdot {\tau}_{2}^{A})={\Lambda}_{1}^{B}({\tau}_{1}^{B}\cdot {\tau}_{2}^{A})+{\Lambda}_{2}^{B}({\tau}_{2}^{B}\cdot {\tau}_{2}^{A})\end{array}\) éq 6.4-1

Ce choix est remis en cause par l’introduction des grands glissements [R5.03.53]: en effet la base de contact change à chaque itération géométrique. Les relations précédentes introduites en dur en utilisant la base de départ n’ont donc plus de sens. Pour résoudre le conflit, on introduit plutôt les deux relations d’égalités sur les composantes :

\(\lbrace \begin{array}{}{\Lambda}_{1}^{A}={\Lambda}_{1}^{B}\\ {\Lambda}_{2}^{A}={\Lambda}_{2}^{B}\end{array}\) éq 6.4-2

En considérant que les bases \(({\tau}_{1}^{A},{\tau}_{2}^{A})\) et \(({\tau}_{1}^{B},{\tau}_{2}^{B})\) sont quasiment identiques du fait de la proximité des points \(A\) et \(B\) , les équations et sont quasiment équivalentes.

Remarques sur les relations imposés par l’algorithme 1 ou 2

Sur les multiplicateurs de contact

Soit la relation linéaire entre les multiplicateurs de Lagrange de contact \({\lambda}_{1}\) , \({\lambda}_{2}\) et \({\lambda}_{3}\) :

\({\lambda}_{2}={\text{αλ}}_{1}+(1-\alpha ){\lambda}_{3}\)

La relation porte sur la valeur de la pression et non pas le vecteur pression de contact. En cas de structure courbe une relation sur les vecteurs pression du type:

:math:`{\lambda}_{2}{n}_{2}={\text{αλ}}_{1}{n}_{1}+(1-\alpha ){\lambda}_{3}{n}_{3}`

n’est pas forcément possible car le vecteur :math:`{n}_{2}`n’est pas une inconnue.

Sur les semi-multiplicateurs de frottement

L’inconnue pour le frottement est vectorielle. On pourrait imaginer lier les vecteurs réactions tangentes entre eux:

\({r}_{{\tau}_{2}}={\mathrm{\alpha r}}_{\text{τ1}}+(1-\alpha ){r}_{\text{τ3}}\)

or

:math:`{r}_{{\tau}_{i}}={\text{μλ}}_{i}{\Lambda}_{i}`

ce qui donne

:math:`{\Lambda}_{2}=\frac{\alpha {\lambda}_{1}{\Lambda}_{1}+(1-\alpha ){\lambda}_{3}{\Lambda}_{3}}{\alpha {\lambda}_{1}+(1-\alpha ){\lambda}_{3}}=\beta {\Lambda}_{1}+(1-\beta ){\Lambda}_{3}`

avec    :math:`\beta =\frac{\alpha {\lambda}_{1}}{\alpha {\lambda}_{1}+(1-\alpha ){\lambda}_{3}}\approx \alpha`        si le maillage est assez fin

Figure 6.5.2-1

Ce choix de relation est impossible car en 3D, si les points 1 et 3 sont en contact glissant, alors le point 2 ne le sera pas! En effet, la norme de :math:`{Lambda}_{2}`sera strictement inférieure à 1 si les directions de glissement ne sont pas colinéaires (voir ).

Figure 6.5.2-1 : Cas d’un point adhérant entre deux points glissants

La solution proposée est d’imposer une relation linéaire entre les normes des vecteurs réactions tangentes, ce qui est équivalent à imposer une relation entre les normes des semi-multiplicateurs de frottement:

\(\parallel {r}_{\text{τ2}}\parallel =\alpha \parallel {r}_{\mathrm{\tau 1}}\parallel +(1-\alpha )\parallel {r}_{\mathrm{\tau 3}}\parallel \iff \parallel {\Lambda}_{2}\parallel =\beta \parallel {\Lambda}_{1}\parallel +(1-\beta )\parallel {\Lambda}_{3}\parallel\) éq. 6.5.2-1

Cette relation est non-linéaire du fait de la norme. La méthode de Newton permet de se ramener à l’imposition successive des relations linéaires. À l’itération de Newton i, la relation est:

:math:`{\mathrm{\delta \Lambda }}_{2}^{i}.\frac{{\Lambda}_{2}^{i-1}}{\parallel {\Lambda}_{2}^{i-1}\parallel }-\text{βδ}{\Lambda}_{1}^{i}.\frac{{\Lambda}_{1}^{i-1}}{\parallel {\Lambda}_{1}^{i-1}\parallel }-(1-\beta ){\mathrm{\delta \Lambda }}_{3}^{i}.\frac{{\Lambda}_{3}^{i-1}}{\parallel {\Lambda}_{3}^{i-1}\parallel }=\beta \parallel {\Lambda}_{1}^{i-1}\parallel +(1-\beta )\parallel {\Lambda}_{3}^{i-1}\parallel -\parallel {\Lambda}_{2}^{i-1}\parallel`


Ce type de relation linéaire n’est actuellement pas disponible dans Code_Aster, où seules les relations linéaires dont les coefficients sont constants pendant toute la durée du calcul sont autorisées.


Actuellement, la relation entre les semi-multiplicateurs de frottement implantée est la suivante:

:math:`{\Lambda}_{2}={\mathrm{\alpha \Lambda }}_{1}+(1-\alpha ){\Lambda}_{3}`


Lorsque le glissement ou l’adhérence sont unidirectionnels, on retrouve bien l’équation [:external:ref:`éq. 4.6-1 <éq. 4.6-1>`] en ayant substitué :math:`\alpha`à :math:`\beta`.

Eléments non simpliciaux (quadrangles, hexaèdres…)

Dans le cas de mailles quaden 2D, pentaou hexaen 3D, on note certaines configurations où des nœuds de la maille n'appartiennent pas à une arête coupée. Sur les deux exemples de la figure , les nœuds rouges ne sont reliés à aucun autre nœud. Afin de satisfaire la LBB, l'idée est de relier ces nœuds aux nœuds connexes mais en essayant de ne pas introduire dans ce cas des relations linéaires.

Figure 6.5.3-1: Les nœuds rouges n’appartiennent pas à une arête coupée.

Pour éviter d'ajouter des relations linéaires, l'idée est d'éliminer de l'interpolation les nœuds qui n'appartiennent pas à une arête coupée. Pour satisfaire la partition de l'unité des contributions de contact, on répartit les fonctions de formes des nœuds éliminés sur les autres nœuds en utilisant la répartition uniforme suivante :


:math:`\tilde{\phi }{\rbrace _{i\in {N}_{\mathrm{actif}}}={\phi }_{i}+{\sum}_{j\in {N}_{\mathrm{elim}}}{\phi }_{j}/{N}_{\mathrm{actif}}`                            éq. 6.5.3-1


Où :math:`{N}_{\mathrm{actif}}`est l'ensemble des nœuds de l'élément directement coupés ou appartenant à une arête coupée. :math:`{N}_{\mathrm{elim}}`est l'ensemble des nœuds à éliminer. Pour le quadrangle de la figure , les fonctions de forme modifiées s'écrivent alors :math:`\tilde{\phi }{\rbrace _{i=1,2,4}={\phi }_{i}+{\phi }_{3}/3`. Pour l'hexaèdre elles s'écrivent:math:`\tilde{\phi }{\rbrace _{i=1,2,4,5,6,7,8}={\phi }_{i}+{\phi }_{3}/7`.

On remarque que d'autres choix de répartition sont possibles: on aurait par exemple pu choisir de répartir un nœud éliminé sur ses nœuds actifs connexes. Pour le quadrangle de la figure , les fonctions de forme modifié s'écriraient alors :math:`\tilde{{\phi }_{1}}={\phi }_{1}`et :math:`\tilde{\phi }{\rbrace _{i=2,4}={\phi }_{i}+{\phi }_{3}/2`.

4.4

L’élimination des degrés de liberté de contact en trop est discutée [§ ]), on choisit la substitution (en mettant 1 sur la diagonale et 0 dans le second membre).

4.3.2

On peut généraliser cette approche pour les fonds de fissure. En effet sur les deux exemples de la figure , l’arête {1-4} est coupée alors que {2-3} ne l’est pas: les nœuds 2 et 3 sont donc à éliminer. En utilisant l’approche décrite précédemment, les fonctions de forme modifiées s’écrivent \(\tilde{\phi }{\rbrace _{i=1,4}={\phi }_{i}+({\phi }_{2}+{\phi }_{3})/2\). Cela revient à faire l’intégration P0 décrite au paragraphe [§] sur les éléments coupées contenant la pointe.

Figure 6.5.3-2 : Les nœuds rouges n’appartiennent pas à une arête coupée.

Bibliographie

1. “Éléments de contacts dérivés d’une formulation hybride continue”, Documentation de Référence de Code_Aster n° [R5.03.52]

2. Ben Dhia H., Vautier I.,“Une formulation pour traiter le contact frottement en 3d dans le Code_Aster », rapport de recherche HI-75/99/007/A, Juin 1999, EDF

3. Alart P., Curnier A., “A mixed formulation for frictional contact problems prone to Newton like solution methods”, Computer Methods in Applied Mechanics and Engineering, vol. 92, pp. 353-375, 1991

4. 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

5. Pellet J., “Dualisation des conditions aux limites”, Documentation de Référence de Code_Aster n° [R3.03.01]

6. Ern A., Guermond J.L.,Theory and practice of finite elements, Springer, 2004

7. Brenner S.C., Scott L.R., The mathematical theory of finite element methods, 2nded., Springer, 2002

8. Laursen T.A., Simo J.C.“A continuum element-based formulation for the implicit solution of multi-body, large deformation frictional contact problem”, International Journal for Numerical Methods in Engineering, vol. 36, pp. 3451-3485, 1993

9. Wriggers P., “Finite element algorithms for contact problem”, Arch. Of Comp. Meth. In Eng., vol. 2, pp. 1-49, 1995

10. Curnier A, He, Q.C., Klarbring A.,“Continuum mechanics modelling of large deformation contact with friction”, Contact mechanics, ed. Plenum Press, 1995

11. Pietrzak G.,“continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems”, Thèse de doctorat, École Polytechnique Fédérale de Lausanne, 1997

12. Alart P., Barboteu M., «Éléments contact, méthode de Newton généralisée et décomposition de domaine» Problèmes non linéaires Appliqués, École CEA – EDF – INRIA 1999

13. Chapelle D., Bathe K. J., “The Inf-sup test”, Computers & Structures, vol. 47, pp. 537-545, 1993

14. Moës N., Béchet E., Tourbier M., “Imposing essential boundary conditions in the X-FEM”, International Journal for Numerical Methods in Engineering , 2006

15. DAUX C., MOES N., DOLBOW J., SUKUMAR N., BELYTSCHKO T., «Arbitrary branched and intersecting cracks with the extended finite element method», International Journal for Numerical Methods in Engineering, 48 (2000), 1741-1760

16. SIAVELIS M., «Modélisation numérique X-FEM de grands glissements avec frottement le longd’un réseau de discontinuités.» Thèse de doctorat, Ecole Centrale de Nantes, 2011

17. 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.

Description des versions

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