r5.03.50 Formulation discrète du contact-frottement

Résumé :

On décrit dans ce document les méthodes numériques utilisées pour traiter les problèmes de contact/frottement en grands déplacements dans l’opérateur STAT_NON_LINEou DYNA_NON_LINEà l’aide d’une formulation discrète. Les algorithmes disponibles permettent de traiter le problème de manière exacte (en contact) ou approchée (en contact et en frottement), sans choix restrictif sur le problème mécanique sous-jacent (cinématique ou lois de comportement non-linéaire).

Ce document contient aussi les éléments théoriques pour décrire la résolution des problèmes unilatéraux généraux (sans appariement comme le suintement par exemple) utilisables par le mot-clef FORMULATION=”LIAISON_UNIL”.

Table des matières

Mécanique du contact et du frottement

Définition du problème et notations

On considère deux solides pouvant entrer en contact frottant, la zone de contact est soit ponctuelle, soit linéique, soit surfacique. Soit \(\left\lbrace n\right\rbrace\) la normale sortante à la surface de l’un des solides en contact et \(\left\lbrace u\right\rbrace\) le vecteur déplacement entre les deux solides. Alors \(g=\langle u\rangle .\left\lbrace n\right\rbrace\) est le déplacement projeté sur cette normale, on l’appellera jeu . A partir de la donnée de la contrainte de Cauchy, on définit la pression [1]_ \(p=\langle n\rangle .\left[\sigma \right].\left\lbrace n\right\rbrace\) et la contrainte de cisaillement tangentiel \(\left\lbrace r\right\rbrace =\left[\sigma \right].\left\lbrace n\right\rbrace -p.\left\lbrace n\right\rbrace\) comme étant exercés par l’une des surfaces sur l’autre.

Figure 2.1-a : définition du repère local de contact.

La force de cisaillement a pour direction dans la zone de contact un vecteur \(\left\lbrace t\right\rbrace\) situé dans le plan tangent \((\left\lbrace {t}_{1}\right\rbrace ,\left\lbrace {t}_{2}\right\rbrace )\) indiqué sur la figure ci-dessus. L’équation (5000) définit la contrainte de cisaillement \(r\) exercée par le solide (5001) sur le solide (5000) par unité de surface de contact.

(2206)\[ \begin{align}\begin{aligned}\left\lbrace r\right\rbrace =\left[\sigma \right].\left\lbrace n\right\rbrace -(\langle n\rangle .\left[\sigma \right].\left\lbrace n\right\rbrace ).\left\lbrace n\right\rbrace =r.\left\lbrace t\right\rbrace\\\textrm{avec}\\r=\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }\end{aligned}\end{align} \]

Conditions de contact de Signorini

On a introduit les deux variables définissant le contact:

• La distance signée entre les deux surfaces ou «gap» \(g\) ;

• La pression de contact \(p\) ;

On définit alors les trois conditions de contact de Hertz-Signorini-Moreau

• Condition géométrique: impénétrabilité de la matière (Signorini-Hertz).

(2207)\[\begin{split}g \ge 0 \iff \left\{ \begin{array}{l} g > 0 \quad \text{pas de contact} \\ g = 0 \quad \text{contact} \end{array} \right.\end{split}\]

• Condition mécanique: intensilité (Signorini-Hertz)

(2208)\[\begin{split}p\ge 0\mathrm{\iff} \left\{ \begin{array}{c} p>0 \quad \text{contact} \\ p=0 \quad \text{pas de contact}\end{array} \right.\end{split}\]

• Condition énergétique: complémentarité/exclusion (Moreau)

(2209)\[\begin{split}p.g=0\iff \left\{ \begin{array}{c} p=0 \quad \text{décollement} \\ g=0 \quad \text{contact}\end{array} \right.\end{split}\]

Graphiquement, les trois conditions sont représentées sur la figure ci-dessous où les zones hachurées représentent les zones exclues et les zones blanches ou les traits noirs représentent les zones autorisées.

Impénétrabilité	Intensilité	Complémentarité
Figure 2.2-a : représentation graphique des conditions de Hertz-Signorini-Moreau.

En combinant les trois conditions on obtient le graphe de la figure ci-dessous.

Figure 2.2-b : graphe de la condition de contact unilatéral.

Le problème de contact ainsi posé introduit une relation non-univoque (\(p\) n’est pas une fonction de \(g\) ), semi-définie positive et non-différentiable en \(p=g=0\) . Il s’agit d’un problème mathématiquement difficile à traiter. Le contact est un phénomène réversible et conservatif pour lequel on peut introduire un potentiel énergétique et dont le résultat ne dépend pas du trajet de chargement, il est semblable au modèle élastoplastique de Hencky. La loi de contact de Signorini s’écrit:

(2210)\[\begin{split}\left\{ \begin{array}{cc}g\ge 0& (a)\\ p\ge 0& (b)\\ p.g=0& (c)\end{array} \right.\end{split}\]

Remarques:

• Le contact est supposé sans adhésion grâce à la deuxième condition (intensilité) .

• La troisième condition permet au problème de contact unilatéral d’être bien posé pour être soluble par des techniques classiques d’optimisation sous contraintes (introduction des conditions de Kuhn & Tucker) commela méthode des contraintes actives.

Formulation du problème de frottement

Définitions

Les critères de frottement choisis sont de la forme:

(2211)\[h(\left\lbrace r\right\rbrace )\le 0\]

où \(h(\left\lbrace r\right\rbrace )\) est une fonction convexe. Le domaine de non glissement est défini par l’intérieur du convexe. Deux critères de frottement de la forme \(h(\left\lbrace r\right\rbrace )\le 0\) sont particulièrement utilisés: le critère de Tresca et le critère de Coulomb.

Le critère de Tresca

Le critère de Tresca est défini par la fonction \(h(\left\lbrace r\right\rbrace )\) suivante:

(2212)\[ \begin{align}\begin{aligned}h(\left\lbrace r\right\rbrace )=\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }-k\le 0\\\textrm{avec } k \textrm{ une constante donnée}\end{aligned}\end{align} \]

On note \(C\) le disque convexe de rayon \(k\) centré à l’origine défini par:

(2213)\[C=\left\lbrace \left\lbrace r\right\rbrace \text{tel que}\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }\le k\right\rbrace\]

La condition de non glissement est alors définie par l’appartenance de \(\left\lbrace r\right\rbrace\) à l’intérieur du disque \(C\) . En cas de glissement, pour \(\left\lbrace r\right\rbrace\) situé sur la frontière de \(C\) , la direction de glissement \(\left\lbrace t\right\rbrace\) de \(\left\lbrace \dot{u}\right\rbrace\) est donnée par la normale au critère en \(\left\lbrace r\right\rbrace\) , comme indiqué sur la .

Figure 2.3.2-a : disque de frottement pour le critère de Tresca.

Ce qui se traduit par une expression de la vitesse de glissement analogue au multiplicateur plastique:

(2214)\[{\dot{u}}_{t}=\lambda .r\]

Avec \(\lambda \ge 0\) , c’est-à-dire que la vitesse de glissement est dans la même direction que la contrainte tangentielle.

Le critère de Coulomb

Le critère de Coulomb est défini par la fonction \(h(\left\lbrace r\right\rbrace )\) suivante:

(2215)\[ \begin{align}\begin{aligned}h(\left\lbrace r\right\rbrace ,\mu ,p)=\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }-k(\mu ,p)\le 0\\\textrm{et}\\k=\mu .\mid p\mid\end{aligned}\end{align} \]

La valeur de \(k\) dépend de la pression de contact \(p\) et du coefficient de frottement de Coulomb \(\mu\) Ainsi défini, \(h\) est un cône. En cas de glissement, pour \(\left\lbrace r\right\rbrace\) situé sur la frontière de \(h\) , la direction de glissement \(\left\lbrace t\right\rbrace\) de \(\dot{u}\) n’est pas donnée par la normale au critère en \(\left\lbrace r\right\rbrace\) , mais par la normale au disque convexe \(C\) de rayon \(k=\mu .\mid p\mid\) . Le critère de Tresca correspond à une «tranche» suivant le plan orthogonal au cône de Coulomb.

Figure 2.3.3-a : cône de frottement pour le critère de Coulomb.

Application au frottement de Coulomb

Écrivons le système d’équations et d’inéquations devant être vérifié par ces grandeurs dans le cas du critère de frottement de Coulomb:

(2216)\[\begin{split} \left\{ \begin{array}{cc}g\ge 0& (a)\\ p\ge 0& (b)\\ p.g=0& (c)\\ \mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }-\mu .\mid p\mid \le 0& (d)\\ \left\lbrace {\dot{u}}_{t}\right\rbrace =\lambda .\left\lbrace r\right\rbrace & (e)\\ \lambda .(\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }-\mu .\mid p\mid )=0& (f)\\ \lambda \ge 0& (g)\end{array} \right.\end{split}\]

Le premier jeu d’équations et d’inéquations ((4948) a-c) correspond à la gestion du contact. Le second lot ((4948) d-g) correspond à la description du frottement obéissant au critère de Coulomb. Il fait intervenir plusieurs champs et les lie entre eux: la pression normale \(p\) , la contrainte de cisaillement \(\left\lbrace r\right\rbrace\) et la vitesse tangente \(\left\lbrace {\dot{u}}_{t}\right\rbrace\) . Il peut se comprendre comme suit:

(2217)\[\begin{split} \left\{ \begin{array}{cc}\text{Si}\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }<\mu .\mid p\mid \text{alors}\lambda =0\text{et}\left\lbrace {\dot{u}}_{t}\right\rbrace =0& (a)\\ \text{Si}\mathrm{\parallel }\left\lbrace r\right\rbrace \mathrm{\parallel }=\mu .\mid p\mid \text{alors}\lambda >0\text{et}\left\lbrace {\dot{u}}_{t}\right\rbrace =\lambda .\left\lbrace r\right\rbrace & (b)\end{array} \right.\end{split}\]

Sur la figure ci-dessous on a représenté le cône de Coulomb. Dans l’espace des contraintes, l’effort de contact frottant ne peut se trouver qu’à l’intérieur du cône de Coulomb: s’il est strictement à l’intérieur, le contact est adhérent; s’il est sur la surface du cône, le contact est glissant.

Figure 2.3.4-a : interprétation du cône de frottement.

On peut donner une autre représentation de ce critère (voir figure ci-dessous).

Figure 2.3.4-b : graphe du frottement de Coulomb.

Le frottement induit la notion de seuil. La relation introduite par le frottement de Coulomb est une relation non-univoque non-différentiable. Contrairement au contact, la relation n’est pas associative [2]_ et ne peut pas dériver directement d’un potentiel de dissipation énergétique. On remarquera également que le frottement lie la vitesse relative des deux surfaces et non le déplacement. Toutefois on peut établir qu’en statique, si l’on résout le problème sous forme incrémentale (ce qui est le cas lorsqu’on utilise Newton comme dans code_aster ), on peut remplacer la vitesse relative par l’incrément de déplacement tangentiel.

Formulation par inclusions différentielles

Le caractère non-différentiable des relations de frottement nous amène à introduire la notion de sous-inclusion différentielle.

On note \(V\) l’ensemble des déplacements cinématiquement admissibles du problème. La relation entre la vitesse de glissement relatif \(\left\lbrace \dot{u}\right\rbrace\) et la contrainte de cisaillement \(r\) traduit les deux états possibles du système: non glissement ou glissement relatif suivant la direction normale en \(\left\lbrace r\right\rbrace\) au disque convexe \(C\) . Pour les trois critères présentés, la fonction \(\left\lbrace {\dot{u}}_{t}\right\rbrace (r)\) et sa réciproque \(\left\lbrace r({\dot{u}}_{t})\right\rbrace\) appartiennent toutes deux aux sous‑différentiels de deux pseudo-potentiels conjugués \({\Psi}_{c}(r)\) et \({\Psi}_{c}^{\text{*}}({\dot{u}}_{t})\) , de telle sorte que l’on peut écrire:

(2218)\[ \begin{align}\begin{aligned}\left\lbrace \dot{u}\right\rbrace \in \partial {\Psi}_{c}(r)\\\textrm{et}\\\left\lbrace r\right\rbrace \in \partial {\Psi}_{c}^{\text{*}}(\dot{u})\end{aligned}\end{align} \]

L’apparition des inclusions différentielles provient du caractère non-différentiable des lois de contact‑frottement. En effet, \({\Psi}_{c}\) désigne la fonction indicatrice du disque convexe \(C\) de rayon \(k\) , centré à l’origine, précédemment défini. Elle est telle que:

(2219)\[\begin{split}{\Psi}_{c}(r)=\lbrace \begin{array}{c}0\text{si}\left\lbrace r\right\rbrace \in C\\ +\infty \text{sinon}\end{array}\end{split}\]

Le sous-différentiel \(\partial {\Psi}_{c}(r)\) de la fonction \({\Psi}_{c}\) en \(r\) se confond avec la normale extérieure à \(C\) en \(\left\lbrace r\right\rbrace\) . \({\Psi}_{c}^{\text{*}}(\dot{u})=k\mathrm{\parallel }\dot{u}\mathrm{\parallel }\) , où \(k\) est le seuil de résistance au frottement, \({\Psi}_{c}^{\text{*}}\) est la conjuguée de Fenchel de la fonction indicatrice \({\Psi}_{c}\) . \({\Psi}_{c}^{\text{*}}\) est positivement homogène de degré un. Cette fonction s’interprète comme la densité de puissance dissipée dans le glissement. Utilisant les notions de sous-différentiel, on peut établir les relations suivantes pour \(\left\lbrace \dot{u}\right\rbrace\) et \(\left\lbrace r\right\rbrace\) associés:

(2220)\[\begin{split}\left\{ \begin{array}{c}\dot{u}\in \partial {\Psi}_{c}(r)\mathrm{\iff }\left\{ \begin{array}{c}{\Psi}_{c}(r)=0\\ \dot{u}.(r'-r)\le 0,\forall r'\in C\end{array}\right.\\ r\in \partial {\Psi}_{c}^{\text{*}}(\dot{u})\mathrm{\iff }\left\{ \begin{array}{c}\dot{u}\in V\\ r.(\dot{v}-\dot{u})\le {\Psi}_{c}^{\text{*}}(\dot{v})-{\Psi}_{c}^{\text{*}}(\dot{u}),\forall \dot{v}\in V\end{array}\right.\\ {\Psi}_{c}^{\text{*}}(\dot{u})+{\Psi}_{c}(r)=\dot{u}.r=k\mathrm{\parallel }\dot{u}\mathrm{\parallel }\end{array}\right.\end{split}\]

Remarques:

• Lesdeux pseudo-potentiels conjugués présentés sont non-différentiables.

• Une fois connue la réaction normale pour le critère de Coulomb, on se ramène localement à un critère de frottement de Tresca dont le seuil vaut \(k=\mu .|p|\) .

• Les critères locaux adoptés ayant une forme circulaire on en déduit que \(\dot{u}\in \partial {\Psi}_{c}(r)\) implique qu’il existe \(\lambda\) réel positif tel que \(\dot{u}=\lambda r\) .

• La formulation du problème en vitesse suggère une résolution numérique incrémentale du problème du frottement. La résolution du problème d’équilibre sera donc présentée sous forme incrémentale.

Résolution du problème d’équilibre

On considère deux solides de volume total \(\Omega ` dont la surface de contact est :math:`{\Gamma}_{c}\) . Pour simplifier, on supposera l’existence d’une énergie de déformation différentiable pour caractériser la réponse des deux solides séparés à des sollicitations externes (en fait, on peut démontrer que les résultats donnés ci-après sont indépendants de cette hypothèse). On note \(V\) l’ensemble des champs de déplacement cinématiquement admissibles, contraints par le respect des conditions de contact et de frottement sur l’interface. L’équilibre des deux solides en l’absence de frottement s’écrit:

Trouver U champ de déplacement cinématiquement admissible tel que

(2221)\[\left\lbrace U=\underset{v\in V}{\text{argmin}}\left[\Phi (\varepsilon (v))-W(v)\right]\right\rbrace \iff \left\lbrace \Phi (U)-W(U)\le \Phi (v)-W(v),\forall v\in V\right\rbrace\]

En élasticité, \(\Phi (v)={\int}_{\Omega}\phi (\epsilon (v)).d\Omega\) est l’énergie de déformation. La fonction \(W(v)\) représente le travail des forces externes. Une condition nécessaire (qui devient suffisante si \(\Phi\) est strictement convexe) pour que cet équilibre soit vérifié est que:

(2222)\[D\Phi (U)-\mathit{DW}\left(U\right)=D\Phi (U)-{L}^{\text{ext}}=0\]

où \(D\) est l’opérateur dérivée de Gâteaux et \({L}^{\text{ext}}\) est la forme linéaire associée aux forces externes.

Avec l’introduction du frottement, le problème doit être abordé sous forme incrémentale. On est conduit (voir [bib2] et [bib3] au problème de minimisation suivant sur l’ensemble \(\overline{V}\) des champs cinématiquement admissibles contraints par le respect des conditions de contact et de frottement sur l’interface:

U connu, trouver \(\Delta U\in \overline{V}\) tel que

(2223)\[U+\Delta U=\underset{\Delta v\in \stackrel{ˉ}{V}}{\text{argmin}}\left[\Phi (\varepsilon (U+\Delta v))+{\Psi}_{c}^{\text{*}}(\Delta {v}_{t})-W(U+\Delta v)\right]\]

\(U+\Delta U\) est ainsi solution de:

(2224)\[\underset{\Delta v\in \stackrel{ˉ}{V}}{\min}\left[\underset{\Omega}{\int}\varphi (\varepsilon (U+\Delta v))d\Omega +\underset{{\Gamma}_{c}}{\int}k\mid \Delta {v}_{t}\mid d{\Gamma}_{c}-W(U+\Delta v)\right]\]

où \(\Delta {v}_{t}\) est la composante tangentielle de l’incrément de déplacement relatif [3] du solide \(2\) par rapport au solide \(1\) le long de la surface de contact, avec les conventions adoptées au r5.03.50-formulation-frottement .

En utilisant les relations \({\Psi}_{c}^{\text{*}}(\Delta {v}_{t})=k\mid \Delta {v}_{t}\mid ` et :math:`{\Psi}_{c}^{\text{*}}(\Delta {v}_{t})\ge r\Delta v\) si \(r\in C\) on en déduit que \(U+\Delta U\) est solution du problème de \(\text{minmax}\) suivant, sur l’espace \(V\) des champs cinématiquement admissibles :

(2225)\[\underset{\Delta v\in V}{\text{Min}}\underset{r}{\text{Max}}\text{}J(U+\Delta v,r)\]

où la fonctionnelle \(J\) vaut:

(2226)\[J(U+\Delta v,r)=\underset{\Omega}{\int}\phi (\varepsilon (U+\Delta v)).d\Omega +\underset{{\Gamma}_{c}}{\int}(r.\Delta {v}_{t}-{\Psi}_{c}(r)).d{\Gamma}_{c}-W(U+\Delta v)\]

La présence de la fonction indicatrice dans cette expression indique que le cisaillement \(r\) sur la surface de contact \({\Gamma}_{c}\) appartient au disque convexe de frottement \(C\) .

Formulation variationnelle

Si \(\phi\) est convexe, le problème de \(\text{minmax}\) à résoudre se met de façon équivalente sous la forme:

Trouver \(\Delta U\in V\) et \(r\in C\) , ensemble de variables indépendantes tel que

(2227)\[\delta J(U+\Delta U,r)\ni 0\]

Ceci revient à résoudre le système d’équations à l’équilibre suivant:

(2228)\[\begin{split}\left\{ \begin{array}{cc}\underset{\Omega}{\int}\frac{\partial \varphi }{\partial \varepsilon }(\varepsilon (U+\Delta U)).\delta \varepsilon .d\Omega +\underset{{\Gamma}_{c}}{\int}r.\delta {v}_{t}.d{\Gamma}_{c}-{L}^{\text{ext}}.\delta v=0& (a)\\ \underset{{\Gamma}_{c}}{\int}(\delta r.\Delta {U}_{t}.d{\Gamma}_{c}-\partial {\Psi}_{c}(r))\ni 0& (b)\end{array}\right.\end{split}\]

ou encore de manière équivalente:

(2229)\[\begin{split}\left\{ \begin{array}{cc}\underset{\Omega}{\int}\frac{\partial \varphi }{\partial \varepsilon }(\varepsilon (U+\Delta U)).\delta \varepsilon .d\Omega +\underset{{\Gamma}_{c}}{\int}r.\delta {v}_{t}.d{\Gamma}_{c}-{L}^{\text{ext}}.\delta v=0& (a)\\ \Delta {U}_{t}\in \partial {\Psi}_{c}(r)& (b)\\ r=({\sigma}_{1}.n).t\text{sur}{\Gamma}_{c}& (c)\end{array}\right.\end{split}\]

Comme dans la section précédente, \({L}^{\text{ext}}\) est la forme linéaire associée aux forces externes. La forme linéaire \({L}^{\text{frot}}\) est associée aux forces de cisaillement exercées par le solide \(2\) sur la surface de contact du solide \(1\) . On notera aussi que la formulation variationnelle permet de retrouver non seulement les équations d’équilibre du système mais aussi l’appartenance de \(\Delta {U}_{t}\) au sous-différentiel de \({\Psi}_{c}\) .

Étapes de la résolution du problème de contact-frottement

La méthode de résolution discrète implantée dans code_aster est fondée sur une écriture des relations d’interpénétration sur les nœuds des maillages en vis-à-vis, ce qui implique:

1. Une description discrète des surfaces de contact (maillage);

2. La recherche de la distance minimale de projection et de la position de cette projection (opération dite d’appariement);

3. L’écriture des relations cinématiques entre les nœuds;

4. Des algorithmes de résolution du problème de contact/frottement;

Nous allons détailler chacun de ces aspects dans la suite du document.

Pour les formulations dites «discrètes» (par opposition à la formulation dite «continue», voir [r5.03.52]), code_aster résout le problème de contact par des méthodes qui appartiennent à la famille des méthodes appelées «méthode des statuts» dans la littérature, avec un découplage de type global-local.

Les algorithmes de contact/frottement agissent en deux temps:

1. Problème géométrique d’appariement

• Repérage: définition des surfaces potentielles de contact (cf. r5.03.50-definition-zones-potentielles) par l’utilisateur [4].

• Appariement: détermination des couples potentiels de contact (cf. r5.03.50-algorithme-appariement) par une méthode de recherche du jeu minimum entre un nœud et une facette (pour l’appariement nœud/facette) ou du nœud le plus proche d’un autre nœud (pour l’appariement nodal).

• Cinématique: écriture de la relation de non pénétration par la détermination de la direction de projection et l’évaluation des coefficients (cf. r5.03.50-relations-cinematiques). La relation est écrite entre le nœud esclave et les nœuds maîtres.

2. Problème mécanique. Plusieurs types d’algorithmes sont implantés dans code_aster :

• Un algorithme basé sur la méthode des contraintes actives [bib1] utilisable en contact sans frottement uniquement. Cet algorithme repose sur les hypothèses que la matrice tangente du problème est symétrique et définie positive, et garantit la convergence en un nombre fini d’itérations dans ce cas [bib1], [bib15], [bib16]. Il ne doit pas être utilisé en dehors de ce cadre, c’est à dire en particulier pour les matrices non-symétrique, par exemple en modélisations THM ou en grandes déformations. C’est celui qui est utilisé par défaut et qui correspond à ALGO_CONT=”CONTRAINTE”.

• Une variante de la méthode des contraintes actives à l’aide d’une résolution itérative par gradient conjugué projeté. ALGO_CONT=”GCP”. Méthode réservée exclusivement aux problèmes de contact sans frottement.

• Un algorithme de résolution par régularisation des conditions de contact et/ou de frottement, que l’on active avec ALGO_CONT=”PENALISATION’et ALGO_FROT=”PENALISATION’en choisissant des coefficients de pénalisation judicieux. Ceciimplique a fortiori une étude paramétrique sur la valeur de ces coefficients: ils ne sont pas homogènes à \(\left\lbrace u\right\rbrace\) , mais de dimension d’un rapport d’une force par un déplacement, et il est difficile de donner une règle générale sur le choix de leurs valeurs. Il est fortement conseillé à l’utilisateur de vérifier a posteriori que la condition unilatérale est correctement prise en compte en visualisant les résultats dans la zone concernée. Si ce n’est pas le cas, il est nécessaire d’augmenter le coefficient de pénalisation de contact jusqu’à ce que cela le soit. Cet algorithme de régularisation peut être choisi dans le cas d’une modélisation impliquant une matrice non symétrique, et ne pouvant donc pas être traité par l’algorithme de contrainte actives.

3. Conditions unilatérales. Plus généralement, il est possible d’utiliser les algorithmes ALGO_CONT=”CONTRAINTE” et ALGO_CONT=”PENALISATION” pour imposer des conditions aux limites d’inégalités portant sur n’importe quel degré de liberté. Celles-ci sont définies par FORMULATION =”LIAISON_UNIL” dans DEFI_CONTACT (voir la documentation [u4.44.11]), et prenne la forme algébrique \(\left[A\right]\left\lbrace u\right\rbrace \le \left\lbrace d\right\rbrace\) . Au lieu d’être construits par appariement, \(\left[A\right]\) et \(\left\lbrace d\right\rbrace\) sont fixes et construits à partir des informations données par l’utilisateur dans DEFI_CONTACT, mais le reste du traitement est sensiblement identique.

[1]

Le terme « pression » est une abréviation abusive de « densité d’effort de contact »

[2]

Par contre le frottement de Tresca est bien une loi associative.

[3]

\(\Delta {v}_{t}\) est la composante tangentielle de l’incrément de déplacement relatif en formulation quasi-statique, elle devient la vitesse relative des deux surfaces en formulation dynamique.

[4]

Il n’y pas de mécanisme automatique de repérage des surfaces en contact dans code_aster

Problème géométrique d’appariement

Introduction

Le premier problème à traiter dans le cas du contact est de définir de manière fiable et robuste quelles parties des deux surfaces risquent d’entrer en contact. C’est l’opération dite d’appariement, commune avec la formulation hybride mixte dite méthode «continue» présentées dans le document [r5.03.52].

Il y a deux types d’appariement disponibles dans code_aster :

• L’appariement nodal;

• L’appariement maître-esclave (ou nœud-facette).

L’appariement nodal (APPARIEMENT=”NODAL”) impose que le déplacement relatif entre un nœud esclave et le nœud maître qui lui est apparié, projeté sur la direction de la normale au nœud esclave, soit inférieur au jeu initial dans cette direction. L’utilisation de cette formulation est déconseillée car elle nécessite d’avoir des maillages compatibles (nœuds «en face») qui le restent au cours de la déformation (hypothèse de petits glissements), et pour lesquels les normales maître et esclave sont à peu près colinéaires. Sans ces hypothèses, l’approximation faite devient hasardeuse et il est préférable d’utiliser la formulation nœud-facette.

L’appariement maître-esclave, choisi par le mot-clé APPARIEMENT=”MAIT_ESCL”, n’accorde pas un rôle équivalent aux deux surfaces: la surface (\({S}_{1}\) ) décrite sous GROUP_MA_MAITou MAILLE_MAIT est appelée la surface maître et la surface (\({S}_{2}\) ) est la surface esclave. Les conditions de non interpénétration expriment que les nœuds de la surface esclave (des étoiles sur la figure ci-dessous) ne pénètrent pas dans les mailles de la surface maître. On peut voir, par contre, qu’il est possible que les nœuds maîtres (des ronds) pénètrent dans la surface esclave.

Figure 3.1-a : surface maître et surface esclave.

Remarque:

• Les nœuds esclaves sont par défaut tous les nœuds appartenant aux mailles de contact définissant la surface esclave. Les mots-clés SANS_NOet SANS_GROUP_NOpermettent de donner, zone par zone, une liste des nœuds qui doivent être enlevés de la liste des nœuds esclaves. Cela permet d’enlever les nœuds soumis à des conditions aux limites de Dirichlet incompatibles avec le contact (voir r5.03.50-compatibilite-Dirichlet ).

Définition des zones potentielles de contact

code_aster ne disposant pas de mécanisme de détection automatique des zones potentielles de contact, c’est donc à l’utilisateur de définir a priori les zones dont il prédit qu’elles vont entrer en contact. Ces zones doivent être suffisamment larges pour ne pas observer d’interpénétration. Il faut noter que la phase d’appariement est une opération peu coûteuse en général (bien moins coûteuse que la résolution du problème mécanique de contact) et que l’on peut donc définir des zones suffisamment étendues sans risquer de pénaliser les performances, sous réserve de respecter certaines précautions pour assurer, en particulier, l’unicité des projections.

On considère les trois solides de la figure ci-dessous, représentés en 2D. On a défini trois zones possibles d’interpénétration entre les solides: une zone entre le solide \(A\) et le solide \(B\) , et deux zones entre le solide \(B\) et le solide \(C\) . L’utilisateur, qui définit ces zones dans le fichier de commande, suppose ici qu’en dehors de ces zones, il n’y a pas de risque d’interpénétration, compte tenu du chargement.

Figure 3.2-a : définition de trois zones de contact.

Chaque zone de contact est définie dans l’opérateur DEFI_CONTACT, mot-clef facteur ZONE. Une zone se compose par définition de deux surfaces dont on cherche à empêcher l’interpénétration: la première est définie sous le mot-clé GROUP_MA_MAIT(ou MAILLE_MAIT), la seconde sous le mot-clé GROUP_MA_ESCL(ou MAILLE_ESCL), c’est-à-dire par la donnée des mailles de bord qui les constituent. Ces mailles sont des SEG2ou des SEG3pour un maillage 2D, des TRIA3, TRIA6, QUAD4, QUAD8ou QUAD9pour un maillage 3D.

Remarques:

• Les mailles de bord nécessaires au contact ne seront pas créées par code_aster à partir des éléments volumiques et doivent donc déjà exister dans le fichier de maillage.

• Le choix des surfaces qui seront maîtres ou esclaves est important. Des informations sur ce sujet sont disponibles dans la documentation [u2.04.04].

Cas particulier du contact pour un câble ou une poutre en 3D

Il est possible en 3D de traiter le contact entre une maille SEG2ou SEG3(modélisant un câble ou une poutre) et une surface. Dans ce cas, il faut impérativement utiliser la méthode d’appariement “MAIT_ESCL’et donner les segments sous le mot-clé GROUP_MA_ESCL(mailles esclaves). La section de la poutre peut être alors prise en compte par l’utilisation du mot-clé DIST_ESCL(cf. r5.03.50-introduction-jeu-fictif).

Cas de l’appariement nodal

On doit choisir de prendre comme surface esclave celle qui comporte le moins de nœuds (un message d’erreur vous arrêtera si votre surface maître contient moins de nœuds que votre surface esclave), afin de maximiser les chances d’avoir un appariement injectif (un nœud maître n’est apparié qu’à un seul nœud esclave). Le nœud maître apparié à chaque nœud esclave est déterminé par un calcul de plus proche voisin expliqué dans le paragraphe r5.03.50-recherche-proche-voisin . On utilise la normale au nœud maître pour écrire la relation de non interpénétration.

Même dans le cas de l’appariement nodal, les surfaces de contact sont définies en termes de mailles. Les nœuds esclaves et maîtres sont alors les nœuds des mailles ainsi définies.

Orientation des normales

Il est impératif que les mailles de contact soient définies de façon à ce que la normale soit sortante : la connectivité des segments doit être définie dans l’ordre \(\mathrm{AB}\) , celle des triangles dans l’ordre \(\mathrm{ABC}\) , et celle des quadrangles dans l’ordre \(\mathrm{ABCD}\) , comme indiqué sur la figure . Pour une meilleure lecture du dessin, on a ici un peu écarté la maille de bord servant au contact de la «face» de l’élément volumique 2D ou 3D sur lequel elle s’appuie. On pourra utiliser pour cela l’opérateur MODI_MAILLAGE, options ORIE_PEAU_2D, ORIE_PEAU_3D ou ORIE_COQUE.

Figure 3.3-a : numérotation des mailles de contact pour avoir une normale sortante.

Algorithme d’appariement

Principe

L’algorithme d’appariement procède en deux temps:

1. Pour chaque nœud esclave on recherche le nœud maître le plus proche;

2. On cherche ensuite la maille maître attachée au nœud maître précédemment déterminé qui soit la plus proche du nœud esclave;

Dans le cas de l’appariement nodal, on ne fait évidemment que la première étape.

Recherche du plus proche voisin d’un nœud

Cette phase est commune aux deux formulations: nœud/facette et nœud/nœud. La méthode utilisée pour rechercher le nœud maître le plus proche d’un nœud esclave est très simple: il suffit de calculer la distance (en géométrie actualisée, cf. r5.03.50-reactualisation-geometrique ) entre le nœud esclave et les nœuds maîtres candidats.

Recherche de la maille maître appariée

L’algorithme pour rechercher la maille maître qui est appariée au nœud esclave est le suivant:

• Connaissant le nœud maître le plus proche du nœud esclave \(P\) (voir r5.03.50-recherche-proche-voisin ), on examine successivement les mailles maîtres contenant ce nœud.

• Pour chaque maille ainsi repérée, on cherche le projeté orthogonal \(M\) du nœud esclave \(P\) sur la maille maître.

• La maille minimisant la distance \(\mathit{PM}\) est choisie pour être appariée au nœud esclave.

Remarques:

• Si le nœud esclave se projette en dehors d’une maille maître, on rabat sa projection sur cette maille sous certaines conditions (voir r5.03.50-traitement-projections ).

• Si au moins une projection a lieu à l’intérieur d’une maille maître, alors celle-ci sera préférée à une projection conduisant à un rabattement (quelle que soit la distance mesurée)

• Si toutes les projections du nœud esclave ont lieu en dehors de leur maille maître respective, le nœud esclave est considéré comme non apparié et sera donc exclus de la phase de résolution

Calcul de la projection du nœud esclave sur la maille maître

Notations

Les parties des surfaces \(\partial {\Omega}^{i}\) susceptibles d’entrer en contact lors de la déformation des deux solides sont notées \({\gamma}_{c}^{i}\) . On suppose l’existence de cartes régulières notées \({\Theta}^{i}\) décrivant les surfaces \({\gamma}_{c}^{i}\) . Ces cartes sont définies comme suit:

(2230)\[\begin{split}\begin{array}{cc}{\Theta}^{i}:& {\omega}^{i}\to \overline{{\Omega}^{i}}\\ & ({\xi}_{1}^{i},{\xi}_{2}^{i})\to {p}^{i}={\Theta}^{i}({\xi}_{1}^{i},{\xi}_{2}^{i})\end{array}\end{split}\]

où \({\omega}^{i}\) est un domaine borné (de référence) contenu dans \({ℝ}^{2}\) (c’est l’espace paramétrique de référence de l’élément fini). Par ailleurs, on désigne par \({\varphi}^{i}\) la transformation du solide \({B}^{i}\) , définie par:

(2231)\[\begin{split}\begin{array}{cc}{\varphi}^{i}:& \overline{{\Omega}^{i}}\to \overline{{\Omega}_{t}^{i}}\\ & \left\lbrace {p}^{i}\right\rbrace \to \left\lbrace {x}^{i}\right\rbrace \end{array}\end{split}\]

On se place toujours à un temps \(t\) fixé. La surface maître est notée \({\gamma}_{c}^{1}\) et la surface esclave est notée \({\gamma}_{c}^{2}\) . On effectue l’appariement en recherchant, pour tout point \(\left\lbrace {x}^{i}\right\rbrace ={\varphi}^{i}({p}^{i},t)\) de la frontière \({\gamma}_{c}^{1}\) le point \(\left\lbrace \tilde{x}\right\rbrace\) de \({\gamma}_{c}^{2}\) le plus proche. Cela revient à résoudre le problème d’optimisation sous contraintes suivant. Pour tout point \(\left\lbrace {x}^{1}\right\rbrace ={\varphi}^{1}({\Theta}^{1}({\xi}_{1}^{1},{\xi}_{2}^{1}),t)\) avec \(({\xi}_{1}^{1},{\xi}_{2}^{1})\in {\omega}^{1}\) , et pour tout \(t\ge 0\) , trouver \(\tilde{{\zeta}_{t}}=({\xi}_{1}^{2},{\xi}_{2}^{2})\in {\omega}^{2}\) tel que:

(2232)\[{\tilde{\zeta}}_{t}=\underset{({\xi}_{1}^{2},{\xi}_{2}^{2})\in {\omega}^{2}}{\text{argmin}}\left\lbrace \frac{1}{2}.∥{\varphi}^{1}({\Theta}^{1}({\xi}_{1}^{1},{\xi}_{2}^{1}),t)-{\varphi}^{2}({\Theta}^{2}({\xi}_{1}^{2},{\xi}_{2}^{2}),t)∥{2}^{}\right\rbrace\]

La solution \(\tilde{{\zeta}_{t}}\) est la position dans l’espace de référence paramétrique de la projection \(M\) du nœud esclave \(P\) sur la maille maître.

Formulation du problème de minimisation

Pour résoudre le problème non-linéaire (4854), on propose d’utiliser un algorithme de minimisation de Newton. Le problème est reformulé de la manière suivante: pour un nœud esclave donné on cherche le point (en coordonnées paramétriques) qui minimise la distance à une maille maître donnée (celle définie au paragraphe r5.03.50-recherche-maille-maitre ).

On note \(f\) la norme \({L}_{2}\) de la distance entre un nœud esclave situé en \(\left\lbrace {x}_{e}\right\rbrace\) et un point \(\left\lbrace {x}_{m}\right\rbrace\) appartenant à la maille maître en vis-à-vis:

(2233)\[f(\left\lbrace {x}_{m}\right\rbrace )={∥\left\lbrace {x}_{e}\right\rbrace -\left\lbrace {x}_{m}\right\rbrace ∥}^{2}\]

On écrit la fonctionnelle dans l’espace paramétrique, \(\left\lbrace {x}_{m}\right\rbrace\) s’exprimant à partir des \({\mathit{nb}}_{\mathit{no}}\) nœuds \(\left\lbrace {x}_{i,m}\right\rbrace\) de la maille maître appariée:

(2234)\[f({\xi}_{1},{\xi}_{2})=f(\left\lbrace \zeta \right\rbrace )={∥\left\lbrace {x}_{e}\right\rbrace -\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}{\phi }_{i}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ∥}^{2}\]

C’est cette quantité qu’il faut minimiser. On applique l’algorithme de Newton aux conditions de stationnarité (Euler-Lagrange) du problème de minimisation ainsi énoncé, c’est-à-dire à la fonctionnelle vectorielle \(\left\lbrace \nabla f(\left\lbrace \zeta \right\rbrace )\right\rbrace\) , expression du développement de Taylor au premier ordre, autour du point \(\left\lbrace {\zeta}_{0}\right\rbrace\) :

(2235)\[\left\lbrace \nabla f(\left\lbrace \zeta \right\rbrace )\right\rbrace \underset{\text{linéarisation}}{=}\left\lbrace \nabla f(\left\lbrace {\zeta}_{0}\right\rbrace )\right\rbrace +\left[H({\zeta}_{0})\right].(\left\lbrace \zeta \right\rbrace -\left\lbrace {\zeta}_{0}\right\rbrace )\]

Où \(\left[H\right]\) est la matrice Hessienne des dérivées secondes de \(f\) et \(\left\lbrace \nabla f(\left\lbrace \zeta \right\rbrace )\right\rbrace\) est son gradient. Un minimum de \(f\) survient lorsque son gradient est nul.

L’algorithme itératif s’écrit alors:

1. Partir du point initial \(\left\lbrace {\zeta}_{i=0}\right\rbrace\) sur la maille maître. Ce point de départ est simplement choisi en \(\left\lbrace {\zeta}_{i=0}\right\rbrace =\left\lbrace 0\right\rbrace\) ;

2. Évaluer le gradient \(\left\lbrace \nabla f\left(\left\lbrace {\zeta}_{i}\right\rbrace \right)\right\rbrace\) et la matrice Hessienne \(\left[H\left({\zeta}_{i}\right)\right]\) en ce point;

3. Calculer la direction de descente \(\left\lbrace {d\zeta }_{i}\right\rbrace =-{\left[H\left({\zeta}_{i}\right)\right]}^{-1}.\left\lbrace \nabla f\left(\left\lbrace {\zeta}_{i}\right\rbrace \right)\right\rbrace\)

4. Calculer le paramètre de recherche linéaire \(\alpha\) (voir r5.03.50-recherche-linéaire-algorithme-projection )

5. Calculer le point suivant tel que \(\left\lbrace {\zeta}_{i+1}\right\rbrace =\left\lbrace {\zeta}_{i}\right\rbrace +\alpha \left\lbrace {d\zeta }_{i}\right\rbrace\) ;

6. Si le processus a convergé, on s’arrête, sinon on boucle en 2.

Remarques:

• Ce problème s’écrit sans contraintes;

• On ne peut pas garantir l’unicité de la solution, le processus de Newton trouvera le premier point réalisant les conditions de stationnarité.

• L’existence du gradient et de la matrice Hessienne est assurée si la maille est suffisamment régulière, ce qui est toujours le cas sur un maillage élément fini pas trop distordu.

• Le processus itératif s’arrête par un critère sur l’incrément de déplacement dans l’espace paramétrique \(\epsilon =\frac{\sqrt{\langle {\zeta}_{i+1}\rangle \left\lbrace {\zeta}_{i+1}\right\rbrace -\langle {\zeta}_{i}\rangle \left\lbrace {\zeta}_{i}\right\rbrace }}{\sqrt{\langle {\zeta}_{i+1}\rangle \left\lbrace {\zeta}_{i+1}\right\rbrace }}\) avec \(\epsilon \le {10}^{-4}\) (valeur non modifiable par l’utilisateur). Cette valeur étant estimée par rapport aux coordonnées paramétriques, il n’y a pas de problèmes avec la taille des mailles et les unités utilisées.

• Le nombre maximum d’itérations de Newton est fixé à \(200\) (valeur non modifiable par l’utilisateur).

• Si le nombre maximum d’itérations est atteint, on sélectionne comme projeté l’itéré qui a minimisé la distance entre le nœud esclave et la maille maître (ce n’est donc pas parfaitement le projeté orthogonal).

• Il existe une variante dans laquelle la direction de recherche n’est pas estimée par l’algorithme de Newton mais une direction fixe donnée par l’utilisateur, ce qui peut être utile dans certains cas difficiles (maille parfaitement convexe qui ne donne pas une projection unique par exemple). On l’active par le mot-clef TYPE_APPA=”FIXE’et le vecteur est donné par DIRE_APPA.

• La projection sur les éléments segments en 3D (cas des poutres ou des câbles) nécessite une définition de la normale par l’utilisateur via les mots-clefs VECT_MAIT/VECT_ESCL.

Expression de la matrice tangente et du résidu

On donne ci-dessous les expressions des termes intervenants dans l’écriture de l’algorithme de Newton décrit ci-dessus. On rappelle la fonctionnelle:

(2236)\[f\left({\xi}_{1},{\xi}_{2}\right)=f\left(\left\lbrace \zeta \right\rbrace \right)={∥\left\lbrace {x}_{e}\right\rbrace -\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}{\varphi}_{i}\left(\left\lbrace \zeta \right\rbrace \right).\left\lbrace {x}_{i,m}\right\rbrace ∥}^{2}\]

Le gradient de cette fonctionnelle:

(2237)\[\begin{split}\left\lbrace \nabla f\left(\left\lbrace \zeta \right\rbrace \right)\right\rbrace =\left\lbrace \begin{array}{c}\frac{\partial f}{\partial {\xi}_{1}}\\ \frac{\partial f}{\partial {\xi}_{2}}\end{array}\right\rbrace\end{split}\]

Les termes s’écrivent explicitement:

(2238)\[\frac{\partial f}{\partial {\xi}_{j}}({\xi}_{1},{\xi}_{2})=2(\left\lbrace {x}_{e}\right\rbrace -\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}{\phi }_{i}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ).(-\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{\partial {\phi }_{i}}{\partial {\xi}_{j}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace )\]

Le hessien de cette fonctionnelle:

(2239)\[\begin{split}\left[H(\zeta )\right]=\left[\begin{array}{cc}\frac{{\partial}^{2}f}{\partial {\xi}_{1}^{2}}& \frac{{\partial}^{2}f}{\partial {\xi}_{1}.\partial {\xi}_{2}}\\ \frac{{\partial}^{2}f}{\partial {\xi}_{1}.\partial {\xi}_{2}}& \frac{{\partial}^{2}f}{\partial {\xi}_{2}^{2}}\end{array}\right]\end{split}\]

Les termes s’écrivent explicitement:

(2240)\[\frac{{\partial}^{2}f}{\partial {\xi}_{1}^{2}}({\xi}_{1},{\xi}_{2})={∥\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{\partial {\phi }_{i}}{\partial {\xi}_{j}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ∥}^{2}+2(\left\lbrace {x}_{e}\right\rbrace -\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}{\phi }_{i}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ).(-\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{{\partial}^{2}{\phi }_{i}}{\partial {\xi}_{j}^{2}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace )\]

Et pour les termes croisés:

(2241)\[\begin{split}\begin{array}{cc}\frac{{\partial}^{2}f}{\partial {\xi}_{j}\partial {\xi}_{k}}({\xi}_{1},{\xi}_{2})& =(\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{\partial {\phi }_{i}}{\partial {\xi}_{j}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ).(\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{\partial {\phi }_{i}}{\partial {\xi}_{k}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace )\\ & +2(\left\lbrace {x}_{e}\right\rbrace -\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}{\phi }_{i}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace ).(-\sum_{i=1}^{{\mathit{nb}}_{\mathit{no}}}\frac{{\partial}^{2}{\phi }_{i}}{\partial {\xi}_{j}.\partial {\xi}_{k}}(\left\lbrace \zeta \right\rbrace ).\left\lbrace {x}_{i,m}\right\rbrace )\end{array}\end{split}\]

Recherche linéaire dans l’algorithme de projection

Pour améliorer la robustesse de l’algorithme de projection, une phase de recherche linéaire est ajoutée dans l’algorithme de Newton (cf. r5.03.50-formulation-problème-minimisation ). Il s’agit, étant donnée une direction de descente \(\left\lbrace {d\zeta }_{i}\right\rbrace\) , de déterminer un paramètre d’avancée \(\alpha\) qui minimise une fonctionnelle \(F\) associée à \(f\) .

L’algorithme de recherche linéaire utilisé (recherche linéaire avec rebroussement et règle d’Armijo) est le même que pour l’intégration implicite des lois de comportement [r5.03.14].

La fonctionnelle \(F\) , associée à la fonctionnelle distance \(f\) définie en , est la suivante:

(2242)\[F(\alpha )=\frac{1}{2}{∥\nabla f\left(\left\lbrace \zeta \right\rbrace +\alpha \left\lbrace {d\zeta }_{i}\right\rbrace \right)∥}^{2}\]

Pour choisir \(\alpha\) , on ne va en fait pas chercher à minimiser exactement la fonctionnelle \(F\) mais une approximation de cette celle-ci (quadratique, cubique). Pour savoir quelle approximation utiliser, on s’appuie sur la règle d’Armijo. Pour plus de détails, on renvoie à [r5.03.14] et [bib14] qui détaillent la mise en œuvre d’un tel algorithme.

Pour la recherche linéaire implémentée dans l’algorithme de projection, les paramètres choisis (en dur) sont les suivants:

Paramètre de la règle d’Armijo	\(\omega =0,0001\)
Borne min pour le rabattement	\({\alpha}_{min}=0,1\)
Borne max pour le rabattement	\({\alpha}_{max}=0,5\)
Nombre maximum d’interpolations cubiques	\({\mathit{it}}_{max}=2\)

Ces paramètres impliquent notamment que si le pas par défaut (\(\alpha =1\) ) ne satisfait pas la règle d’Armijo, c’est-à-dire que la direction de descente choisie ne rapproche pas du projeté orthogonal, alors \(\alpha\) sera au plus égal à \(0,5\) . Ce choix permet de favoriser la robustesse (éviter de suivre une direction de recherche erronée) plutôt que la performance (valeurs de \(\alpha\) supérieures à \(1\) ).

Pour des maillages de mauvaise qualité, avec des mailles distordues, la recherche linéaire a montré son intérêt, trouvant le bon projeté orthogonal sur des cas difficiles là où un algorithme sans recherche linéaire échouait.

Traitement des projections hors maille

Il existe des projections dont le résultat est sensible à des paramètres purement numériques ou dont l’existence et l’unicité mathématiques ne sont pas garanties. Dans certaines conditions, on peut détecter du contact entre deux surfaces alors qu’il y en a pas. Le problème vient d’abord d’une définition incorrecte et imparfaite des surfaces susceptibles d’entrer en contact. Prenons le cas du contact en 2D (les surfaces de contact sont donc des segments) où un nœud esclave doit se projeter en dehors de la surface maître (voir figure ):

Figure 3.4.5-a : projection hors d’une maille maître.

L’utilisateur peut choisir de «re-projeter» ce nœud esclave sur la maille maître prolongée ou de ne pas le faire (voir figure ).

Figure 3.4.5-b : principe de la zone de tolérance pour la projection sur une maille maître.

La valeur limite de cette re-projection est fixée par le mot-clef TOLE_PROJ_EXT qui prend pour argument la valeur (rapportée à l’élément de référence) de l’extension de la maille maître dans laquelle on autorise la re-projection. Par défaut, cette valeur est fixée à \(0.50\) . Ce qui signifie que tout nœud esclave se projetant à une distance supérieure à la moitié de la longueur de la maille maître ne sera pas re-projeté. Pour interdire complètement la re-projection, il suffit de fixer TOLE_PROJ_EXT à zéro. Cet opérateur est valable en 2D et en 3D. Dans le cas 3D, il s’agit de l’extension d’une maille surfacique de contact, et, comme on raisonne dans l’espace paramétrique, la courbure des bords des éléments est bien prise en compte.

De même si les surfaces de contact sont trop étendues, tout nœud esclave situé derrière la maille maître est apparié. La notion de «devant-derrière» est donnée par l’orientation des normales donnée par l’utilisateur (voir r5.03.50-orientation-normales ). On peut restreindre cette recherche de maille appariée avec le paramètre TOLE_APPA qui spécifie la distance maximale de recherche des mailles «appariables» avec le nœud esclave.

Cas de la projection dans une direction donnée

Lorsqu’on utilise l’option DIRE_APPA, l’utilisateur peut forcer l’algorithme de projection non pas à chercher la distance minimale entre un nœud esclave et sa projection orthogonale sur la surface maître, mais à simplement trouver le point de la surface maître lié au nœud esclave par un vecteur donné. Dans ce cas l’algorithme est le suivant:

1. Partir du point initial \(\left\lbrace {\zeta}_{i=0}\right\rbrace\) sur la maille maître. Ce point de départ est simplement choisi en \(\left\lbrace {\zeta}_{i=0}\right\rbrace =\left\lbrace 0\right\rbrace\) ;

2. Construire la matrice de Newton \(\left[H({\zeta}_{i})\right]\) par l’utilisation de la direction donné par DIRE_APPA;

3. Construire le second membre de Newton \(\left\lbrace R\left(\left\lbrace {\zeta}_{i}\right\rbrace \right)\right\rbrace\) par la différence entre le nœud esclave et son projeté sur la maille maître;

4. Calculer la direction de descente \(\left\lbrace {d\zeta }_{i}\right\rbrace =-{\left[H\left({\zeta}_{i}\right)\right]}^{-1}.\left\lbrace R\left(\left\lbrace {\zeta}_{i}\right\rbrace \right)\right\rbrace\)

5. Calculer le point suivant tel que \(\left\lbrace {\zeta}_{i+1}\right\rbrace =\left\lbrace {\zeta}_{i}\right\rbrace +\left\lbrace {d\zeta }_{i}\right\rbrace\) ;

6. Si le processus a convergé, on s’arrête, sinon on boucle en 2.

Il n’y a pas d’algorithme de recherche linéaire contrairement au cas précédent.

Relations cinématiques

Définition de la matrice de contact

On effectue une modélisation idéalisée du phénomène de contact, en ce sens qu’il suppose les frontières des corps parfaitement définies par une ligne ou une surface : on écrit alors une condition de non interpénétration discrète et linéarisée [bib2] . Soit \(P\) un nœud esclave, \(M\) sa projection sur la maille maître qui a été déterminée lors de l’appariement. En 2D, cette maille maître a deux nœuds (SEG2) ou trois nœuds (SEG3). En 3D, elle peut en avoir trois, quatre, six, huit ou neuf (TRIA3, QUAD4, TRIA6, QUAD8, QUAD9). Le déplacement du point \(M\) est une combinaison linéaire des déplacements des nœuds de l’élément fini, avec pour coefficients les valeurs des fonctions de forme \(\Phi\) en \(M\) . Plaçons-nous dans le cas où la maille maître est un SEG2 pour simplifier l’exposé. On a alors:

(2243)\[\left\lbrace {u}_{M}\right\rbrace =\left[{\Phi}_{A}\left(M\right)\right].\left\lbrace {u}_{A}\right\rbrace +\left[{\Phi}_{B}\left(M\right)\right].\left\lbrace {u}_{B}\right\rbrace\]

Dans un premier temps, on choisitde prendre comme direction \(\left\lbrace N\right\rbrace\) la normale sortante de la maille maître (cf. ).

Figure 3.5.1-a : projection d’un nœud esclave sur une maille SEG2.

Le jeu normal s’écrit comme la différence entre le nœud esclave \(P\) et sa projection \(M\) sur la maille maitre:

(2244)\[\langle N\rangle .\left[\left\lbrace {u}_{P}\right\rbrace -\sum_{j=1}^{{n}_{m}}\left[{\Phi}_{{B}_{j}}(M)\right].\left\lbrace {u}_{{B}_{j}}\right\rbrace \right]=d\]

Si l’on écrit une telle relation pour tous les couples de contact, on obtient les conditions géométriques de non pénétration sous forme matricielle:

(2245)\[\left[{A}^{\text{c}}\right].\left\lbrace u\right\rbrace =d\]

La matrice \(\left[{A}^{\text{c}}\right]\) , appelée matrice de contact, contient une ligne par couple de contact, et autant de colonnes que de degrés de liberté physiques du problème. Supposons que l’on ait deux mailles de contact de type SEG2, selon le schéma de la .

Figure 3.5.1-b : écriture de la matrice de contact A sur un exemple.

Si l’on note par exemple \({u}_{B}\) le déplacement du nœud \(B\) selon la direction \(x\) , \({v}_{B}\) le déplacement du nœud \(B\) selon la direction \(y\) , et \({d}^{1}\) et \({d}^{2}\) les jeux courants pour les deux couples:

(2246)\[\begin{split}\left[{A}^{\text{c}}\right].\left\lbrace \begin{array}{c}{u}_{P}\\ {v}_{P}\\ {u}_{Q}\\ {v}_{Q}\\ {u}_{B}\\ {v}_{B}\\ {u}_{C}\\ {v}_{C}\\ {u}_{D}\\ {v}_{D}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{d}^{1}\\ {d}^{2}\end{array}\right\rbrace\end{split}\]

Avec la matrice de contact \(\left[{A}^{\text{c}}\right]\) :

(2247)\[\begin{split}\left[{A}^{\text{c}}\right]=\left[\begin{array}{cccccccccc}{N}_{x}^{1}& {N}_{y}^{1}& 0& 0& -\frac{1}{2}.{N}_{x}^{1}& -\frac{1}{2}.{N}_{y}^{1}& -\frac{1}{2}.{N}_{x}^{1}& -\frac{1}{2}.{N}_{y}^{1}& 0& 0\\ 0& 0& {N}_{x}^{2}& {N}_{y}^{2}& 0& 0& -\frac{3}{4}.{N}_{x}^{2}& -\frac{3}{4}.{N}_{y}^{2}& -\frac{1}{4}.{N}_{x}^{2}& -\frac{1}{4}.{N}_{y}^{2}\end{array}\right]\end{split}\]

On n’a ici considéré que les degrés de liberté des nœuds impliqués dans le contact; la matrice \(\left[{A}^{\text{c}}\right]\) devrait être plus creuse. Mais en pratique, on réduit toujours la matrice de contact sur les degrés de liberté actifs. On ne stocke donc que les coefficients non nuls.

Définition de la matrice de frottement

La notion de matrice de contact s’étend au cas des glissements tangentiels, sur le plan tangent. C’est la matrice des relations cinématiques de frottement \(\left[{A}^{\text{f}}\right]\) .

Choix de la normale

Dans le paragraphe précédent, on a choisi de prendre comme direction \(\left\lbrace N\right\rbrace\) la normale sortante de la maille maître. C’est le comportement par défaut dans code_aster . Toutefois, il est possible de choisir d’autres normales:

• La normale esclave NORMALE=”ESCL”;

• Une moyenne entre la normale maître et la normale esclave NORMALE=”MAIT_ESCL”;

Il est également possible de demander à utiliser des normales lissées (LISSAGE=”OUI”), c’est à dire qu’au lieu d’utiliser la normale maître au point de projection, on peut prendre une normale issue de l’interpolation linéaire entre les normales aux nœuds de la maille maître.

De même, le calcul des normales se fait toujours via les fonctions de forme de l’élément, c’est ce qu’on définit comme la «vraie» normale ou normale «automatique». Mais il est possible d’imposer une normale sur la maille maître, la maille esclave ou sur les deux de manière différente:

• Directement (VECT_MAIT ou VECT_ESCL =”FIXE”)

• Indirectement par utilisation d‘un trièdre (VECT_MAIT ou VECT_ESCL=”VECT_Y”). Dans ce dernier cas, la normale utilisée sera le vecteur issu du produit vectoriel entre la tangente à la maille et le vecteur VECT_Y donné.

Remarques:

• S’agissant d’un appariement de type nœud-facette, la normale esclave est calculée elle-même par lissage, l’option de LISSAGE n’a donc pas d’effet si on choisit NORMALE=”ESCL”;

• L’utilisation de normales prédéfinies est obligatoire dans le cas des poutres;

• Le choix d’une normale autre que la maille maître doit être limité à des cas exceptionnels comme lorsque le maillage ou les problèmes de compatibilités du contact avec les conditions aux limites ont imposé à l’utilisateur de prendre pour maître une surface maillée grossièrement;

• L’utilisation de TYPE_APPA=”FIXE” pour la recherche de la maille maître la plus proche (voir r5.03.50-calcul-projection-esclave-maitre) ne préjuge pas du choix de la normale dans l’écriture de la relation de non-pénétration, qui reste au choix de l’utilisateur. Mais il est plus cohérent de choisir une normale fixe ( VECT_MAIT=”FIXE”).

Coefficients de la matrice de contact

Éléments standards

Les valeurs des fonctions de forme \({\Phi}_{{B}_{j}}(M)\) des nœuds maîtres au point \(M\) pour les différentes mailles de contact sont standards (voir [r1.01.01]).

Mailles QUAD8

Les mailles de type quadrilatère à huit nœuds présentent un défaut. En effet, les fonctions de formes classiques ne sont pas positives sur tout le domaine et conduisent à des résultats aberrants lorsqu’elles sont utilisées dans le contact. Le symptôme principal lié à l’utilisation des fonctions de forme classiques est l’apparition de forces de contact négatives non-physiques qui provoquent des oscillations (entre les nœuds sommets et les nœuds milieux).

Pour éviter ce phénomène, code_aster procède à la modification de l’élément QUAD8en imposant des relations linéaires entre nœuds milieux et nœuds sommets et en utilisant finalement les fonctions de forme du QUAD4.

De manière générale, il est préférable d’éviter d’utiliser ce genre d’élément et de lui préférer des éléments quadratiques complets comme les QUAD9, car on considère les QUAD8comme des éléments linéaires et on introduit en plus des relations linéaires qui peuvent être gênantes.

Éléments de COQUE_3D

Les éléments de type COQUE_3D sont des éléments finis mixtes non isoparamétriques. Ils se basent sur des mailles quadratiques de type QUAD9 (respectivement TRIA7) mais le nœud milieu ne portant pas de degré de liberté de translation, on procède à une projection sur un QUAD8 (respectivement TRIA6) et on retombe donc sur les défauts évoqués au paragraphe précédent. Ces éléments finis sont donc déconseillés dans les problèmes de contact.

Introduction d’un jeu fictif

On peut vouloir modéliser le contact entre des structures a certaines particularités («trou» ou «bosse») que l’on ne souhaite pas mailler (voir figure ).

Figure 3.6-a : trous et bosses.

Une solution consiste à mailler la surface sans ces défauts et à y ajouter une distance donné par l’utilisateur (voir figure ).

Figure 3.6-b : surface maillée sans défaut.

On corrige la valeur du jeu:

(2248)\[\langle N\rangle .(\left\lbrace {u}_{P}\right\rbrace -\left\lbrace {u}_{M}\right\rbrace ).=d-({d}_{\text{e}}+{d}_{\text{m}})\]

où \({d}_{\text{e}}\) et \({d}_{\text{m}}\) sont donnés par l’utilisateur respectivement sous les mots-clés DIST_MAIT et DIST_ESCL pour chaque zone de contact. Ces distances sont signées: elles représentent la translation à appliquer au nœud du maillage dans la direction de la normale \(\left\lbrace n\right\rbrace\) sortante pour obtenir le point de la structure réelle. Ces mots-clés permettent de rendre compte également du contact entre coques dont seules les surfaces moyennes sont maillées: \({d}_{\text{e}}\) et \({d}_{\text{m}}\) valent alors la demi-épaisseur des coques (valeurs positives).

Remarques:

• Si l’on utilise DIST_MAIT et DIST_ESCL, il faut prendre garde à l’interprétation visuelle des résultats. Si \(({d}_{\text{e}}+{d}_{\text{m}})>0\) , le code pourra annoncer du contact alors que la visualisation montrera un espacement des deux maillages. Si \(({d}_{\text{e}}+{d}_{\text{m}})<0\) , le code pourra annoncer du contact alors que la visualisation montrera deux maillages interpénétrés;

• Pour se souvenir des signes, penser à:

• \(({d}_{\text{e}}+{d}_{\text{m}})>0\) : ajout de matière par rapport au maillage

• \(({d}_{\text{e}}+{d}_{\text{m}})<0\) : ablation de matière par rapport au maillage

Les options DIST_POUTRE et DIST_COQUE s’appuient sur les caractéristiques élémentaires définies dans l’opérateur AFFE_CARA_ELEM pour ajouter le jeu fictif correspondant à l’épaisseur (dans le cas des coques) ou au rayon (dans le cas des poutres à section circulaire).

Réactualisation géométrique

Dans le cadre de la modélisation du contact en grands déplacements, l’évolution de la géométrie des surfaces joue un rôle fondamental. En effet, c’est elle qui conditionne le calcul des normales aux faces potentiellement en contact et donc qui conditionne l’appariement réalisé.

La réactualisation géométrique est définie par le mot-clé REAC_GEOM de la commande DEFI_CONTACT. Son fonctionnement est le suivant:

• Si REAC_GEOM=”SANS”: il n’y a pas de réactualisation géométrique. Tout le calcul s’effectue sur la configuration initiale avec l’appariement initial.

• Si REAC_GEOM=”CONTROLE”: il faut renseigner NB_ITER_GEOM. Au sein du même pas de charge, on effectue NB_ITER_GEOM fois le cycle itération jusqu’à convergence, réactualisation géométrique, appariement .

• Si REAC_GEOM=”AUTOMATIQUE”: la décision de refaire un appariement géométrique est faite automatiquement par le logiciel. Le critère est le suivant:

(2249)\[\frac{{\mathrm{\parallel }\Delta {u}^{+}-\Delta {u}^{-}\mathrm{\parallel }}_{\infty}}{{\mathrm{\parallel }\Delta {u}^{-}\mathrm{\parallel }}_{\infty}}<\text{RESI\_GEOM}\]

RESI_GEOM vaut \(1\text{\%}\) par défaut. Si la norme infinie du déplacement entre deux instants, divisée par la norme infinie du déplacement obtenu depuis l’équilibre est supérieur à RESI_GEOM, alors on réactualise.La norme infinie est définie par:

(2250)\[{\parallel \Delta u\parallel }_{\infty}=\max\sqrt{{u}_{x}^{2}+{u}_{y}^{2}+{u}_{z}^{2}}\]

Remarques:

• A la première itération ou en cas de mouvements de corps rigide (en dynamique), \({\mathrm{\parallel }\Delta {u}^{-}\mathrm{\parallel }}_{\infty}=0\) . Pour éviter de diviser par zéro, on force la réactualisation;

• La norme infinie est évaluée sur tous les nœuds du maillage (et non pas seulement sur les nœuds en contact) portant des degrés de liberté de déplacement.

On peut tout d’abord remarquer que l’appariement est soumis à la phase de réactualisation géométrique. En outre, le fait de réaliser plusieurs fois au sein du même pas de charge le cycle itération jusqu’à convergence, réactualisation géométrique, appariement permet de suivre l’évolution de la géométrie de la structure. Il faut en effet souligner que cette évolution géométrique est une des composantes non linéaires d’un calcul de contact en grands déplacements.

Dans la pratique, on peut conseiller les choix suivants pour le mot clé REAC_GEOM:

• Pour un calcul en petits déplacements, REAC_GEOM=”SANS”. On travaille sur la configuration initiale.

• Pour des calculs en grands déplacements, soit utiliser REAC_GEOM=”AUTOMATIQUE”(valeur par défaut), soit utiliser REAC_GEOM=”CONTROLE’et une valeur pour NB_ITER_GEOM dépendant de l’importance de l’évolution géométrique des surfaces.

Dans le cas où l’utilisateur ne laisse pas à code_aster la possibilité de gérer automatiquement la réactualisation géométrique, le code l’avertira par une alarme si le critère automatique (valant 1%)n’est pas assuré du fait du choix de l’utilisateur.

Problème mécanique de contact/frottement

La prise en compte du contact-frottement dans un problème mécanique a deux conséquences:

• La modification de l’équation d’équilibre pour prendre en compte les réactions de contact-frottement aux interfaces;

• L’application de lois supplémentaires régissant le contact et le frottement (voir r5.03.50-mecanique-contact-frottement) pour calculer ces réactions mais aussi imposer des conditions sur la cinématique.

Pour le contact sans frottement, les conditions de Signorini se ramènent à un problème d’optimisation sous contrainte classique (Kuhn & Tucker). En revanche, pour le frottement de Coulomb, on ne peut pas écrire de problème d’optimisation équivalent sans faire d’hypothèse(s) supplémentaire(s) .

Dans code_aster , les méthodes discrètes de résolution du problème de contact/frottement sont fondées sur une approche découplée entre l’équilibre et la loi de contact/frottement. Après la résolution du problème de mécanique sans contact-frottement, on corrige la solution (cinématique et réactions) en appliquant la loi de Signorini-Coulomb. Cette stratégie permet de ne faire aucune autre hypothèse sur la nature du problème mécanique, que ce soit la cinématique ou les relations de comportement. Toutefois, il est indispensable que la matrice de rigidité du problème sans contact soit symétrique .

Remarques:

• L’expression «faire un calcul avec contact» veut dire que l’on écrit les relations de non-pénétration, mais n’implique pas qu’il y ait contact effectif pour le chargement considéré. Or c’est la résolution du problème de contact-frottement qui coûte le plus cher.

• Le contact agissant comme correction sur les résultats issus d’un calcul mécanique classique, il est indispensable que le problème sans contact soit mécaniquement bien posé et numériquement soluble. En particulier, les éventuels mouvements de corps rigides doivent être supposés éliminés hors résolution du problème de contact.

Problème mécanique sans contact/frottement

La résolution d’un problème non-linéaire dans l’opérateur STAT_NON_LINE (ou DYNA_NON_LINE) est décrite en détail dans le document [r5.03.01]. À chaque pas de temps \(i\) , on cherche à vérifier l’équilibre global de la structure:

(2251)\[\begin{split}\lbrace \begin{array}{c}\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace +{\left[B\right]}^{T}.\left\lbrace {\lambda}_{i}\right\rbrace =\left\lbrace {L}_{i}^{\text{méca}}({u}_{i})\right\rbrace \\ \left[B\right].\left\lbrace {u}_{i}\right\rbrace =\left\lbrace {u}_{i}^{d}\right\rbrace \end{array}\end{split}\]

Afin d’éviter de surcharger les équations, on fait l’hypothèse queles conditions limites sont éliminées et où l’on n’a donc pas de matrice \(\left[B\right]\) , on cherche donc à résoudre:

(2252)\[\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{méca}}({u}_{i})\right\rbrace\]

Ce problème non-linéaire est résolu par une méthode itérative de type Newton-Raphson qui a pour caractéristiques:

• Un découpage a priori du chargement en «pas de temps» (notés par l’indice \(i\) );

• Une linéarisation du problème d’équilibre par la méthode de Newton (les itérations étant notées par l’indice \(n\) ).

Les inconnues sont calculées de façon incrémentale. À partir de \(\left\lbrace {u}_{i-1}\right\rbrace\) , solution satisfaisant l’équilibre en \({t}_{i-1}\) , on détermine \(\left\lbrace \Delta {u}_{i}\right\rbrace\) qui permet d’obtenir la solution en \({t}_{i}\) :

(2253)\[\left\lbrace {u}_{i}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}\right\rbrace\]

L’incrément \(\left\lbrace \Delta {u}_{i}\right\rbrace\) estd’abord estimé en linéarisant le problème par rapport au temps autour de \((\left\lbrace {u}_{i-1}\right\rbrace ,{t}_{i-1})\) (phase de prédiction ou d’Euler). Puis on utilise une méthode de Newton ou une de ses variantes pour résoudre l’équation (59) de manière itérative: on calcule une suite \(\left\lbrace \delta {u}_{i}^{n}\right\rbrace\) où l’exposant \(n\) est le numéro de l’itération. Pour simplifier, on ne ferapas de distinction entrela phase de prédiction etla phase de correction de Newton. On écrit finalement:

(2254)\[\left\lbrace {u}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace\]

On se place au temps \({t}_{i}\) et à l’itération de Newton \(n\) . On a utilisé les notations suivantes:

• \(\left\lbrace {u}_{i}^{n}\right\rbrace\) : déplacementsà l’instant \({t}_{i}\) et à l’itération de Newton \(n\) ;

• \(\left\lbrace {u}_{i-1}\right\rbrace\) : déplacementà l’instant \({t}_{i-1}\) . Cette solution respecte la condition d’équilibre de la structure;

• \(\left\lbrace \delta {u}^{n}\right\rbrace\) : incrément des déplacementspourl’itération de Newton \(n\) ;

• \(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace\) : incrément desdéplacementscumulésdepuis le début du pas de temps, avant l’itération de Newton \(n\) ;

Après linéarisation de (59), pour l’itération de Newton \(n\) , on introduit [1]_ la matrice tangente mécanique \(\left[{K}^{\text{m},n-1}\right]\) eton doit alors trouver \(\left\lbrace \delta {u}^{n}\right\rbrace\) tel que:

(2255)\[\left[{K}^{\text{m},n-1}\right].\left\lbrace \delta {u}^{n}\right\rbrace =\left\lbrace {L}_{i}^{\text{méca},n-1}\right\rbrace -\left\lbrace {L}_{i}^{int,n-1}\right\rbrace\]

Modification de l’équation d’équilibre avec contact/frottement

Pour assurer l’équilibre de la structure, il faut introduire des forces de contact/frottement. De la manière la plus générale possible, on écrit que le chargement extérieur mécanique est modifié par les efforts de contact/frottement:

(2256)\[\left\lbrace {L}_{i}^{\text{méca}}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{c}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{f}}({u}_{i})\right\rbrace\]

Où:

• \(\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace\) est le vecteur des forces extérieures (conditions de Neumann);

• \(\left\lbrace {L}_{i}^{\text{c}}({u}_{i})\right\rbrace\) est le vecteur des forces de contact;

• \(\left\lbrace {L}_{i}^{\text{f}}({u}_{i})\right\rbrace\) est le vecteur des forces de frottement;

A priori, toutes ces quantités sont non-linéaires car elles dépendent du vecteur déplacement \(\left\lbrace {u}_{i}\right\rbrace\) (on parle de chargements «suiveurs»). En injectant (63) dans (59), on cherche donc à résoudre:

(2257)\[\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{c}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{f}}({u}_{i})\right\rbrace\]

A chaque itération de Newton \(n\) , on linéarise l’équation (64) par rapport à \(\left\lbrace {u}_{i}^{n}\right\rbrace\) .Ce processus introduit la matrice tangente \(\left[{K}^{\text{m},n-1}\right]\) qui contiendra les contributions issues de la linéarisation des efforts intérieurs et extérieurs et la matrice \(\left[{K}^{\text{cf},n-1}\right]\) qui contiendra les contributions issues de la linéarisation des forces de contact/frottement, on doit alors trouver \(\left\lbrace \delta {u}^{n}\right\rbrace\) tel que:

(2258)\[(\left[{K}^{\text{m},n-1}\right]+\left[{K}^{\text{cf},n-1}\right]).\left\lbrace \delta {u}^{n}\right\rbrace =\left\lbrace {L}_{i}^{\text{ext},n-1}\right\rbrace -\left\lbrace {L}_{i}^{int,n-1}\right\rbrace -\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace -\left\lbrace {L}_{i}^{\text{f},n-1}\right\rbrace\]

Lois de contact/frottement

Les lois de Signorini et de Coulomb impliquent des inégalités et des égalités. La formulation discrète disponible dans code_aster consiste à modifier les relations d’inégalité en relations d’égalité. Pour ce faire, il y a deux méthodes possibles:

1. Dualiser les lois de contact/ frottement;

2. Régulariser les lois de contact/ frottement .

Cinématique

Matrices de contact/frottement

Dans le paragraphe r5.03.50-definition-matrice-contact , nous avons vu comment écrire les conditions cinématiques de contact en introduisant la matrice de contact \(\left[{A}^{\text{c}}\right]\) . De la même manière, mais en considérant les conditions cinématiques de frottement (sur le plan tangent), on évaluera les matrices de glissement \(\left[{A}^{\text{g}}\right]\) et d’adhérence \(\left[{A}^{\text{a}}\right]\) . Ces matrices quasi-pleines et rectangulaires sont calculées par l’évaluation des relations cinématiques entre le nœud esclave et la maille maître appariée. La construction des relations cinématiques se fait lors de la procédure d’appariement géométrique qui intervient au début de chaque pas de temps et à chaque réactualisation géométrique. Pour des raisons de performances, on n’utilise et on ne construit que la sous-partie des ces matrices correspondant à l’activation des différents seuils (contact ou glissement), au sein des algorithmes de résolution du contact/frottement.

Ces matrices dépendent de la réactualisation géométrique.

Au niveau global, les matrices ne dépendant pas de la solution sans contact du problème mécanique.

Conditions unilatérales de contact

La relation de non pénétration consiste à dire que déplacement relatif selon une direction donnée ne peut pas dépasser le jeu initial \(\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace\) , mesuré sur le maillage,dans cette direction. La condition unilatérale de contact s’écrit (voir r5.03.50-definition-matrice-contact ):

(2259)\[\left[{A}^{\text{c}}\right].\left\lbrace {u}_{i}\right\rbrace \le \left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace\]

Cette équation traduit le fait que tout mouvement de la structure doit se faire avec le respect de la condition de non-pénétration, ou encore que le déplacement des nœuds de la surface de contact est inférieur au jeu initial \(\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace\) . Après résolution du problème mécanique par la méthode de Newton, la condition de non-pénétration devient:

(2260)\[\left[{A}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )\le \left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace\]

Cette condition peut s’écrire de manière itérative :

(2261)\[\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace \le \left\lbrace {d}^{\text{c},n-1}\right\rbrace\]

Avec \(\left\lbrace {d}^{\text{c},n-1}\right\rbrace\) le jeu évalué avant l’itération de Newton courante \(n\) :

(2262)\[\left\lbrace {d}^{\text{c},n-1}\right\rbrace =\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace )\]

Conditions d’adhérence

En adhérence, les nœuds ne bougent pas sur le pas de temps, c’est-à-dire:

(2263)\[\left[{A}^{\text{a}}\right].(\left\lbrace {u}_{i}\right\rbrace -\left\lbrace {u}_{i-1}\right\rbrace )=0\]

\(\left[{A}^{\text{a}}\right]\) est la matrice des nœuds en contact adhérent, c’est-à-dire la sous-partie de la matrice de frottement \(\left[{A}^{\text{f}}\right]\) (voir r5.03.50-definition-matrice-contact) appliquée aux nœuds en contact adhérent.

Conditions de glissement

La contrainte de cisaillement \(\left\lbrace {r}_{i}\right\rbrace\) est colinéaire à la direction tangente de glissement, soit:

(2264)\[\left[{A}^{\text{g}}\right].(\left\lbrace {u}_{i}\right\rbrace -\left\lbrace {u}_{i-1}\right\rbrace )=\lambda .\left\lbrace {r}_{i}\right\rbrace\]

\(\left[{A}^{\text{g}}\right]\) est la matrice des nœuds en contact glissant , c’est-à-dire la sous-partie de la matrice de frottement \(\left[{A}^{\text{f}}\right]\) (voir r5.03.50-definition-matrice-contact) appliquée aux nœuds en contact glissant . On a:

(2265)\[\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )=\lambda .\left\lbrace {r}_{i}(\left\lbrace {u}_{i}\right\rbrace )\right\rbrace\]

L’équation ( ) dépend du déplacement final \(\left\lbrace {u}_{i}\right\rbrace\) . Les conditions de glissement sont introduites avec un multiplicateur de Lagrange \(\left\lbrace {\mu}^{\text{g}}\right\rbrace\) tel que:

(2266)\[\left\lbrace {\mu}^{\text{g}}({u}_{i})\right\rbrace =\mu .\left[∣{\mu}^{\text{c}}({u}_{i})∣\right].\left\lbrace {t}^{n}\right\rbrace\]

Avec \({t}^{n}\) le vecteur unitaire de la direction de glissement qui vaut en 3D:

(2267)\[\left\lbrace t\right\rbrace =\frac{\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )}{∥\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )∥}\]

Soit encore:

(2268)\[ \begin{align}\begin{aligned}\left\lbrace {\mu}^{\text{g}}({u}_{i})\right\rbrace =\left[{k}^{\text{g}}({u}_{i})\right].\left\lbrace {t}^{n}\right\rbrace\\\textrm{avec}\\\left[{k}^{\text{g}}({u}_{i})\right]=\mu .\left[∣{\mu}^{\text{c}}({u}_{i})∣\right]\end{aligned}\end{align} \]

Le multiplicateur de Lagrange est non linéaire et dépend de la solution \(\left\lbrace {u}_{i}\right\rbrace\) .En 2D, c’est aussinon-linéaire mais le multiplicateur ne dépend pas de la direction de glissement, on a donc:

(2269)\[{\left\lbrace {\mu}^{\text{g}}({u}_{i})\right\rbrace }_{\text{2D}}=\mu .\left[∣{\mu}^{\text{c}}({u}_{i})∣\right]=\left[{k}^{\text{g}}({u}_{i})\right]\]

Dualisation

Dans le cas de la dualisation, on utilise un algorithme essai-erreur qui postule a priori l’état d’une liaison et qui vérifie son état après application de la loi de Signorini-Coulomb. Pour prendre en compte les contraintes (on parle des contraintes dans le sens de conditions de «restriction» et non dans le sens mécanique du terme) portant sur le champ de déplacements ou sur les efforts, on les fait intervenir dans les équations au travers de multiplicateurs de Lagrange (comme cela peut être fait dans code_aster pour les conditions aux limites cinématiques). On introduit deux ensembles de multiplicateurs de Lagrange:

• \(\left\lbrace {\mu}^{\text{c}}\right\rbrace\) portant sur les conditions de contact;

• \(\left\lbrace {\mu}^{\text{g}}\right\rbrace\) portant sur les conditions de glissement.

En écrivant l’équilibre de la structure, on donne l’interprétation suivante des multiplicateurs de Lagrange:

• \(\left\lbrace {L}_{i}^{\text{c}}({u}_{i})\right\rbrace =[{A}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) représente les forces nodales de contact;

• \(\left\lbrace {L}_{i}^{\text{g}}({u}_{i})\right\rbrace =[{A}^{\text{g}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{g}}\right\rbrace\) représente les forces nodales de glissement.

Même si l’algorithme de frottement utilise la régularisation, on verra que l’introduction des «multiplicateurs de Lagrange» est nécessaire pour linéariser les forces nodales de glissement.

Régularisation

Le principe de la régularisation est de modifier les lois de contact-frottement pour en obtenir de plus faciles à manipuler, par exemple, que les relations deviennent univoques et dérivables (pour pouvoir appliquer la méthode de Newton par exemple). On introduit donc des hypothèses qui rendent le modèle non-exact par rapport aux lois de Signorini-Coulomb.

Régularisation des conditions de contact

Le principe de la régularisation est de modifier le graphe de la loi de contact afin de supprimer le caractère non-univoque de la relation de contact.

Figure 4.3.3.1-a : régularisation de la condition de contact (forme continue) .

On rappelle que la condition de non-interpénétration s’écrit sous forme itérative:

(2270)\[\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace \le \left\lbrace {d}^{\text{c},n-1}\right\rbrace\]

L’idée de la régularisation est de pénaliser la situation pour laquelle il y a interpénétration. L’interpénétration \(\left\lbrace {h}^{n}\right\rbrace\) vaut:

(2271)\[\left\lbrace {h}^{n}\right\rbrace =\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace >0\]

En d’autres termes, on écrit que la force de contact est d’autant plus grande que l’interpénétration \(\left\lbrace {h}^{n}\right\rbrace\) est importante, d’où la forme régularisée de la force de contact:

(2272)\[ \begin{align}\begin{aligned}{\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}={E}_{N}.[{A}^{\text{c}}{]}^{T}.{\left\lbrace {h}^{n}\right\rbrace }^{+}\\\textrm{avec}\\{E}_{N}>0\end{aligned}\end{align} \]

\({E}_{N}\) est le coefficient de régularisation (ou de pénalisation, de dimension d’un rapport d’une force par un déplacement)du contact et on note \({\square}^{+}\) la partie positive d’une quantité. Lecoefficient de pénalisation de contact s’interprète comme le ressort de raideur \({E}_{N}\) qui s’oppose à la pénétration du nœud esclave dans la surface maître en lui appliquant une force nodale correspondante. Plus il est grand, moins il y a d’interpénétration et plus la force de rappel est élevée.

Régularisation de la condition d’adhérence

Le principe de la régularisation de l’adhérenceest de modifier le graphe de la loi de Coulomb. En premier lieu,onsupprime le caractère non-univoque de la loi dans la partie adhérente.

Figure 4.3.3.2-a : régularisation de la condition d’adhérence – Modèle de Tresca (forme continue) .

On établit la forme régularisée de la force d’adhérence:

(2273)\[ \begin{align}\begin{aligned}{\left\lbrace {L}_{i}^{\text{a},n}\right\rbrace }_{\text{regu}}={E}_{T}.{\left[{A}^{\text{a}}\right]}^{T}.\left[{A}^{\text{a}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )\\\textrm{avec}\\{E}_{T}>0\end{aligned}\end{align} \]

\({E}_{T}\) est le coefficient de régularisation (ou de pénalisation, de dimension d’un rapport d’une force par un déplacement)du frottement. La notion d’adhérence à proprement parler va donc disparaître, tous les nœuds glissent. On voit sur la figure ci-dessus que cette régularisation n’est pas satisfaisante car on n’a fait que transformer le problème de Coulomb en problème de Tresca: la force est proportionnelle au glissement relatif. Pour s’approcher du modèle de Coulomb, il faut ajouter une inéquation supplémentaire:

(2274)\[\Vert {E}_{T}.\left[{A}^{\text{a}}\right].\left(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace \right)\Vert \le \mu .\left|p\left(u\right)\right|\]

Cette inéquation supplémentaire modifie le graphe (voir figure ( )). Il en résulte une forme approchée du modèle de Coulomb, dans laquelle on définira les nœuds adhérents par rapport au seuil de contact, produit du coefficient de Coulomb \(\mu\) et de la pression de contact \(p\left(u\right)\) . Le caractère approché de la loi de Coulomb se traduit par le fait que les nœuds dits «adhérents» glisseront d’autant plus que le coefficient de régularisation sera faible. L’interprétation physique est donc moins directe que dans le cas du contact, ce qui explique en grande partie la difficulté à trouver en pratique une valeur satisfaisante de ce coefficient. Toutefois l’erreur commise («fausse» détection du seuil de glissement) est souvent peu importante par rapport aux hypothèses de modélisation. Il reste une non-linéarité qui concerne la partie «glissante» de la loi. Le traitement de cette non-linéarité est reporté au paragraphe r5.03.50-modelisation-glissement.

Figure 4.3.3.2-b : régularisation de la condition d’adhérence – Modèle de Coulomb (forme continue) .

Modélisation du glissement

Les conditions de glissement sont toujours introduites avec un multiplicateur de Lagrange \(\left\lbrace {\mu}^{\text{g}}\right\rbrace\) quis’exprime (en 3D) à partir de (75):

(2275)\[\left\lbrace {\mu}^{\text{g}}({u}_{i})\right\rbrace =\left[{k}^{\text{g}}({u}_{i})\right].\left\lbrace {t}^{n}\right\rbrace\]

Ce qui fait apparaître les deux inconnues qui dépendent toutes les deux de la solution: la pression de contact dans la matrice des seuils de glissement \(\left[{k}^{\text{g}}\right]\) et la direction de glissement \(\left\lbrace {t}^{\text{n}}\right\rbrace\) . Dans la pratique, cependant, on considère que la connaissance du seuil de glissement \(\left[{k}^{\text{g},n-1}\right]\) est acquise à l’itération précédente \((n-1)\) , ce qui revient à se ramener à un critère de Tresca pour chaque itération. A convergence le seuil est évidemment fixe: il n’y a donc plus de différences entre les seuils au cours des itérations. \(\left\lbrace {\mu}^{\text{g},n}\right\rbrace\) est donc approximé par rapport à l’état de l’itération précédente \((n-1)\) :

(2276)\[\left\lbrace {\mu}^{\text{g},n}\right\rbrace \approx \left[{k}^{\text{g},n-1}\right].\frac{\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )}{∥\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )∥}=\left[{k}^{\text{g},n-1}\right].\left\lbrace {t}^{n}\right\rbrace\]

On a également le seuil de glissement à l’itération précédente:

(2277)\[\left[{k}^{\text{g},n-1}\right]=\mu .\left[∣{\mu}^{\text{c},n-1}∣\right]\]

Linéarisations

Pour résoudre le problème non-linéaire, il est nécessaire de linéariser les quantités qui dépendent de \(\left\lbrace {u}_{i}^{n}\right\rbrace\) dans l’équation d’équilibre (64). L’écriture des forces de contact/frottement introduit des quantités non-linéaires qu’il convient donc de linéariser pour pouvoir appliquer l’algorithme de Newton. Dans le tableau suivant, on recense les cas où la linéarisation est nécessaire suivant le type d’algorithme.

	Dualisation	Régularisation
\(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace\)	Pas de linéarisation	Linéarisation
\(\left\lbrace {L}_{i}^{\text{a},n}\right\rbrace\)		Linéarisation
\(\left\lbrace {L}_{i}^{\text{g},n}\right\rbrace\)		Linéarisation

Linéarisation des forces de contact

On régularise la condition de contact à partir de l’expression (79) de la force de contact, quand on est dans la phase d’interpénétration, c’est-à-dire en supposant \(\left\lbrace {h}^{n}\right\rbrace >0\) :

(2278)\[{\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}={E}_{N}.[{A}^{\text{c}}{]}^{T}.(\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace )\]

En linéarisant, on obtient:

(2279)\[{\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}\underset{\text{linéarisation}}{\to }{\left\lbrace {\stackrel{ˆ}{L}}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}=[{K}_{i}^{\text{c},n-1}].\left\lbrace \delta {u}^{n}\right\rbrace +\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace\]

\([{K}_{i}^{\text{c},n-1}]\) apporte une nouvelle contribution à la matrice tangente du problème, c’est la matrice tangente «de contact». Elle vaut évidemment:

(2280)\[[{K}_{i}^{\text{c},n-1}]={E}_{N}.[{A}^{\text{c},n-1}{]}^{T}.[{A}^{\text{c},n-1}]\]

Et \(\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace\) contribue au second membre, il vaut:

(2281)\[\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace =-{E}_{N}.[{A}^{\text{c},n-1}{]}^{T}.\left\lbrace {d}^{\text{c},n-1}\right\rbrace\]

On rappelle que le coefficient de pénalisation de contact s’interprète comme le ressort de raideur \({E}_{N}\) qui s’oppose à la pénétration du nœud esclave dans la surface maître. Plus il est grand, moins il y a d’interpénétration et plus la force de rappel est élevée, mais ce coefficient intervient aussi dans la matrice tangente, ce qui modifie son conditionnement et rend la résolution du système linéaire plus difficile. Quand le contact n’est pas activé (c’est-à-dire si \(\left\lbrace {h}^{n}\right\rbrace \le 0\) ), la matrice et le vecteur second membre sont nuls.

Linéarisation des forces d’adhérence

On régularise la condition d’adhérence à partir de l’expression (4442) de la force d’adhérence, quand on est dans la phase d’adhérence (dans le sens défini dans le paragraphe r5.03.50-regularisation-condition-adherence), c’est-à-dire en supposant que l’inégalité (4443) est strictement respectée:

(2282)\[{\left\lbrace {L}_{i}^{\text{a},n}\right\rbrace }_{\text{regu}}={E}_{T}.{\left[{A}^{\text{a}}\right]}^{T}.\left[{A}^{\text{a}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )\]

En linéarisant, on obtient:

(2283)\[{\left\lbrace {L}_{i}^{\text{a},n}\right\rbrace }_{\text{regu}}\underset{\text{linéarisation}}{\to }{\left\lbrace {\stackrel{ˆ}{L}}_{i}^{\text{a},n}\right\rbrace }_{\text{regu}}=[{K}_{i}^{\text{a},n-1}].\left\lbrace \delta {u}^{n}\right\rbrace +\left\lbrace {L}_{i}^{\text{a},n-1}\right\rbrace\]

\([{K}_{i}^{\text{a},n-1}]\) apporte une nouvelle contribution à la matrice tangente du problème, c’est la matrice tangente «d’adhérence». Elle vaut évidemment:

(2284)\[[{K}_{i}^{\text{a},n-1}]={E}_{T.}[{A}^{\text{a},n-1}{]}^{T.}[{A}^{\text{a},n-1}]\]

Et \(\left\lbrace {L}_{i}^{\text{a},n-1}\right\rbrace\) contribue au second membre, il vaut:

(2285)\[\left\lbrace {L}_{i}^{\text{a},n-1}\right\rbrace ={E}_{T.}[{A}_{\text{a}}^{n-1}{]}^{T.}\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace\]

On rappelle que lecoefficient de pénalisation de contact s’interprète comme le ressort de raideur \({E}_{T}\) (de dimension d’un rapport d’une force par un déplacement) qui s’oppose au glissement relatifdes deux surfaces. Plus il est grand,plus la force de rappel est élevée, mais ce coefficient intervient aussi dans la matrice tangente, ce qui modifie son conditionnement et rend la résolution du système linéaire plus difficile.

Linéarisation des forces de glissement

Nous avons vu dans le paragraphe r5.03.50-modelisation-glissement que les forces de glissement dépendaient du déplacement final \(\left\lbrace {u}_{i}\right\rbrace\) par le seuil de Tresca, mais que nous faisions l’hypothèse que les informations à l’itération \((n-1)\) suffisaient pour résoudre la loi de Coulomb [2]_ en écrivant:

(2286)\[\left\lbrace {\mu}^{\text{g},n}\right\rbrace =\left[{k}^{\text{g},n-1}\right].\left\lbrace {t}^{n}\right\rbrace\]

Néanmoins, il faut préparer ce qui se passera à l’itération suivante \((n+1)\) si le seuil n’est pas correct et nécessitera une nouvelle itération. A la prochaine itération de Newton, le déplacement solution sera:

(2287)\[\left\lbrace {u}_{i}^{n+1}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n}\right\rbrace +\left\lbrace \delta {u}^{n+1}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace +\left\lbrace \delta {u}^{n+1}\right\rbrace\]

Il faut donc linéariser la quantité (4455) en \((n+1)\) :

(2288)\[\left\lbrace {\mu}^{\text{g},n+1}\right\rbrace =\left[{k}^{\text{g},n}\right].\left\lbrace t\right\rbrace\]

L’équation (4457) est de la forme \(\frac{\left\lbrace h({x}^{n+1})\right\rbrace }{∥\left\lbrace h({x}^{n+1})\right\rbrace ∥}\) avec \(\left\lbrace {x}^{n+1}\right\rbrace =\left\lbrace {x}^{n}\right\rbrace +\left\lbrace {\delta x}^{n+1}\right\rbrace\) . La linéarisation s’écrit:

(2289)\[\left[k\right].\frac{\left\lbrace h({x}^{n+1})\right\rbrace }{∥\left\lbrace h({x}^{n+1})\right\rbrace ∥}\underset{\text{linéarisation}}{\to }\left[k\right].(\frac{\left\lbrace h({x}^{n})\right\rbrace }{∥\left\lbrace h({x}^{n})\right\rbrace ∥}+\frac{1}{∥\left\lbrace h({x}^{n})\right\rbrace ∥}.(\left[I\right]-\frac{\left\lbrace h({x}^{n})\right\rbrace \langle \left\lbrace h({x}^{n})\right\rbrace \rangle }{{∥\left\lbrace h({x}^{n})\right\rbrace ∥}^{2}}).\left\lbrace \delta {x}^{n+1}\right\rbrace )\]

On applique (4458) sur l’expression (4457) avec:

(2290)\[ \begin{align}\begin{aligned}\left\lbrace {x}^{n}\right\rbrace \leftarrow \left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace\\\textrm{et}\\\left\lbrace {x}^{n+1}\right\rbrace \leftarrow \left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace +\left\lbrace \delta {u}^{n+1}\right\rbrace\\\left\lbrace {\delta x}^{n+1}\right\rbrace \leftarrow \left\lbrace \delta {u}^{n+1}\right\rbrace\\\left\lbrace h({x}^{n})\right\rbrace \leftarrow \left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )\\\textrm{et}\\\left\lbrace h({x}^{n+1})\right\rbrace \leftarrow \left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace +\left\lbrace \delta {u}^{n+1}\right\rbrace )\\\left[k\right]\leftarrow \left[{k}^{\text{g},n}\right]\end{aligned}\end{align} \]

On linéarise (4457) :

(2291)\[\left\lbrace {\mu}^{\text{g},n+1}\right\rbrace \underset{\text{linéarisation}}{\to }\left\lbrace {\stackrel{ˆ}{\mu}}^{\text{g},n+1}\right\rbrace\]

On commence par noter \(\left\lbrace {g}^{\text{t},n}\right\rbrace\) le glissement tangentiel relatif à l’itération \(n\) :

(2292)\[\left\lbrace {g}^{\text{t},n}\right\rbrace =\left[{A}^{\text{g}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace )\]

On obtient:

(2293)\[\left\lbrace {\stackrel{ˆ}{\mu}}^{\text{g},n+1}\right\rbrace =\left[{k}^{\text{g},n}\right].(\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace }{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}+\frac{1}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}.(\left[I\right]-\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace \langle {g}^{\text{t},n}\rangle }{{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}^{2}}).\left\lbrace \delta {u}^{n+1}\right\rbrace )\]

Si on note:

(2294)\[\left\lbrace {\rho}^{\text{g},n}\right\rbrace =\left[{k}^{\text{g},n}\right].\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace }{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}\]

Et:

(2295)\[\left[{B}^{\text{g},n}\right]=\frac{\left[{k}^{\text{g},n}\right]}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}.(\left[I\right]-\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace \langle {g}^{\text{t},n}\rangle }{{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}^{2}})\]

On peut écrire (4462) sous la forme condensée suivante:

(2296)\[\left\lbrace {\stackrel{ˆ}{\mu}}^{\text{g},n+1}\right\rbrace =\left[{B}^{\text{g},n}\right].\left\lbrace \delta {u}^{n+1}\right\rbrace +\left\lbrace {\rho}^{\text{g},n}\right\rbrace\]

La force de glissement s’exprime en fonction du multiplicateur de Lagrange:

(2297)\[\left\lbrace {L}_{i}^{\text{g},n+1}\right\rbrace ={\left[{A}^{\text{g}}\right]}^{T}.\left\lbrace {\mu}^{\text{g},n+1}\right\rbrace\]

La matrice \(\left[{A}^{\text{g}}\right]\) ne dépendant du déplacement, la force de glissement linéarisée peut s’écrire:

(2298)\[\left\lbrace {L}_{i}^{\text{g},n+1}\right\rbrace \underset{\text{linéarisation}}{\to }\left\lbrace {\stackrel{ˆ}{L}}_{i}^{\text{g},n+1}\right\rbrace ={\left[{A}^{\text{g}}\right]}^{T}.\left\lbrace {\stackrel{ˆ}{\mu}}^{\text{g},n+1}\right\rbrace =\left[{K}^{\text{g},n}\right].\left\lbrace \delta {u}^{n+1}\right\rbrace +\left\lbrace {L}^{\text{g},n}\right\rbrace\]

Avec la matrice de glissement \(\left[{K}^{\text{g},n}\right]\) telle que:

(2299)\[\left[{K}^{\text{g},n}\right]={\left[{A}^{\text{g}}\right]}^{T}.\left[{B}^{\text{g},n}\right]\]

Et le vecteur:

(2300)\[\left\lbrace {L}^{\text{g},n}\right\rbrace ={\left[{A}^{\text{g}}\right]}^{T}.\left\lbrace {\rho}^{\text{g},n}\right\rbrace\]

\(\left[{K}^{\text{g},n}\right]\) apporte une nouvelle contribution à la matrice tangente du problème, c’est la matrice tangente «de glissement». Elle vaut évidemment:

(2301)\[\left[{K}^{\text{g},n}\right]=\frac{{\left[{A}^{\text{g}}\right]}^{T}.\left[{k}^{\text{g},n}\right]}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}.(\left[I\right]-\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace \langle {g}^{\text{t},n}\rangle }{{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}^{2}})\]

La deuxième partie de l’expression est précédée du signe \(-\) , l’effet de cette contribution est particulièrement déstabilisant pour le comportement global de la matrice tangente au système, plus particulièrement lorsque l’on est loin de l’équilibre et donc au début de la résolution à chaque nouveau pas de temps. On décide donc de ne le prendre en compte que partiellement en l’affectant d’un coefficient \(\theta \in [0,1]\) que l’on peut modifier via le mot-clef COEF_MATR_FROT:

(2302)\[\left[{\tilde{K}}_{\theta}^{\text{g},n}\right]=\frac{{\left[{A}^{\text{g}}\right]}^{T}.\left[{k}^{\text{g},n}\right]}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}.(\left[I\right]-\theta .\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace \langle {g}^{\text{t},n}\rangle }{{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}^{2}})\]

On conseille d’utiliser une valeur initiale de \(0.5\) pour ce coefficient et de le diminuer si la convergence n’est pas obtenue. Dans le cas où \(\theta =0\) la convergence semble toujours être obtenue mais est particulièrement lente. Lorsque l’on est proche de la solution, il est par contre très utile d’avoir une valeur de ce coefficient égale à \(1\) de façon à accélérer la convergence. Cela est fait automatiquement dans le code lorsque le résidu RESI_GLOB_RELA est inférieur à \({10}^{-3}\) . On remplace \(\left[{K}^{\text{g},n}\right]\) par \(\left[{\tilde{K}}_{\theta}^{\text{g},n}\right]\) .

Le second membre vaut finalement:

(2303)\[\left\lbrace {L}^{\text{g},n}\right\rbrace ={\left[{A}^{\text{g}}\right]}^{T}.\left[{k}^{\text{g},n}\right].\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace }{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}\]

En 2D, le multiplicateur de Lagrange pourle glissement est approximé par:

(2304)\[{\left\lbrace {\mu}^{\text{g}}({u}_{i})\right\rbrace }_{\text{2D}}=\mu .\left[∣{\mu}^{\text{c}}({u}_{i})∣\right]=\left[{k}^{\text{g}}({u}_{i})\right]\]

On peut démontrer \({\left[{K}^{\text{g},n}\right]}_{\text{2D}}=0\) . Et donc:

(2305)\[{\left\lbrace {\stackrel{ˆ}{L}}_{i}^{\text{g},n+1}\right\rbrace }_{\text{2D}}={\left[{A}^{\text{g}}\right]}^{T}.\left\lbrace {\rho}^{\text{g},n}\right\rbrace\]

Algorithmes généraux

Les méthodes discrètes de résolution du problème de contact/frottement sont fondées sur une approche découplée entre l’équilibre et le contact/frottement.Lecontact/frottement (noté \(C\) dans la deuxième colonne) est traité après chaque itération de Newton du problème global (noté \(G\) dans la deuxième colonne).

Sans prise en compte du contact, on notera les vecteurs solutionsavec un \(\text{~}\) , par exemple:

(2306)\[\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\]

Cas de la dualisation

La procédure générale pour le cas dualiséest la suivante:

A l’itération de Newton \(n\)
1	G	Résolution du problème d’équilibre sans contact, équation (65)	\(\to \left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\)
2	G	Mise à jour des déplacements sans contact	\(\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\)
3	C	Modification des déplacements pour respecter les conditions de contact	\(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace \to \left\lbrace \delta {u}^{n}\right\rbrace\)
4	C	Mise à jour des déplacements avec prise en compte du contact	\(\left\lbrace {u}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace\)
5	C	Calcul des efforts de contact	\(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace\) et \(\left\lbrace {L}_{i}^{\text{f},n}\right\rbrace\)
6	C	Calcul éventuels de matrices (cas du glissement)	\(\left[{K}^{\text{g},n+1}\right]\)
7	G	Calcul des efforts intérieurs et extérieurs avec les déplacements modifiés	\(\left\lbrace {L}_{i}^{\text{ext},n}\right\rbrace\) et \(\left\lbrace {L}_{i}^{int,n}\right\rbrace\)
8	G	Vérification de l’équilibre

La matrice \(\left[{K}^{\text{cf},n+1}\right]\) n’est calculée que dans le cas du glissement. Elle ne servira pas nécessairement, tout dépendra de la convergence de Newton et de l’algorithme de contact/frottement (seuil de Tresca convergé sur le problème de Coulomb).

Cas de la régularisation

La procédure générale pour le cas régularisé (pénalisation)est la suivante:

A l’itération de Newton \(n\)
1	G	Résolution du problème d’équilibre sans contact, équation (65)	\(\to \left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\)
2	G	Mise à jour des déplacements sans contact	\(\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\)
3	C	Calcul des efforts de contact/frottement	\(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace\) et \(\left\lbrace {L}_{i}^{\text{f},n}\right\rbrace\)
4	C	Calcul des matrices modifiées	\(\left[{K}^{\text{cf},n+1}\right]\)
5	G	Calcul des efforts intérieurs et extérieurs	\(\left\lbrace {L}_{i}^{\text{ext},n}\right\rbrace\) et \(\left\lbrace {L}_{i}^{int,n}\right\rbrace\)
6	G	Vérification de l’équilibre

L’algorithme de contact/frottement ne modifie pas les déplacements, c’est le système global qui va les modifier. Ce qui implique que pour les méthodes pénalisées, on aura toujours besoin d’au moins deux itérations de Newton, même dans le cas élastique: une itération pour résoudre le problème d’équilibre sans contact/frottement et une seconde itération pour intégrer les conditions de contact-frottementqui vont modifierlesystème global.

[1]

Il existe plusieurs manières de choisir et calculer cette matrice. Pour plus de détails voir r5.03.01

[2]

L’idée de transformer le problème de Coulomb en une succession de problèmes de Tresca est très largement justifiée dans la littérature (voir [bib13] par exemple) et majoritairement utilisée dans les codes commerciaux. C’est aussi la stratégie retenue dans la version « point fixe » de l’algorithme en formulation continue (voir r5.03.52).

Résolution algorithmique – Contact sans frottement

Liaisons de contact

Chaque nœud esclave potentiellement en contact a un statut dont les algorithmes devront déterminer la nature. Un nœud esclave est appelé «liaison potentielle de contact» ou plus simplement «liaison». Le terme liaison faisant allusion au fait que la condition de contact est le résultat de l’imposition d’une relation cinématique sur les degrés de liberté de déplacement du nœud esclave. Ces liaisons sont réunies dans différents ensembles:

• :math:`` est l’ensemble des liaisons possibles (actives et non actives);

• \({\Xi}^{\text{nc}}\) est l’ensemble des nœuds esclaves qui ne sont pas en contact (liaisons non actives);

• \({\Xi}^{\text{c}}\) est l’ensemble des nœuds effectivement en contact (liaisons actives).

On a donc les relations suivantes entre les ensembles:

• \({\Xi}^{\text{c}}\cap {\Xi}^{\text{nc}}=\mathrm{\varnothing }\) car les nœuds sont en contact ou non;

• \({\Xi}^{l}={\Xi}^{\text{c}}\oplus {\Xi}^{\text{nc}}\) car les nœuds potentiellement en contact le sont ou pas;

Les méthodes dualisées en contact

Équilibre de la structure en présence de contact

On rappelle que l’équation d’équilibre en présence de contact s’écrit:

(2307)\[\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{c}}\right\rbrace\]

Dans le cas du contact dualisé, la force de contact s’écrit en fonction du multiplicateur de Lagrange de contact:

(2308)\[\left\lbrace {L}_{i}^{\text{c}}\right\rbrace ={\left\lbrace {L}_{i}^{\text{c}}\right\rbrace }_{\text{dual}}=[{A}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\]

Après linéarisation de l’équation d’équilibre (4476), on introduit la matrice tangente \(\left[{K}^{\text{m},n-1}\right]\) qui contiendra les contributions issues de la linéarisation des efforts intérieurs et extérieurs et la matrice \(\left[{K}^{\text{c},n-1}\right]\) pourles forces de contact, on trouve \(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\) , incrément de solution du problème d’équilibre à l’itération \(n\) mais sans application de la loi de contact:

(2309)\[(\left[{K}^{\text{m},n-1}\right]+\left[{K}^{\text{c},n-1}\right]).\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace =\left\lbrace {L}_{i}^{int,n-1}\right\rbrace -\left\lbrace {L}_{i}^{\text{ext},n-1}\right\rbrace +\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace\]

Nous avons vu que les forces de contact dans le cas dualisé ne dépendentpas du déplacement (voir r5.03.50-linearisation). On a donc \(\left[{K}^{\text{c},n-1}\right]=0\) . Ce qui nous permet d’exprimer l’équilibre de la structure avec prise en compte des forces de contact (en notant \(\left[K\right]=\left[{K}^{\text{m},n-1}\right]\) et \(\left\lbrace F\right\rbrace =\left\lbrace {L}_{i}^{\text{ext},n-1}\right\rbrace -\left\lbrace {L}_{i}^{int,n-1}\right\rbrace -\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace\) pour alléger):

(2310)\[\left[K\right].\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace =\left\lbrace F\right\rbrace\]

La solution \(\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace\) est obtenue après équilibre et avant application des conditions de contact. Elle s’écrit:

(2311)\[\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\]

Pour obtenir la solution \(\left\lbrace {\tilde{u}}_{i}^{n}\right\rbrace\) , nous n’avons pas appliqué les conditions de contact (loi de Signorini) pour l’itération de Newton courante \(n\) . Par contre, les réactions de contact calculées à l’itération de Newton précédente \(n-1\) sont bien prises en compte dans \(\left\lbrace F\right\rbrace\) . Pour résoudre complètement le problème d’équilibre avec contact/frottement, à l’itération \(n\) , il faut appliquer la loi de Signorini, qui s’exprime par le système suivant:

(2312)\[\begin{split}\lbrace \begin{array}{cc}\left[{A}^{\text{c}}\right].\left\lbrace {u}_{i}^{n}\right\rbrace \le \left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace & (a)\\ \left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \ge 0& (b)\\ \langle {\mu}_{i}^{\text{c}}\rangle .(\left[{A}^{\text{c}}\right].\left\lbrace {u}_{i}^{n}\right\rbrace )=0& (c)\end{array}\end{split}\]

On rappelle l’interprétation desconditions de Signorini:

• L’équation (a) représente les conditions géométriques de non interpénétration, l’inégalité s’entendant composante par composante (chaque ligne est relative à un couple potentiel de contact).

• L’équation (b) exprime l’absence d’opposition au décollement (les surfaces de contact ne peuvent connaître que des compressions), c’est la condition dite d’intensilité.

• L’équation (c) est la condition de compatibilité. Lorsque pour une liaison donnéele multiplicateur de Lagrange est non nul, on a contact et donc le jeuest nul. Lorsque le jeu est non nul (les deux surfaces ne sont pas en contact), le multiplicateur associé doit être nul (pas de compression).

L’application de la loi de Signorini va modifier l’incrément de déplacement \(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\) (qui devient \(\left\lbrace \delta {u}^{n}\right\rbrace\) ), d’où la solution obtenue après équilibre et après application de la loide contact qui s’écrit:

(2313)\[\left\lbrace {u}_{i}^{n}\right\rbrace =\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}^{n}\right\rbrace\]

Système réduit aux liaisons actives

Nous allons écrire le système permettant de résoudre complètement le problème d’équilibre avec prise en compte de la loi de Signorini. L’idée est de transformer les inégalités du système (4481) en égalités. On commence par évaluerle jeu donné par le déplacement calculé avant la correction du contact \(\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace\) , pour toutes les liaisons [1]_ :

(2314)\[{⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{l}}={⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace )⌋}_{{\Xi}^{l}}\]

On dit qu’une liaison \(J\) est active si sonjeu \({⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{\text{J}\in {\Xi}^{l}}\) est négatif, ce qui indique une interpénétration (la condition de contact n’est pas respectée), cette liaison devient donc une liaison de contact et permet donc de définir l’ensemble initialement postulé des liaisons de contact \({\Xi}^{\text{c}}\) :

(2315)\[{\Xi}^{\text{c}}=\left\lbrace J\in {\Xi}^{l}\mid {⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{J}<0\right\rbrace\]

On postule que, pour cesliaisons actives, le jeu effectif sera nul, et que donc l’inégalité \({⌊\left[{A}^{\text{c}}\right].\left\lbrace {u}_{i}^{n}\right\rbrace ⌋}_{{\Xi}^{l}}\le {⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace ⌋}_{{\Xi}^{l}}\) devient une égalité pour l’ensemble les liaisons actives:

(2316)\[{⌊\left[{A}^{\text{c}}\right].\left\lbrace {u}_{i}^{n}\right\rbrace =\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\]

Si on utilise \(\left\lbrace {d}^{\text{c},n-1}\right\rbrace\) , le jeu évalué avant l’itération de Newton courante:

(2317)\[\left\lbrace {d}^{\text{c},n-1}\right\rbrace =\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace )\]

Avec:

(2318)\[\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace =\left\lbrace {d}^{\text{c},n-1}\right\rbrace -\left[{A}^{\text{c}}\right].\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\]

L’égalité (4485) s’écrit finalement:

(2319)\[{⌊\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace =\left\lbrace {d}^{\text{c},n-1}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\]

Le système «mixte» inégalité/égalité (contenant l’équation d’équilibre et les conditions de Signorini)se transforme finalement en système simple qui ne traite que d’égalités, sur la base des liaisons de contact postulées \({\Xi}^{\text{c}}\) :

(2320)\[\begin{split}\lbrace \begin{array}{c}\left[K\right].\left\lbrace \delta {u}^{n}\right\rbrace +[{A}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace F\right\rbrace \\ {⌊\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace =\left\lbrace {d}^{\text{c},n-1}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\end{array}\end{split}\]

On cherche les inconnues \(\left\lbrace \delta {u}^{n}\right\rbrace\) et \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) qui sont solutions du système suivant:

(2321)\[\begin{split}\left[\begin{array}{cc}[K]& [{A}^{\text{c}}{]}^{T}\\ [{A}^{\text{c}}]& [0]\end{array}\right].\left\lbrace \begin{array}{c}\left\lbrace \delta {u}^{n}\right\rbrace \\ \left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \end{array}\right\rbrace =\left\lbrace \begin{array}{c}\left\lbrace F\right\rbrace \\ \left\lbrace {d}^{\text{c},n-1}\right\rbrace \end{array}\right\rbrace \mathrm{\iff }\left[{K}^{\text{c}}\right].\left\lbrace \begin{array}{c}\left\lbrace \delta {u}^{n}\right\rbrace \\ \left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \end{array}\right\rbrace =\left\lbrace \begin{array}{c}\left\lbrace F\right\rbrace \\ \left\lbrace {d}^{\text{c},n-1}\right\rbrace \end{array}\right\rbrace\end{split}\]

On remarquera que la formulation dualisée de la condition de non-interpénétration traduit la stationnarité du Lagrangien \(L\) ainsi défini:

(2322)\[L(\left\lbrace \delta {u}^{n}\right\rbrace ,\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )=\frac{1}{2}.\langle \delta {u}^{n}\rangle .\left[K\right].\left\lbrace \delta {u}^{n}\right\rbrace -\langle \delta {u}^{n}\rangle .\left\lbrace F\right\rbrace +\langle {\mu}_{i}^{\text{c}}\rangle .(\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace )\]

Cette stationnarité s’exprime sous la forme du problème de point-selle:

(2323)\[\underset{\left\lbrace \delta {u}^{n}\right\rbrace }{\min}\underset{\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \ge 0}{\max}L(\left\lbrace \delta {u}^{n}\right\rbrace ,\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )\]

La méthode des contraintes actives – ALGO_CONT=”CONTRAINTE”

Il s’agit de généraliser la démarche du paragraphe précédent en adoptant un algorithme itératif qui travaille en trois phases:

1. Faire des hypothèses sur l’état des liaisons (transformations inégalités → égalité);

2. Résoudre le nouveau système ainsi créé;

3. Vérifier les hypothèses initiales et boucler (éventuellement) en 1.

On pourra trouver une description complète de la méthode avec les justifications théoriques nécessaires dans [bib1] et [bib2] . Le principe est le suivant: on postule un ensemble de contraintes dites actives, qui correspondent à un jeu nul (la relation inégalité devient une égalité) ; on résout le système d’équations obtenu dans ce sous-espace, et on regarde si le postulat de départ était justifié. Si l’ensemble choisi était trop petit (on avait oublié des liaisons actives), on rajoute à l’ensemble la liaison violant le plus la condition de non-interpénétration ; si l’ensemble choisi était trop grand (des liaisons supposées actives ne le sont en fait pas), on enlève de l’ensemble la liaison la plus improbable i.e. celle dont le multiplicateur de Lagrange viole le plus la condition d’intensilité. Le fait d’enlever ou de rajouter une seule liaison à chaque itération de la méthode garantit la convergence en un nombre finid’itérations inférieur ou égal à deux fois le nombre maximum de liaisons, pour peu que la matrice \(\left[K\right]\) soit symétrique et définie positive [bib1], [bib15], [bib16].

En élasticité, à la fin des itérations de contraintes actives, on a un résultat convergé au sens de Newton. En plasticité ou si on réactualise la géométrie, ce n’est pas le cas car plusieurs itérations de Newton sont nécessaires pour obtenir l’équilibre. Après chaque itération de Newton, on lance l’algorithme de contraintes actives pour satisfaire les conditions de contact. Ainsi, en élasticité, on convergera nécessairement pour chaque pas en une itération si REAC_GEOM=”SANS”.

Écriture du problème itératif

On part de l’incrément obtenu sans traiter le contact \(\left\lbrace \delta {\stackrel{~}{u}}^{n}\right\rbrace\) et on effectue les itérations de contraintes actives \(k\) jusqu’à convergence propre de cet algorithme. La convergence au sens des contraintes actives est obtenue lorsqu’aucune liaison ne viole la condition cinématique (inégalité ((4481) a))et lorsque les multiplicateurs de Lagrange associés sont tous positifs (inégalité ((4481) b)).

On note \(k\) les itérations de contraintes actives. La solution de départ sans correction du contact est \(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\) et l’incrément ajouté par la nouvelle itération de contact est \(\left\lbrace {\delta}_{k}\right\rbrace\) . On note:

(2324)\[ \begin{align}\begin{aligned}\left\lbrace \delta {u}_{0}^{n}\right\rbrace =\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\\\left\lbrace {\delta}_{0}\right\rbrace =\left\lbrace 0\right\rbrace\\\left\lbrace \delta {u}_{k}^{n}\right\rbrace =\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\left\lbrace {\delta}_{k}\right\rbrace\end{aligned}\end{align} \]

On se place à l’itération de contact \(k\) . On cherche à résoudre le système (4489):

(2325)\[\begin{split}\left\{ \begin{array}{c}\left[K\right].\left\lbrace \delta {u}_{k}^{n}\right\rbrace +[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace F\right\rbrace \\ {⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {u}_{k}^{n}\right\rbrace =\left\lbrace {d}^{\text{c},n-1}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\end{array} \right.\end{split}\]

On note que le multiplicateur de Lagrange du contact \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) n’est pas résolu de manière incrémentale mais de manière totale (sur le pas de temps). En injectant (4493) dans ce système, on obtient:

(2326)\[\begin{split}\left\{ \begin{array}{c}\left[K\right].\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\left[K\right].\left\lbrace {\delta}_{k}\right\rbrace +[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace F\right\rbrace \\ {⌊\left[{A}_{k}^{\text{c}}\right].(\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\left\lbrace {\delta}_{k}\right\rbrace )=\left\lbrace {d}^{\text{c},n-1}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\end{array} \right.\end{split}\]

La première équation peut se ré-écrire:

(2327)\[\left\lbrace {\delta}_{k}\right\rbrace ={\left[K\right]}^{-1}.\left\lbrace F\right\rbrace -\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace -{\left[K\right]}^{-1}.[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\]

En tenant compte que \(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace ={\left[K\right]}^{-1}.\left\lbrace F\right\rbrace\) , on simplifie:

(2328)\[{\delta}_{k}=\delta {\tilde{u}}^{n}-\delta {u}_{k-1}^{n}-{K}^{-1}.{[{A}_{k}^{c}]}^{T}.{\mu}_{i}^{c}\]

A partir de l’état postulé des liaisons actives:

(2329)\[{⌊\left[{A}_{k}^{\text{c}}\right].(\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\left\lbrace {\delta}_{k}\right\rbrace )=\left\lbrace {d}^{\text{c},n-1}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\]

Et en reprenant l’expression de \(\left\lbrace {d}^{\text{c},n-1}\right\rbrace\) donnée par (69):

(2330)\[{⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace ⌋}_{{\Xi}_{C}^{k}}={⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace )⌋}_{{\Xi}_{k}^{\text{c}}}\]

On utilise l’expression de \(\left\lbrace {\delta}_{k}\right\rbrace\) donnée par (2482):

(2331)\[{⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].{\left[K\right]}^{-1}.[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}={⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace )⌋}_{{\Xi}_{k}^{\text{c}}}\]

Au final, les multiplicateurs de Lagrange \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) pour les liaisons de contact sont solutions du système suivant:

(2332)\[{⌊-\left[{A}_{k}^{\text{c}}\right].{\left[K\right]}^{-1}.[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace {\tilde{d}}_{\text{c},n}^{}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\]

On remarque que (2487) est le complément de Schur \({\left[K\right]}_{\text{schur}}=\left[{A}_{k}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}_{k}^{\text{c}}\right]}^{T}\) de la matrice \(\left[{K}^{\text{c}}\right]\) . On peut finalement calculer les incréments de déplacement \(\left\lbrace {\delta}_{k}\right\rbrace\) avec:

(2333)\[\left\lbrace {\delta}_{k}\right\rbrace =\sum_{l=0}^{k-1}\left\lbrace {\delta}_{l}\right\rbrace -{\left[K\right]}^{-1}.{[{A}_{k}^{c}]}^{T}.{\left\lbrace {\mu}_{i}\right\rbrace }^{c}\]

La résolution de (2487) est la partie la plus coûteuse en temps calcul de l’algorithme. Toute l’efficacité de la stratégie consiste à utiliser le fait qu’on résout ce système uniquement sur l’ensemble des liaisons actives \({\Xi}_{k}^{\text{c}}\) , on utilise donc deux propriétés.

• Le calcul du complément de Schur \({\left[K\right]}_{\text{schur}}\) fait toujours appel à une factorisation de type \({\mathit{LDL}}^{T}\) , ce qui permet de gagner du temps, car une telle factorisation possède la propriété remarquable d’être incrémentale, c’est-à-dire que l’ajout d’une liaisonn’oblige pasà reconstruire la factorisée depuisle début mais seulement la partie qu’on modifie.

• Le résolution (descente/remontée) ne se fait aussi que sur l’ensemble des liaisons actives \({\Xi}_{k}^{\text{c}}\) .

Validité de l’ensemble de liaisons actives choisi

A la fin de chaque itération de contrainte active, il convient de vérifier si l’ensemble \({\Xi}_{k}^{\text{c}}\) est bien correct. Soit la liaison \(J\in {\Xi}^{l}\) , trois situations sont possibles:

1. Le déplacement relatif compense le jeu initial \({⌊\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}_{k}^{n}\right\rbrace )⌋}_{J\in {\Xi}^{l}}={⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace ⌋}_{J\in {\Xi}^{l}}\) ;

2. Le déplacement relatif est inférieur au jeu initial \({⌊\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}_{k}^{n}\right\rbrace )⌋}_{J\in {\Xi}^{l}}<{⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace ⌋}_{J\in {\Xi}^{l}}\) ;

3. Le déplacement relatif est supérieur au jeu initial \({⌊\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}_{k}^{n}\right\rbrace )⌋}_{J\in {\Xi}^{l}}>{⌊\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace ⌋}_{J\in {\Xi}^{l}}\) .

La situation \((3)\) est interdite car elle correspond à une violation de la condition de non-interpénétration. La situation \((1)\) correspond à une liaison dite active, la situation \((2)\) à une liaison non-active.

Au début de l’itération \(k\) de l’algorithme, on avait postulé un ensemble de liaisons actives \({\Xi}_{k}^{\text{c}}\) . On a trouvé un incrément possible \(\left\lbrace {\delta}_{k}\right\rbrace\) des inconnues sous ces hypothèses eton va maintenant vérifier que cet incrément est compatible avec les hypothèses. En pratique, cela consiste à faire deux vérifications:

1. L’ensemble \({\Xi}_{k}^{\text{c}}\) est-il trop petit? On vérifie que liaisons supposées non actives ne violent pas la condition de non-interpénétration, sinon on en active une, celle qui viole le plus la condition.

2. L’ensemble \({\Xi}_{k}^{\text{c}}\) est-il trop grand? On vérifie que les liaisons supposées actives sont associées à des multiplicateurs de contact \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) positifs ou nuls, sinon on en désactive celle qui viole le plus la condition.

L’ensemble \({\Xi}_{k}^{\text{c}}\) est-il trop petit?

Pour vérifier que l’ensemble \({\Xi}_{k}^{\text{c}}\) n’est pas trop petit, on va calculer pour toutes les liaisons supposées non-actives, la quantité suivante:

(2334)\[{\rho}_{J}={⌊\frac{\left\lbrace {d}_{\text{ini}}^{\text{c}}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].(\left\lbrace {u}_{i-1}\right\rbrace +\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {u}_{k}^{n}\right\rbrace )}{\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace }⌋}_{J\in {\Xi}_{k}^{\text{nc}}}={⌊\frac{\left\lbrace {d}^{\text{c},n-1}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {u}_{k}^{n}\right\rbrace }{\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace }⌋}_{J\in {\Xi}_{k}^{\text{nc}}}\]

Il y a deux cas:

1. Si \({⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace ⌋}_{J\in {\Xi}_{k}^{\text{nc}}}<0\) , le jeu pour la liaison \(J\) va augmenter, et donc la liaison supposée non-active reste dans cet état, \({\rho}_{J}\) est strictement supérieur à \(1\) ;

2. Si \({⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace ⌋}_{J\in {\Xi}_{k}^{\text{nc}}}>0\) , le jeu pour la liaison \(J\) va diminuer, et donc la liaison supposée non-active va être activée, \({\rho}_{J}\) est inférieur ou égal à \(1\) ;

On examine donc \(\stackrel{ˉ}{\rho}=\underset{J}{\text{Min}}{\rho}_{j}\) sur l’ensemble des liaisons \(J\) déclarées non actives. Si \(\stackrel{ˉ}{\rho}<1\) , cela indique qu’une liaison au moins est violée (situation \((3)\) ): on rajoute alors à la liste des liaisons actives le numéro de la liaison dont l’interpénétration est la plus grande, c’est-à-dire celle qui réalise le minimum de \({\rho}_{J}\) et on écrit \(\left\lbrace \delta {u}_{k+1}^{n}\right\rbrace =\left\lbrace \delta {u}_{k}^{n}\right\rbrace +\stackrel{ˉ}{\rho}.\left\lbrace {\delta}_{k}\right\rbrace\) (ce quicorrespond à un jeu nul pour la liaison rajoutée).

On présente l’algorithme utilisé:

Si \({\Xi}_{k}^{\text{c}}={\Xi}^{l}\)

\(\stackrel{ˉ}{\rho}=1\)

Sinon

\(\stackrel{ˉ}{\rho}=1\)

Boucle sur les liaisons \(J\in {\Xi}_{k}^{\text{nc}}\)

Calcul de \({\alpha}_{J}={⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace {\delta}_{k}\right\rbrace ⌋}_{J}\)

Si \({\alpha}_{J}<0\)

\({\rho}_{J}=\frac{{⌊\left\lbrace {d}^{\text{c},n-1}\right\rbrace -\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {u}_{k}^{n}\right\rbrace ⌋}_{J}}{{\alpha}_{J}}\)

\(\stackrel{ˉ}{\rho}=\min({\rho}_{J},\stackrel{ˉ}{\rho})\)

\({J}_{min}=J\)

Fin Boucle

En sortie de cet algorithme, on aura le numéro de la liaison qui viole le plus la condition de non-pénétration \({J}_{min}\) et la valeur \(\stackrel{ˉ}{\rho}\) : \(\left\lbrace {J}_{min},\stackrel{ˉ}{\rho}\right\rbrace =\mathit{IsPetit}({\Xi}_{k}^{\text{c}})\)

L’ensemble \({\Xi}_{k}^{\text{c}}\) est-il trop grand?

La seconde vérification consiste à se demander si l’ensemble des liaisons actives est trop grand. On se place maintenant dans le cas où \(\stackrel{ˉ}{\rho}\ge 1\) , c’est-à-dire que l’on sait que \({\Xi}_{k}^{\text{c}}\) n’est pas trop petit. Alors:

• Si aucune liaison n’est active, la méthode a convergé vers un état sans contact;

• S’il y a des liaisons supposées actives:

• Si tous les multiplicateurs de Lagrange \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) sont positifs ou nuls, on a également convergé vers un état avec contact effectif;

• S’il existe des multiplicateurs de Lagrange \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) négatifs, les liaisons correspondantes ne devraient pas être actives: on retire de l’ensemble des liaisons actives la liaison dont le multiplicateur négatif est le plus grand en valeur absolue.

En sortie de cet algorithme, on aura le numéro de la liaison qui viole le plus la condition d’intensilité: \(\left\lbrace {J}_{max}\right\rbrace =\mathit{IsGrand}({\Xi}_{k}^{\text{c}})\) .

Algorithme

Au final, l’algorithme pour les contraintes actives est le suivant:

\(\mathit{Ini}\) \({\mathit{\delta u}}_{0}^{n}=\delta {\tilde{u}}^{n}\)

Calcul de \({⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{l}}\)

Évaluation de \({\Xi}_{0}^{\text{c}}=\left\lbrace J\in {\Xi}^{l}\mid {⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{J}<0\right\rbrace\) et \(\left[{A}_{0}^{\text{c}}\right]\)

\({B}_{k}\) Boucle sur les contraintes actives \(k=1,{\mathit{Iter}}_{max}\)

\(\left\lbrace {\delta}_{k-1}\right\rbrace =\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace -\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace\)

\({\Xi}_{k}^{\text{c}}\leftarrow {\Xi}_{k-1}^{\text{c}}\) et \(\left[{A}_{k}^{\text{c}}\right]\leftarrow \left[{A}_{k-1}^{\text{c}}\right]\)

Si \({\Xi}_{k}^{\text{c}}\mathrm{\ne }\mathrm{\varnothing }\)

Calcul de \({\left[{K}_{k}\right]}_{\text{schur}}=\left[{A}_{k}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}_{k}^{\text{c}}\right]}^{T}\)

Factorisation de \({⌊{\left[{K}_{k}\right]}_{\text{schur}}⌋}_{{\Xi}_{k}^{\text{c}}}\)

Résolution de \({⌊{\left[{K}_{k}\right]}_{\text{schur}}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace {\tilde{d}}_{\text{c},n}^{}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\to {⌊\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\)

\(\left\lbrace {\delta}_{k}\right\rbrace =\left\lbrace {\delta}_{k-1}\right\rbrace -{\left[K\right]}^{-1}.[{A}_{k}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\)

\(\left\lbrace {J}_{min},\stackrel{ˉ}{\rho}\right\rbrace =\mathit{IsPetit}({\Xi}_{k}^{\text{c}})\)

\(\left\lbrace \delta {u}_{k}^{n}\right\rbrace =\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace +\stackrel{ˉ}{\rho}.\left\lbrace {\delta}_{k}\right\rbrace\)

Si \(\stackrel{ˉ}{\rho}<1\)

\({\Xi}_{k}^{\text{c}}\leftarrow {\Xi}_{k}^{\text{c}}+\left\lbrace {J}_{min}\right\rbrace\)

Sinon

Si \({\Xi}_{k}^{\text{c}}=\mathrm{\varnothing }\)

Goto Fin \(F\)

\(\left\lbrace {J}_{max}\right\rbrace =\mathit{IsGrand}({\Xi}_{k}^{\text{c}})\)

Si \(\left\lbrace {J}_{max}\right\rbrace =0\)

Goto Fin \(F\)

Sinon

\({\Xi}_{k}^{\text{c}}\leftarrow {\Xi}_{k}^{\text{c}}-\left\lbrace {J}_{max}\right\rbrace\)

\({B}_{k}\) Boucle \(k=k+1\)

\(F\) FIN

Calcul de \(\left\lbrace {L}_{i}^{\text{c}}\right\rbrace =[{A}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\)

La méthode du Gradient Conjugué Projeté – ALGO_CONT=”GCP”

La méthode de résolution présentée dans cette partie est une application de l’algorithme du Gradient Conjugué Projeté. Il s’agit très précisément de la version itérative de la méthode des Contraintes Actives présentée dans la partie précédente.

Reformulation du problème de contact

On rappelle que le système à résoudre à chaque itération de Newton \(n\) est le suivant:

(2335)\[\begin{split}\lbrace \begin{array}{c}\left[K\right].\left\lbrace \delta {u}^{n}\right\rbrace +[{A}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace =\left\lbrace F\right\rbrace \\ \left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace \le \left\lbrace {d}^{\text{c},n-1}\right\rbrace \end{array}\end{split}\]

Il provient de la dualisation des conditions de contact et il traduit la stationnarité du Lagrangien \(L\) ainsi défini:

(2336)\[L(\left\lbrace \delta {u}^{n}\right\rbrace ,\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )=\frac{1}{2}.\langle \delta {u}^{n}\rangle .\left[K\right].\left\lbrace \delta {u}^{n}\right\rbrace -\langle \delta {u}^{n}\rangle .\left\lbrace F\right\rbrace +\langle {\mu}_{i}^{\text{c}}\rangle .(\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace )\]

Cette stationnarité peut s’exprimer sous la forme du problème de point-selle:

(2337)\[\underset{\left\lbrace \delta {u}^{n}\right\rbrace }{\min}\underset{\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \ge 0}{\max}L(\left\lbrace \delta {u}^{n}\right\rbrace ,\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )\]

Or le minimum de \(L\) par rapport à \(\delta {u}_{i}\) est connu; il a pour expression:

(2338)\[ \begin{align}\begin{aligned}\left\lbrace \delta {u}^{n}\right\rbrace ={\left[K\right]}^{-1}.(\left\lbrace F\right\rbrace -{\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )\\\textrm{et}\\\langle \delta {u}^{n}\rangle =(\langle F\rangle -\langle {\mu}_{i}^{\text{c}}\rangle .\left[{A}^{\text{c}}\right]).{\left[K\right]}^{-1}\end{aligned}\end{align} \]

On n’a donc plus qu’à faire la maximisation de \(L\) par rapport à \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \ge 0\) . Sachant que maximiser une fonctionnelle \(J\) équivaut à minimiser \(-J\) , on se ramène au problème de minimisation suivant:

(2339)\[\underset{\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace \ge 0}{\min}H(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )\]

Où la fonctionnelle \(H\) s’écrit:

(2340)\[H(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace )=\frac{1}{2}.\langle {\mu}_{i}^{\text{c}}\rangle .\left[{A}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace +\langle {\mu}_{i}^{\text{c}}\rangle .(\left\lbrace {d}^{\text{c},n-1}\right\rbrace -\left[{A}^{\text{c}}\right].\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace )+\frac{1}{2}.\langle F\rangle .{\left[K\right]}^{-1}.\left\lbrace F\right\rbrace\]

Cette expression est la forme duale du problème de contact: elle fait intervenir le champ de multiplicateurs de Lagrange \(\left\lbrace {\mu}_{i}^{\text{c}}\right\rbrace\) et plus du tout le champ de déplacement \(\left\lbrace \delta {u}^{n}\right\rbrace\) . Le problème à résoudre est maintenant une minimisation sous contrainte de positivité de l’inconnue. La méthode du Gradient Conjugué Projeté est une méthode simple et efficace pour ce type de problème. On notera:

• \(\left\lbrace Z\right\rbrace\) la direction de recherche

• \(\left\lbrace r\right\rbrace\) le sous-gradient et \(\left\lbrace {r}^{p}\right\rbrace\) sa version pré-conditionnée

Recherche linéaire

Dans un algorithme de gradient conjugué, il est nécessaire d’estimer un pas d’avancement. Deux variantes de recherche linéaire sont disponibles, admissible ou non. Elles sont représentées graphiquement sur la figure . La variante admissible, qui impose de rester dans le domaine convexe admissible induit une unique résolution par itération; elle conduit à une convergence assez régulière. La variante non-admissible, qui autorise à sortir du domaine admissible pour s’y reprojeter ensuite, induit deuxrésolutions par itération; elle conduit à une convergence assez erratique mais généralement plus rapide que la méthode admissible.

Figure 5.2.4.2-a: variantes de la recherche linéaire.

Voici l’algorithme de la recherche linéaire:

\(\mathit{RechLine}\)

\(\mathit{Ini}\) État de contact \(\left[{A}^{\text{c}}\right]\leftarrow \left[{A}_{k}^{\text{c}}\right]\) et \({\Xi}^{\text{c}}\leftarrow {\Xi}_{k}^{\text{c}}\)

Calcul du second membre \(\left\lbrace F\right\rbrace ={\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace Z\right\rbrace\)

Résolution de \(\left[K\right].\left\lbrace \delta a\right\rbrace =\left\lbrace F\right\rbrace \to \left\lbrace \delta a\right\rbrace\)

Calcul de \(\alpha =\frac{\langle r\rangle .\left\lbrace {r}^{p}\right\rbrace }{\langle \delta a\rangle .\left\lbrace F\right\rbrace }\)

Si ADMISSIBLE

\(\forall J\in {\Xi}^{l}:\text{Si}{\left\lbrace Z\right\rbrace }_{J}<0\text{alors}\alpha =\underset{J\in {\Xi}^{l}}{\min}(\alpha ,-\frac{{\left\lbrace {\mu}^{\text{c}}\right\rbrace }_{J}}{{\left\lbrace Z\right\rbrace }_{J}})\)

\(\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace =\left\lbrace {\mu}_{k}^{\text{c}}\right\rbrace +\alpha .\left\lbrace Z\right\rbrace\)

\(\left\lbrace \delta {u}_{k+1}\right\rbrace =\left\lbrace \delta {u}_{k}\right\rbrace -\alpha .\left\lbrace \delta a\right\rbrace\)

Sinon

\(\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace =\left\lbrace {\mu}_{k}^{\text{c}}\right\rbrace +\alpha .\left\lbrace Z\right\rbrace\)

\(\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace \leftarrow \max(\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace ,0)\)

Calcul du second membre \(\left\lbrace F\right\rbrace ={\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace {\mu}^{\text{c}}\right\rbrace\)

Résolution de \(\left[K\right].\left\lbrace \delta a\right\rbrace =\left\lbrace F\right\rbrace \to \left\lbrace \delta a\right\rbrace\)

\(\left\lbrace \delta {u}_{k+1}\right\rbrace =\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace -\left\lbrace \delta a\right\rbrace\)

\(F\) FIN

En sortie de cet algorithme, on aura le déplacement et les multiplicateurs de Lagrange après calcul du pas d’avancement: \(\left\lbrace \left\lbrace \delta {u}_{k+1}\right\rbrace ,\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace \right\rbrace =\mathit{RechLine}(\left\lbrace Z\right\rbrace ,\left\lbrace r\right\rbrace ,\left\lbrace {r}^{p}\right\rbrace )\) . Cet algorithme comprend deux résolutions mais en utilisant la matrice \(\left[K\right]\) déjà factorisée (c’est la matrice globale du problème d’équilibre), on ne fait donc qu’une descente-remontée, ce qui coûte peu cher.

La variante non-admissible est obligée de recalculer un incrément de déplacement compatible avec la contrainte de positivité sur les multiplicateurs de contact.

Pré-conditionnement

Dans l’algorithme du GCP, il est fait mention d’un appel optionnel à un préconditionneur. Le but de ce préconditionneur est d’accélérer la convergence de la méthode. Sa définition provient des considérations d’analyse fonctionnelle suivantes: dans la phase de mise à jour \(\left\lbrace {\mu}_{k+1}^{\text{c}}\right\rbrace =\left\lbrace {\mu}_{k}^{\text{c}}\right\rbrace +\alpha .\left\lbrace {Z}_{k}\right\rbrace\) et si l’on n’a pas fait appel au préconditionneur, le champ \(\left\lbrace {\mu}_{k}^{\text{c}}\right\rbrace\) appartient à \({H}^{-1/2}({\Gamma}_{\text{c}})\) tandis que \(\left\lbrace {Z}_{k}\right\rbrace\) appartient à \({H}^{1/2}({\Gamma}_{\text{c}})\) . Cette somme n’a donc pas de sens mathématique. Pour lui en donner un, il faut envoyer \(\left\lbrace {Z}_{k}\right\rbrace\) dans \({H}^{-1/2}({\Gamma}_{\text{c}})\) , c’est l’opération de préconditionnement. Sachant que \(\left\lbrace {Z}_{k}\right\rbrace\) est obtenu par l’expression \(\left\lbrace {Z}_{k}\right\rbrace =\left\lbrace {r}_{k}^{p}\right\rbrace +{\gamma}_{k}.\left\lbrace {Z}_{k-1}\right\rbrace\) , c’est sur \(\left\lbrace {r}_{k}^{p}\right\rbrace\) que l’on va opérer. C’est ce qui est fait par le préconditionneur en résolvant le problème auxiliaire suivant:

(2341)\[\begin{split}\left\{ \begin{array}{c}{⌊\left[K\right].\left\lbrace a\right\rbrace +{\left[{A}_{k}^{\text{c}}\right]}^{T}.\left\lbrace {r}_{k}^{p}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}=0\\ {⌊\left[{A}_{k}^{\text{c}}\right].\left\lbrace a\right\rbrace =\left\lbrace {r}_{k}^{p}\right\rbrace ⌋}_{{\Xi}_{k}^{\text{c}}}\end{array}\right.\end{split}\]

On résout donc un problème de déplacement imposé sur la partie où le contact est effectif et l’on récupère les réactions à l’encastrement \(\left\lbrace {r}_{k}^{p}\right\rbrace\) qui appartiennent bien à \({H}^{-1/2}({\Gamma}_{\text{c}})\) . Dans la terminologie de la décomposition de domaine, il s’agit là d’un préconditionneur de Dirichlet.

Pour résoudre ce problème auxiliaire, on utilise aussi un algorithme de gradient conjuguémais sans projection car aucune contrainte de positivité n’y apparaît. Cette approche itérative autorise une résolution approchée de manière à économiser du temps de calcul.

• \(\left\lbrace X\right\rbrace\) la direction de recherche;

• \(\left\lbrace s\right\rbrace\) le sous-gradient.

Voici l’algorithme du pré-conditionneur:

\(\mathit{PreCond}\)

\(\mathit{Ini}\) \(\left\lbrace {a}_{0}\right\rbrace =0\)

État de contact \(\left[{A}^{\text{c}}\right]\leftarrow \left[{A}_{k}^{\text{c}}\right]\) et \({\Xi}^{\text{c}}\leftarrow {\Xi}_{k}^{\text{c}}\)

Si \({\Xi}^{\text{c}}=\mathrm{\varnothing }\)

Goto Fin \(F\)

\({B}_{p}\) Boucle sur \(p=1,{\mathit{Iter}}_{max}\)

Calcul du gradient \(\left\lbrace {s}_{p}\right\rbrace =\left[{A}^{\text{c}}\right].\left\lbrace {a}_{p-1}\right\rbrace -\left\lbrace r\right\rbrace\)

Calcul du résidu \(\varepsilon =\underset{J\in {\Xi}^{\text{c}}}{max}\left\lbrace {s}_{p}\right\rbrace\)

Si \(\varepsilon <{\varepsilon}_{\mathit{precond}}\)

Goto Fin \(F\)

Si \(p=1\) ou réactualisation

\(\left\lbrace {X}_{p}\right\rbrace =\left\lbrace {s}_{p}\right\rbrace\)

Sinon, conjugaison

\(\beta =\frac{\langle {s}_{p}\rangle .\left\lbrace {s}_{p}\right\rbrace }{\langle {s}_{p-1}\rangle .\left\lbrace {s}_{p-1}\right\rbrace }\)

\(\left\lbrace {X}_{p}\right\rbrace =\left\lbrace {s}_{p}\right\rbrace +\beta .\left\lbrace {X}_{p-1}\right\rbrace\)

Calcul du second membre \(\left\lbrace F\right\rbrace ={⌊{\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace {X}_{p}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\)

Résolution de \(\left[K\right].\left\lbrace \delta a\right\rbrace =\left\lbrace F\right\rbrace \to \left\lbrace \delta a\right\rbrace\) | Calcul du complément de Schur \({\left\lbrace \delta a\right\rbrace }_{\text{schur}}={⌊\left[{A}^{\text{c}}\right].\left\lbrace \delta a\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\)

Calcul du pas d’avancement \(\alpha =\frac{\langle {s}_{p}\rangle .\left\lbrace {s}_{p}\right\rbrace }{\langle F\rangle .{\left\lbrace \delta a\right\rbrace }_{\text{schur}}}\)

Nouveau sous-gradient \(\left\lbrace {s}_{p}\right\rbrace =\left\lbrace {s}_{p-1}\right\rbrace +\alpha .\left\lbrace {X}_{p}\right\rbrace\)

Nouveau déplacement \(\left\lbrace {a}_{p}\right\rbrace =\left\lbrace {a}_{p-1}\right\rbrace -\alpha .\left\lbrace \delta a\right\rbrace\)

\({B}_{k}\) Fin Boucle \(k=k+1\)

\(F\) FIN

\(\left\lbrace {r}_{k}^{p}\right\rbrace =\left\lbrace {s}_{p}\right\rbrace\)

Projection: \(({max}_{J}({r}_{k}^{p},0)\mathit{si}{\mu}_{J}<0)\)

En sortie de cet algorithme, on aura sous-gradient pré-conditionné : \(\left\lbrace {r}_{k}^{p}\right\rbrace =\mathit{PreCond}(\left\lbrace {r}_{k}\right\rbrace )\) .

Algorithme

Voici l’algorithme global pour la méthode GCP:

\(\mathit{Ini}\) Si \(\langle {\mu}_{i-1}^{\text{c}}\rangle \cdot \left\lbrace {\mu}_{i-1}^{\text{c}}\right\rbrace \mathrm{\ne }0\) alors \([K].\left\lbrace v\right\rbrace =[{A}_{i-1}^{\text{c}}{]}^{T}.\left\lbrace {\mu}_{i-1}^{\text{c}}\right\rbrace \to \left\lbrace v\right\rbrace\)

Si \(\langle {\mu}_{i-1}^{\text{c}}\rangle \cdot \left\lbrace {\mu}_{i-1}^{\text{c}}\right\rbrace =0\) alors \(\left\lbrace v\right\rbrace =\left\lbrace 0\right\rbrace\)

Calcul de \({\mathit{\delta u}}_{0}^{n}=\delta {\stackrel{̃}{u}}^{n}-v\)

Évaluation de \({\Xi}_{0}^{\text{c}}=\left\lbrace J\in {\Xi}^{l}\mid {⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{J}<0\right\rbrace\)

\({B}_{k}\) Boucle sur \(k=1,{\mathit{Iter}}_{max}\)

\({\Xi}_{k}^{\text{c}}\leftarrow {\Xi}_{k-1}^{\text{c}}\)

Calcul du sous-gradient \(\left\lbrace {r}_{k}\right\rbrace =(\left[{A}_{k}^{\text{c}}\right].\left\lbrace \delta {u}_{k-1}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace )\)

Projection du sous-gradient \(\left\lbrace {r}_{k}\right\rbrace =\underset{J\in {\Xi}^{l}}{max}(\left\lbrace {r}_{k}\right\rbrace ,0)\)

Calcul du résidu \(\varepsilon =\underset{J\in {\Xi}^{l}}{max}\left\lbrace {r}_{k}\right\rbrace\)

Préconditionnement \(\left\lbrace {r}_{k}^{p}\right\rbrace =\mathit{PreCond}(\left\lbrace {r}_{k}\right\rbrace )\)

Conjugaison

Calcul du coefficient \({\gamma}^{k}=\frac{\langle {r}_{k}\rangle .\left\lbrace {r}_{k}^{p}\right\rbrace -\langle {r}_{k}\rangle .\left\lbrace {r}_{k-1}^{p}\right\rbrace }{\langle {r}_{k-1}\rangle .\left\lbrace {r}_{k-1}^{p}\right\rbrace }\)

\(\left\lbrace {Z}_{k}\right\rbrace =\left\lbrace {r}_{k}\right\rbrace +{\gamma}^{k}.\left\lbrace {Z}_{k-1}\right\rbrace\)

Recherche linéaire \(\left\lbrace \left\lbrace \delta {u}_{k+1}\right\rbrace ,\left\lbrace {\mu}^{\text{c}}\right\rbrace \right\rbrace =\mathit{RechLine}(\left\lbrace Z\right\rbrace ,\left\lbrace r\right\rbrace ,\left\lbrace {r}^{p}\right\rbrace )\)

\({B}_{k}\) Boucle \(k=k+1\)

\(F\) FIN

On peut faire les commentaires suivants sur l’algorithme:

• Lors de la recherche linéaire, la résolution du système \(\left[K\right].\left\lbrace \delta a\right\rbrace ={\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace Z\right\rbrace\) et le calcul du terme \(\langle \delta a\rangle .{\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace Z\right\rbrace\) , on perçoit que l’algorithme présenté est la version itérative de la méthode des contraintes actives. En effet, si l’on explicite le terme \(\langle \delta a\rangle .{\left[{A}^{\text{c}}\right]}^{T}.\left\lbrace Z\right\rbrace\) , on obtient \(\langle Z\rangle .(\left[{A}_{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}_{\text{c}}\right]}^{T}).\left\lbrace Z\right\rbrace\) . On retrouve le complément de Schur \({\left[K\right]}_{\text{schur}}=\left[{A}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}^{\text{c}}\right]}^{T}\) qui est explicitement construit dans la méthode des contraintes actives. Comme dans toutes les méthodes itératives, le Gradient Conjugué Projeté ne construit pas cet opérateur mais utilise son effet sur un vecteur (ici sur \(\left\lbrace Z\right\rbrace\) );

• Pour être vraiment efficace, cet algorithme doit être utilisé avec un solveur direct ou un solveur itératif préconditionné par la méthode LDLT_SP.

  • Avec un solveur direct, la matrice du système étant déjà factorisée, chaque résolution se résume à une descente-remontée:

  • Avec un solveur itératif préconditionné par la méthode LDLT_SP, chaque résolution se résume à quelques itérations du solveur itératif.

• On conjugue à chaque itération.

La méthode pénalisée en contact –ALGO_CONT=”PENALISATION”

La méthode pénalisée est particulièrement appropriée lorsque la matrice \(\left[K\right]\) est non symétrique. Néanmoins, elle nécessite une étude de sensibilité au coefficient de pénalisation (de dimension d’un rapport d’une force par un déplacement) et a minima une vérification sur les résultats que le contact est correctement pris en compte (les structures ne s’interpénètrent pas).

Équilibre de la structure en présence de contact

On rappelle que l’équation d’équilibre en présence de contact s’écrit:

(2342)\[\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{c}}\right\rbrace\]

Dans le cas du contact pénalisé, la force de contact s’écrit:

(2343)\[\left\lbrace {L}_{i}^{\text{c}}\right\rbrace ={\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}={E}_{N}.[{A}^{\text{c}}{]}^{T}.(\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace )\]

Après linéarisation de l’équation d’équilibre (4476), on introduit la matrice tangente \(\left[{K}^{\text{m},n-1}\right]\) qui contiendra les contributions issues de la linéarisation des efforts intérieurs et extérieurs et la matrice \(\left[{K}^{\text{c},n-1}\right]\) pourles forces de contact:

(2344)\[[{K}_{i}^{\text{c},n-1}]={E}_{N}.[{A}^{\text{c},n-1}{]}^{T}.[{A}^{\text{c},n-1}]\]

Le second membre valant:

(2345)\[\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace =-{E}_{N.}[{A}^{\text{c},n-1}{]}^{T.}\left\lbrace {d}^{\text{c},n-1}\right\rbrace\]

Concernant l’algorithme, nous avons vu au paragraphe r5.03.50-cas-regularisation qu’il agit en correction du problème d’équilibre sans contact.

Quand nous sommes dans l’algorithme de contact, à l’itération de Newton \(n\) , nous évaluons \(\left[{K}_{i}^{\text{c},n+1}\right]\) et \(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace\) , qui ne seront utilisés qu’à l’itération de Newton suivante. Il y a un décalage d’une itération: il faut nécessairementau moins deux itérations de Newton pour résoudre un problème avec contactpénalisé.

Algorithme

\(\mathit{Ini}\)	Calcul de \({⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{l}}\)
	Évaluation de \({\Xi}^{\text{c}}=\left\lbrace J\in {\Xi}^{l}\mid {⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{J}<0\right\rbrace\) et \(\left[{A}^{\text{c},n}\right]\)
	Pressions de contact \({⌊\left\lbrace {\mu}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}=-{E}_{N}.{⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\)
	Second membre \(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace ={\left[{A}^{\text{c},n}\right]}^{T}.{⌊\left\lbrace {\mu}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\)
	Matrice \(\left[{K}^{\text{c},n+1}\right]={E}_{N}.{\left[{A}^{\text{c},n}\right]}^{T}.\left[{A}^{\text{c},n}\right]\)

Remarques:

• Les quantités \(\left[{K}_{i}^{\text{c}}\right]\) et \(\left\lbrace {L}_{i}^{\text{c}}\right\rbrace\) ne sont évaluées que pour les liaisons actives. Pour le reste des inconnues, on initialise à zéro;

• La modification du système n’introduit pas de nouvelles variables par rapport au problème sans contact;

• Il n’y a pas d’indice \(k\) dans les différentes quantités (en particulier l’ensemble des liaisons actives ou la matrice de contact) car ce n’est pas un algorithme itératif;

• Pour le calcul de la matrice, on procède en deux temps: évaluations des matrices élémentaires, puis assemblage. Mais ces deux opérations sont spécifiques au contact discret, on n’utilise pas le mécanisme calculs élémentaires/assemblage générique de code_aster.

• Le choix du coefficient de pénalisation est crucial: trop faible, des interpénétrations seront observées, trop fort, le conditionnement de la matrice tangente se dégrade.

[1]

On note par \({⌊⌋}_{\square}\) les opérations qu’on fait sur un sous-ensemblede liaisons tel que défini dans le paragraphe r5.03.50-liaisons-frottement-sans.

Résolution algorithmique – Contact avec frottement

Liaisons de frottement

On ajoute deux ensembles par rapport à ceux définis dans le paragraphe r5.03.50-liaisons-frottement-sans pour prendre en compte le statut de frottement d’un nœud esclave:

• \({\Xi}^{\text{a}}\) est l’ensemble des nœuds de contact adhérents;

• \({\Xi}^{\text{g}}\) est l’ensemble des nœuds de contact glissants.

On a donc les relations suivantes entre les ensembles:

• \({\Xi}^{\text{c}}={\Xi}^{\text{a}}\oplus {\Xi}^{\text{g}}\) car les nœuds en contact sont soit glissants, soit adhérents;

• \({\Xi}^{\text{a}}\mathrm{\subseteq }{\Xi}^{\text{c}}\) et \({\Xi}^{\text{g}}\mathrm{\subseteq }{\Xi}^{\text{c}}\) car seuls les nœuds en contact peuvent être adhérents ou glissants.

Cinématique et jeu tangent

Nous avons vu dans le paragraphe r5.03.50-cinematique que le frottement utilise des quantités cinématiques analogues au cas du contact seul, en introduisant la notion de jeu «tangent» \(\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace\) (glissement relatif des surfaces dans le cas de non-adhésion), ce jeu tangent se définit par:

(2346)\[\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace =\left[{A}^{\text{f}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace )\]

On utilise la matrice de frottement \(\left[{A}^{\text{f}}\right]\) (voir r5.03.50-definition-matrice-frottement). L’utilisation de la quantité \(\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace\) à la place de \(\left\lbrace \delta {u}^{n}\right\rbrace\) se justifie du fait qu’on résout un problème de Tresca, donc sans modifier le jeu pendant le processus de résolution du problème de glissement (voir r5.03.50-modelisation-glissement).

Remarque:

• Dans les algorithmes présentés dans la suite du document, on ne fait généralement pas de distinction entre la partie glissante et la partie adhérente des nœuds au niveau des écritures, pour alléger les notations.

La méthode pénalisée en contact et en frottement

Équilibre de la structure en présence de contact

On rappelle que l’équation d’équilibre en présence de contact s’écrit:

(2347)\[\left\lbrace {L}_{i}^{int}({u}_{i})\right\rbrace =\left\lbrace {L}_{i}^{\text{ext}}({u}_{i})\right\rbrace -\left\lbrace {L}_{i}^{\text{c}}\right\rbrace\]

Dans le cas du contact pénalisé, la force de contact s’écrit:

(2348)\[\left\lbrace {L}_{i}^{\text{c}}\right\rbrace ={\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace }_{\text{regu}}={E}_{N}.[{A}^{\text{c}}{]}^{T}.\left(\left[{A}^{\text{c}}\right].\left\lbrace \delta {u}^{n}\right\rbrace -\left\lbrace {d}^{\text{c},n-1}\right\rbrace \right)\]

Après linéarisation de l’équation d’équilibre (2502) on introduit la matrice tangente \(\left[{K}^{\text{m},n-1}\right]\) qui contiendra les contributions issues de la linéarisation des efforts intérieurs et extérieurs et la matrice \(\left[{K}^{\text{c},n-1}\right]\) pourles forces de contact:

(2349)\[[{K}_{i}^{\text{c},n-1}]={E}_{N.}[{A}^{\text{c},n-1}{]}^{T.}[{A}^{\text{c},n-1}]\]

Le second membre valant:

(2350)\[\left\lbrace {L}_{i}^{\text{c},n-1}\right\rbrace =-{E}_{N.}[{A}^{\text{c},n-1}{]}^{T.}\left\lbrace {d}^{\text{c},n-1}\right\rbrace\]

Concernant l’algorithme, nous avons vu au paragraphe r5.03.50-cas-regularisation qu’il agit en correction du problème d’équilibre sans contact.

Quand nous sommes dans l’algorithme de contact, à l’itération de Newton \(n\) , nous évaluons \(\left[{K}_{i}^{\text{c},n+1}\right]\) et \(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace\) , qui ne seront utilisés qu’à l’itération de Newton suivante. Il y a un décalage d’une itération: il faut nécessairementau moins deux itérations de Newton pour résoudre un problème avec contactpénalisé.

Algorithme

Cette méthode est la plus simple à mettre en œuvre. Nous avons vu dans le paragraphe r5.03.50-linearisation-forces-contact que la régularisation des conditions de contact ajoute deux nouvelles contributions:

• \(\left[{K}_{i}^{\text{c}}\right]\) pour la matrice tangente globale;

• \(\left\lbrace {L}_{i}^{\text{c}}\right\rbrace\) pour le second membre.

Pour le frottement, il en est de même pour les conditions d’adhérence (cf. r5.03.50-linearisation-forces-adherence):

• \(\left[{K}_{i}^{\text{a}}\right]\) pour la matrice tangente globale;

• \(\left\lbrace {L}_{i}^{\text{a}}\right\rbrace\) pour le second membre.

Mais aussi pour les conditions de glissement (cf. r5.03.50-linearisation-forces-glissement ):

• \(\left[{\tilde{K}}_{\theta}^{\text{g}}\right]\) pour la matrice tangente globale;

• \(\left\lbrace {L}_{i}^{\text{g}}\right\rbrace\) pour le second membre.

Concernant l’algorithme, nous avons vu au paragraphe r5.03.50-cas-regularisation qu’il agit en correction du problème d’équilibre sans contact/ frottement . Quand nous sommes dans l’algorithme de contact/ frottement , à l’itération de Newton \(n\) , nous évaluons les quantités imposant les conditions de contact/frottement qui seront utilisés qu’à l’itération de Newton suivante. Ce qui nous donne l’algorithme suivant, avec \(n\) l’itération de Newton :

\(\mathit{Ini}\) Calcul de \({⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{l}}\) et \({⌊\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ⌋}_{{\Xi}^{l}}\) |

Évaluation de \({\Xi}^{\text{c}}=\left\lbrace J\in {\Xi}^{l}\mid {⌊\left\lbrace {\tilde{d}}^{\text{c},n}\right\rbrace ⌋}_{J}<0\right\rbrace\) |

Pressions de contact \({\lfloor \left\lbrace {\mu}^{\text{c},n}\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}=-{E}_{N.}{\lfloor \left\lbrace {\stackrel{~}{d}}^{\text{c},n}\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}\) |

Second membre \(\left\lbrace {L}_{i}^{\text{c},n}\right\rbrace ={\left[{A}^{\text{c},n}\right]}^{T.}{\lfloor \left\lbrace {\mu}^{\text{c},n}\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}\) |

Si \({\Xi}^{\text{c}}\ne \mathrm{\varnothing }\) |

Matrice de contact \(\left[{K}^{\text{c},n+1}\right]={E}_{N.}{\left[{A}^{\text{c},n}\right]}^{T.}\left[{A}^{\text{c},n}\right]\) |

Calcul de la norme du glissement \({\lfloor \Vert \left\lbrace {\stackrel{~}{g}}^{\text{t},n}\right\rbrace \Vert \rfloor }_{{\Xi}^{\text{c}}}\) |

Calcul de \({\lfloor \left\lbrace {\lambda}^{\text{f}}\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}=max\left(0,\mu .{\lfloor \left\lbrace {\mu}^{\text{c},n}\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}\right)\) |

Calcul de \({\lfloor \left\lbrace a\right\rbrace \rfloor }_{{\Xi}^{\text{c}}}={\lfloor \left(\Vert \left\lbrace {\stackrel{~}{g}}^{\text{t},n}\right\rbrace \Vert -\frac{1}{{E}_{T}}.\left\lbrace {\lambda}^{\text{f}}\right\rbrace \right)\rfloor }_{{\Xi}^{\text{c}}}\) |

Évaluation de \({\Xi}^{\text{a}}=\left\lbrace J\in {\Xi}^{\text{c}}|{\lfloor \left\lbrace a\right\rbrace \rfloor }_{J}\le 0\right\rbrace\) |

Évaluation de \({\Xi}^{\text{g}}=\left\lbrace J\in {\Xi}^{\text{c}}|{\lfloor \left\lbrace a\right\rbrace \rfloor }_{J}>0\right\rbrace\) |

Pour les nœuds adhérents \({⌊\left\lbrace {\mu}^{\text{f},n}\right\rbrace ⌋}_{J\in {\Xi}^{\text{a}}}=\sqrt{{E}_{T}}\) |

Pour les nœuds glissants \({⌊\left\lbrace {\mu}^{\text{f},n}\right\rbrace ⌋}_{J\in {\Xi}^{\text{g}}}=\sqrt{\frac{{⌊\left\lbrace {\lambda}^{\text{f}}\right\rbrace ⌋}_{J}}{{⌊∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥⌋}_{J}}}\) |

Calcul vecteur \({⌊\left\lbrace v\right\rbrace ⌋}_{{\Xi}^{\text{c}}}={⌊{\left[{A}^{\text{f},n}\right]}^{T}.\left\lbrace {\mu}^{\text{f},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\) et \({⌊\left\lbrace v\right\rbrace ⌋}_{{\Xi}^{\text{nc}}}=0\) |

Calcul première partie matrice de frottement \(\left[{K}^{{\text{f}}_{1}}\right]=\left\lbrace v\right\rbrace .\langle v\rangle\) |

Calcul du second membre de frottement \(\left\lbrace {L}_{i}^{\text{f},n}\right\rbrace =\left[{K}^{{\text{f}}_{1}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace )\) |

Calcul de la norme du glissement \({⌊∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥⌋}_{{\Xi}^{\text{c}}}\) |

Calcul de \({⌊\left\lbrace {\lambda}^{\text{f}}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}=max(0,\mu .{⌊\left\lbrace {\mu}^{\text{c},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}})\) |

Calcul de \({⌊\left\lbrace a\right\rbrace ⌋}_{{\Xi}^{\text{c}}}={⌊(∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥-\frac{1}{{E}_{T}}.\left\lbrace {\lambda}^{\text{f}}\right\rbrace )⌋}_{{\Xi}^{\text{c}}}\) |

Évaluation de \({\Xi}^{\text{a}}=\left\lbrace J\in {\Xi}^{\text{c}}\mid {⌊\left\lbrace a\right\rbrace ⌋}_{J}\le 0\right\rbrace\) |

Évaluation de \({\Xi}^{\text{g}}=\left\lbrace J\in {\Xi}^{\text{c}}\mid {⌊\left\lbrace a\right\rbrace ⌋}_{J}>0\right\rbrace\) |

Pour les nœuds adhérents \({\beta}_{J\in {\Xi}^{\text{a}}}=0\) |

Pour les nœuds glissants \({\beta}_{J\in {\Xi}^{\text{c}}}=\sqrt{\frac{1}{{⌊\left\lbrace {\lambda}^{\text{f}}\right\rbrace .∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥⌋}_{J\in {\Xi}^{\text{c}}}}}\) |

Choix de \(\theta\) |

Calcul vecteur \({⌊\left\lbrace w\right\rbrace ⌋}_{{\Xi}^{\text{c}}}=\sqrt{\theta}.{⌊{\beta}_{J\in {\Xi}^{\text{c}}}.\left[{A}^{\text{f},n}\right].\left\lbrace {L}_{i}^{\text{f},n}\right\rbrace ⌋}_{{\Xi}^{\text{c}}}\) et \({⌊\left\lbrace w\right\rbrace ⌋}_{{\Xi}^{\text{nc}}}=0\) |

Calcul seconde partie matrice de frottement \(\left[{K}^{{\text{f}}_{2}}\right]=\left\lbrace w\right\rbrace .\langle w\rangle\) |

Matrice résultante \(\left[{K}^{\text{f}}\right]=\left[{K}^{{\text{f}}_{1}}\right]-\left[{K}^{{\text{f}}_{2}}\right]\) |

On utilise pas mal d’astuces de calcul dans cet algorithme. On donne une valeur différente de \(\left\lbrace {\mu}^{\text{f}}\right\rbrace\) pour l’adhérence et le glissement mais on utilise le même vecteur \(\left\lbrace v\right\rbrace\) , ce qui permet de retrouver les deux matrices (glissement et adhérence). En effet, pour l’adhérence, on a:

(2351)\[{⌊\left[{K}^{{\text{f}}_{1}}\right]⌋}_{{\Xi}^{\text{a}}}=\left\lbrace v\right\rbrace .\langle v\rangle ={\left[{A}^{\text{f}}\right]}^{T}.\sqrt{{E}_{T}}.\sqrt{{E}_{T}}.\left[{A}^{\text{f}}\right]\]

On retrouve donc l’expression de la matrice d’adhérence ( ) pour les nœuds adhérents:

(2352)\[{⌊\left[{K}^{{\text{f}}_{1}}\right]⌋}_{{\Xi}^{\text{a}}}=[{K}_{i}^{\text{a}}]={E}_{T}.[{A}^{\text{a}}{]}^{T}.[{A}^{\text{a}}]\]

Pour le glissement:

(2353)\[{⌊\left[{K}^{{\text{f}}_{1}}\right]⌋}_{{\Xi}^{\text{g}}}=\left\lbrace v\right\rbrace .\langle v\rangle ={\left[{A}^{\text{f}}\right]}^{T}.\sqrt{\frac{\left\lbrace {\lambda}^{\text{f}}\right\rbrace }{∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥}}.\sqrt{\frac{\langle {\lambda}^{\text{f}}\rangle }{∥\left\lbrace {\tilde{g}}^{\text{t},n}\right\rbrace ∥}}.\left[{A}^{\text{f}}\right]\]

On retrouve l’expression de la première partie de la matrice de glissement ( ) pour les nœuds glissants :

(2354)\[{⌊\left[{K}^{{\text{f}}_{1}}\right]⌋}_{{\Xi}^{\text{g}}}=\frac{{\left[{A}^{\text{g}}\right]}^{T}.\left[{k}^{\text{g},n}\right]}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}\]

Avec la matrice des seuils de Tresca:

(2355)\[\left[{k}^{\text{g}}\right]=\left\lbrace {\lambda}^{\text{f}}\right\rbrace .\langle {\lambda}^{\text{f}}\rangle\]

Pour retrouver la second partie de la matrice de glissement, on utilise l’expression d u vecteur \(\left\lbrace w\right\rbrace\) :

(2356)\[\left\lbrace w\right\rbrace =\sqrt{\theta}.\beta .\left[{A}^{\text{f}}\right].\left[{K}^{{\text{f}}_{1}}\right].(\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace +\left\lbrace \delta {\tilde{u}}^{n}\right\rbrace )\]

Avec l’expression de la matrice \({⌊\left[{K}^{{\text{f}}_{1}}\right]⌋}_{{\Xi}^{\text{g}}}\) dans ( ) et la définition du glissement tangentiel dans ( ), le produit (tensoriel) de \(\left\lbrace w\right\rbrace\) par lui même permet de retrouver la seconde partie de la matrice de glissement ( ) pour les nœuds glissants :

(2357)\[{⌊\left[{K}^{{\text{f}}_{2}}\right]⌋}_{{\Xi}^{\text{g}}}=\left\lbrace w\right\rbrace .\langle w\rangle =\theta .\frac{{\left[{A}^{\text{g}}\right]}^{T}.\left[{k}^{\text{g},n}\right]}{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}.\frac{\left\lbrace {g}^{\text{t},n}\right\rbrace \langle {g}^{\text{t},n}\rangle }{{∥\left\lbrace {g}^{\text{t},n}\right\rbrace ∥}^{2}}\]

On a utilisé le produit ( ), ce qui explique le multiplicateur \(\beta\) utilisé dans le vecteur \(\left\lbrace w\right\rbrace\) .

Remarques:

• La modification du système n’introduit pas de nouvelles variables par rapport au problème sans contact/frottement;

• Il n’y a pas d’indice \(k\) dans les différentes quantités (en particulier l’ensemble des liaisons actives ou lesmatricesde contact/frottement) car ce n’est pas un algorithme itératif;

• Pour le choix de \(\theta\) , c’est l’utilisateur qui le choisit (paramètre COEF_MATR_FROT ), mais ce coefficient est mis à zéro tant que le résidu d’équilibre ( RESI_GLOB_RELA ) est inférieur à \({10}^{-3}\) ;

Résultats théoriques de convergence

Pour les problèmes sans frottement, on trouvera dans [bib1] une démonstration de la convergence pour la méthode des contraintes actives.

Pour les problèmes avec frottement, des résultats de convergence avec unicité de la solution pour le problème discrétisé sont établis dans [bib1] pour de faibles valeurs du coefficient de frottement de Coulomb. Les résultats sont obtenus en utilisant un algorithme de point fixe associé à une méthode de multiplicateurs de Lagrange. Pour chaque problème de contact résolu, on étudie le problème de frottement associé. Une fois celui-ci résolu, on résout un nouveau problème de contact et ainsi de suite. Ces méthodes sont cependant différentes de celles présentées ici et l’on ne peut donc pas présenter de résultats de convergence théorique pour ces dernières.

Redécoupage du pas de temps

Sur le plan théorique, la convergence de la méthode des contraintes actives est assurée en un nombre fini d’itérations. Dans la pratique, certains artefacts numériques peuvent rendre cette convergence délicate. Aussi une stratégie a-t-elle été mise au point pour assurer la robustesse de l’algorithme.

Lors de calculs de contact, notamment si les pas de charge effectués sont trop grands, des phénomènes indésirables peuvent apparaître:

• La matrice de contact est singulière,

• Oscillation de la méthode des contraintes actives: un nœud est détecté alternativement «collé» puis «décollé».

Pour pallier ces difficultés, la stratégie suivante a été adoptée. Si:

• La matrice de contact est singulière,

• Le nombre d’itérations de contraintes actives est supérieur à une limite qui dépend du nombre de liaisons potentielles. Ce nombre est fixé à deux fois le nombre total de nœuds esclaves pour la méthode des contraintes actives, et à ITER_CONT_MULT fois le nombre total de nœuds esclaves pour les autres méthodes.

Alors on redécoupe le pas de temps i.e. on revient au pas de charge précédent et au lieu d’essayer d’atteindre le niveau de chargement suivant en un pas comme on vient de le faire, on en fait plusieurs (Pour plus de précisions sur cette fonctionnalité de l’opérateur STAT_NON_LINE, voir la documentation [U4.51.03]).

Compatibilité avec les conditions aux limites de Dirichlet

Dans le cas des méthodes avec multiplicateurs de Lagrange, on peut observer des incompatibilités avec le fait d’imposer des conditions aux limites de type Dirichlet. En effet, il faut que physiquement le problème ait un sens. On ne peut pas traiter un problème de contact dans la direction de l’axe \(z\) si tous les points ont un déplacement nul suivant \(z\) . Comme nous allons le voir, traiter un tel problème conduit à une singularité des matrices du type \(\left[{A}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}^{\text{c}}\right]}^{T}\) avec le traitement des conditions aux limites de Dirichlet par double lagrange de code_aster .

Écriture des conditions aux limites

En s’inspirant de la documentation de référence [r5.03.01] de STAT_NON_LINE, la dualisation des conditions aux limites de Dirichlet conduit au système d’équations suivant à résoudre:

(2358)\[\begin{split}\lbrace \begin{array}{c}\left[K\right].\left\lbrace \delta u\right\rbrace +{\left[B\right]}^{T}.\left\lbrace \delta \lambda \right\rbrace =\left\lbrace {L}_{i}^{int}\right\rbrace -\left\lbrace {L}_{i}^{\text{ext}}\right\rbrace \\ \left[B\right].\left\lbrace {u}_{i}\right\rbrace =\left\lbrace {u}_{i}^{d}\right\rbrace -\left[B\right].\left\lbrace {u}_{i-1}\right\rbrace \end{array}\end{split}\]

On note alors \(\left[C\right]\) la matrice de rigidité du système telle que:

(2359)\[\begin{split}\left[C\right]=\left[\begin{array}{cc}\left[K\right]& {\left[B\right]}^{T}\\ \left[B\right]& \left[0\right]\end{array}\right]\end{split}\]

Cette matrice possède un inverse de la forme:

(2360)\[\begin{split}{\left[C\right]}^{-1}=\left[\begin{array}{cc}\left[E\right]& \left[F\right]\\ {\left[F\right]}^{T}& \left[G\right]\end{array}\right]\end{split}\]

tel que: \(\left[E\right].{\left[B\right]}^{T}=\left[0\right]\) . On vérifie ainsi que pour chaque condition aux limites \(l\) on a la propriété \(\left[E\right].{\left[{B}_{l}\right]}^{T}=\left[0\right]\) .

Retour au problème de contact

La matrice \(\left[{A}^{\text{c}}\right].{\left[K\right]}^{-1}.{\left[{A}^{\text{c}}\right]}^{T}\) peut aussi s’écrire \(\left[{A}^{\text{c}}\right].\left[E\right].{\left[{A}^{\text{c}}\right]}^{T}\) puisque les vecteurs de liaison \(\left[A\right]\) ne font intervenir que les degrés de liberté de déplacement.

• Il en résulte alors que si un vecteur de liaison \(J\) de la matrice \(\left[A\right]\) est une combinaison linéaire des conditions aux limites de type Dirichlet il vérifie la propriété suivante: \(\left[E\right].{\left[{A}_{J}\right]}^{T}=\left[0\right]\) . La matrice \(\left[{A}^{\text{c}}\right].\left[E\right].{\left[{A}^{\text{c}}\right]}^{T}\) est alors singulière car elle possède une colonne de zéros. Dans la pratique, sans traitement particulier, on finit dans le code sur un message d’arrêt du type arrêt sur matrice de contact-frottement singulière. La détection de ces colonnes singulières a été mise en œuvre dans le code de façon à éliminer des relations de contact-frottement ce type de relations et éviter l’arrêt précédemment décrit.

• Il en résulte alors que si un vecteur de liaison \(J\) de la matrice \(\left[A\right]\) contient une combinaison linéaire des conditions aux limites de type Dirichlet et s’écrit \(\left[{A}_{J}\right]=\sum{\alpha}_{i}.\left[{B}_{i}\right]+\left[{\stackrel{ˉ}{A}}_{J}\right]\) , il vérifie la propriété suivante: \(\left[E\right].{\left[{A}_{J}\right]}^{T}=\left[E\right].{\left[\stackrel{ˉ}{{A}_{J}}\right]}^{T}\) . On peut alors avoir une matrice \(\left[{A}^{\text{c}}\right].\left[E\right].{\left[{A}^{\text{c}}\right]}^{T}\) singulière car elle possède deux lignes identiques. Cette détection n’est pour le moment pas disponible dans le code et on finit dans le code sur un message d’arrêt du type arrêt sur matrice de contact‑frottement singulière.

Remarque:

Ce problème de compatibilité entre le contact-frottement et les conditions aux limites n’apparaît pas avec les méthodes régularisées dans la mesure où l’on rajoute de la rigidité à la rigidité globale et que l’on ne fait pas d’élimination comme dans le calcul des m ultiplicateurs de L agrange.

Conclusion

Des modélisations discrètes du contact-frottement avec surfaces de glissement 1D et 2D ont été implantées dans code_aster . Ces modélisations utilisables avec STAT_NON_LINE et DYNA_NON_LINE sont accessibles sous DEFI_CONTACT.

Les modélisations proposées s’appuient sur les maillages des surfaces venant en contact et permettent de retranscrire nœud à nœud les conditions de contact frottement entre les surfaces après discrétisation de la formulation variationnelle correspondante. La méthode s’étend alors sans difficulté des petits déplacements au cas des grands déplacements. En effet, l’absence d’utilisation d’éléments finis, entre les surfaces pouvant venir en contact, évite la grande distorsion de ces derniers, dans le cas de grands déplacements. On peut alors utiliser soit des conditions de liaisons nœuds à nœuds directes pour des maillages initialement compatibles, soit des conditions de liaisons nœuds à nœuds pondérées suivant une approche par projection de type maître-esclave pour des maillages incompatibles.

Dans le cas de surfaces de glissement 1D on a pu développer un algorithme n’utilisant que des multiplicateurs de Lagrange. La convergence finie de ce type d’algorithme est prouvée pour le contact unilatéral sans frottement et dans le cas avec frottement pour de faibles valeurs du coefficient de frottement de Coulomb Dans le cas de surfaces de glissement 2D, le contact frottant est traité par régularisation.

Bibliographie

[bib1] (1,2,3,4,5,6)

«La méthode des contraintes actives appliquées au contact unilatéral», Note interne EDF n° HI-75/93/016.

[bib2] (1,2,3)

«Algorithme de contraintes actives et contact unilatéral sans frottement», Dumont G., Revue européenne des éléments finis, Vol. 4 n°1/1995, pp. 55-73.

[bib3]

«Quelques méthodes pour traiter les problèmes de contact unilatéral en présence de grands glissements», Vautier I., Note interne EDF n° HI-75/97/013.

[bib4] «Évaluation des difficultés de modélisation du contact unilatéral pour des maillages 3D quadratiques», Vautier I., Compte-rendu interne EDF n° MMN/97/023.

[bib5] «Exemple d’utilisation des fonctionnalités de contact en grands déplacements dans le Code_Aster », Vautier I., Note interne EDF n° HI-75/97/034/A.

[bib6] «Un algorithme de Gradient Conjugué Projeté Préconditionné pour le traitement du contact unilatéral», Tardieu N., 8ème Colloque National en Calcul des Structures, Giens, 2007.

[bib7] «Différents algorithmes pour des problèmes de contact frottement 2D et 3D en grands déplacements», Ben Dhia H., Massin P., Tardieu N., Zarroug M., 4ème Colloque National en Calcul des Structures, Giens, 2001.

[bib8] «Les inéquations en mécanique et en physique», Duvaut G., Lions J.L., Dunod, Paris, 1972.

[bib9] «Analyse convexe et problèmes variationnels», Ekeland I., Temam R., Bordas, 1974.

[bib10] «Remarks on a numerical method for unilateral contact including friction», Licht C., Pratt E., Raous M., International Series of Numerical Mathematics, Vol. 101, 1991, pp. 129-144.

[bib11] «2D and 3D algorithms for frictional problems with small displacements», Massin P., Ben Dhia H., Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2000, Barcelona, 11-14 septembre 2000.

[bib12] «Finite Element Procedures for Contact-Impact Problems», Zhong Z., Oxford University Press, p.146-148, 1993.

[bib13]

«Computational Contact Mechanics», Wriggers P., Springer, Berlin, 2006

[bib14]

«Accélération de la convergence des méthodes de type Newton pour la résolution des systèmes non-linéaires», Ziani M., thèse de doctorat, 2008.

[bib15] (1,2)

Ronald W. Hoppe. Optimization theory, University of Houston, 2006.

[bib16] (1,2)

Elisabeth Wong. Active set methods for quadratic programming, PhD thesis, University of California, San Diego, 2011.