r5.03.52 Éléments de contact dérivés d’une formulation hybride continue#
Résumé:
Ce document décrit la manière dont des éléments de contact frottant sont dérivés d’une formulation hybride continue de problèmes de contact entre solides (2D ou 3D) en grandes transformations et précise la stratégie de résolution utilisée.
On utilise une formulation originale appelée formulation Lagrangienne stabilisée qui permet de retrouver les cas classiques de la littérature (lagrangien, lagrangien augmenté, pénalisation) par un choix judicieux de ses paramètres.
Cette formulation est disponible dans la commande DEFI_CONTACT sous le nom “CONTINUE”.
r5.03.52 Table des matières#
r5.03.52 Table des illustrations#
Cinématique#
Cinématique des solides#
On considère deux solides \({B}^{i}(i=1,2)\) déformables en contact frottant. Ces deux solides occupent dans leur configuration initiale l’adhérence de deux domaines \({O}_{0}^{1}\) et \({\Omega}_{0}^{2}\) de \({\mathrm{ℝ}}^{3}\) et, dans leur configuration courante à l’instant \(t\) , l’adhérence de \({O}_{t}^{1}\) et \({O}_{t}^{2}\) (également compris dans \({\mathrm{ℝ}}^{3}\) ), respectivement. On suppose que, dans leur configuration initiale, ces deux solides sont dans un état naturel, c’est-à-dire sans contraintes résiduelles ou pré-déformations. Considérons une particule matérielle \({M}^{i}\) du solide \({B}^{i}\) , repérée par le vecteur \({p}^{i}\) dans la configuration initiale. On désigne par \({\varphi}_{t}^{i}\) la transformation au temps \(t\) du solide \({B}^{i}\) , cette application transporte la particule \({M}^{i}\) vers la position \({x}_{t}^{i}\) dans la configuration courante. On a:
Au cours de leur mouvement, les solides peuvent entrer en contact, comme indiqué sur la . La frontière \({\Gamma}^{i}=\partial {\Omega}^{i}\) de chaque solide \({B}^{i}\) est décomposée en parties \({\Gamma}_{u}^{i}\) , \({\Gamma}_{g}^{i}\) et \({\Gamma}_{c}^{i}\) dans la configuration initiale, dont les intersections sont vides deux à deux:
Ces frontières se déforment dans la configuration courante en \({\gamma}_{u}^{i}\) , \({\gamma}_{g}^{i}\) et \({\gamma}_{c}^{i}\) . Le solide \({B}^{i}\) est encastré sur \({\Gamma}_{u}^{i}\) et soumis à une densité nominale de forces surfaciques notée \({g}^{i}\) sur la partie \({\Gamma}_{g}^{i}\) . Par ailleurs, on note \({f}^{i}\) le champ de densité volumique d’efforts appliqués sur les solides \({B}^{i}\) . Les parties des surfaces \({\Gamma}^{i}\) susceptibles d’entrer en contact lors de la déformation des deux solides sont notées \({\Gamma}_{c}^{i}\) . Ces surfaces ne sont pas connues a priori (c’est une des non-linéarités du problème à résoudre comme on le verra plus loin).
|
Figure 1 : Description des différentes configurations du problème mécanique |
On suppose l’existence de cartes régulières notées \({\Phi}^{i}\) décrivant les surfaces \({\Gamma}^{i}\) . Ces cartes sont définies comme suit:
où \(\omega\) est un domaine borné (de référence) contenu dans \({ℝ}^{2}\) , voir .
|
Figure 2 : Description des surfaces |
Les coordonnées des surfaces de contact \({\Gamma}_{c}^{i}\) sont transportées par application de la déformation, au temps \(t\) . Elles sont donc construites par composition des deux fonctions:
On définit le tenseur gradient de la transformation \(\underline{\underline{{F}_{t}^{i}}}\) par:
Si on désigne par \({u}^{i}\) le champ des déplacements du solide \({B}^{i}\) :
Alors le tenseur gradient de la transformation se ré-écrit:
En utilisant la fonction \({\varphi}_{t}^{i}\) , on peut écrire:
Où \(\otimes ` est le produit tensoriel et :math:`{\nabla}_{{p}^{i}}\) l’opérateur gradient (en espace). On introduit le jacobien de la transformation:
On notera \(\overline{{\Omega}^{i}}\) l’intérieur strict de \({\Omega}^{i}\) et \({\Gamma}^{i}\) sa frontière:
Construction des bases locales#
Les deux applications \({\Phi}^{i}\) et \({\varphi}_{t}^{i}\) permettent de définir le repère local sur le solide \(i\) . Ce repère local est établi par construction du plan tangent au point considéré et donc, par la donnée des deux vecteurs tangents \({t}_{\alpha =1,2}^{i}\) sur la configuration déformée de la surface \({\gamma}^{i}\) :
On peut définir de la même façon les vecteurs de base de la surface initiale \({\Gamma}_{c}^{i}\) :
En utilisant la définition du tenseur gradient de la transformation (), on exprime explicitement la relation entre \({T}_{\alpha}^{i}\) et \({t}_{\alpha}^{i}\) :
Il n’y a aucune raison que les \({t}_{\alpha}^{i}\) soient unitaires et que le repère \(({t}_{1}^{i},{t}_{2}^{i})\) soit orthogonal. Nous obtenons les normales par simple produit vectoriel:
Dans Code_Aster , la normale est toujours obtenue suivant l’algorithme ci-dessous:
Calcul des vecteurs tangents sur la configuration actuelle;
Modification des vecteurs tangents (suivant les options VECT_MAIT, VECT_ESCL);
Calcul du vecteur normal par produit vectoriel au point considéré.
Par construction, les vecteurs tangents \({t}_{\alpha}^{i}\) sont orthogonaux à la normale \({n}^{i}\) .
Régularisation des bases locales#
Il est parfois nécessaire d’utiliser un champ de normales «lissées» pour minimiser les risques d’instabilités du problème de contact lorsque la surface est décrite par des champs polynomiaux (typiquement le cas des éléments finis). Pour cela, une stratégie très simple a été implémentée dans le code (activable via l’option LISSAGE). L’algorithme est décrit ci-dessous:
Pour chaque élément \(e\) de la surface, pour chaque nœud \(s\) de cet élément partagé par \(k(s)\) élements, on calcule la normale \({n}^{e}(s)\) unitaire en ce nœud \(s\) ;
Pour chaque nœud \(s\) , on fait la moyenne \({n}_{L}(s)\) entre les normales de chaque élément partageant ce nœud commun:
On lisse la normale obtenue par une interpolation conforme au degré de la surface (linéaire ou quadratique). Si on désigne par \({N}_{e}^{s}\) les \(\mathit{nno}(e)\) fonctions de formes de l’élément \(e\) , alors la normale lissée \({n}_{\text{lis}}(\overline{\zeta})\) en un point \(\overline{\zeta}\) de l’élément s’écrit simplement:
Remarque: compte-tenu que l’on fait une moyenne arithmétique non-pondérée entre les normales, si deux éléments adjacents sont de taille et de normales très différentes, le lissage obtenu peut être de mauvaise qualité.
Quelques identités remarquables#
Dans ce paragraphe, nous allons présenter un certain nombre d’identités remarquables relatives aux propriétés de la base \(({t}_{1},{t}_{2},n)\) . Nous rappelons d’abord la définition des vecteurs tangents:
Les vecteurs \({t}_{\alpha}\) et \(n\) sontorthogonaux entre eux, donc leur produit scalaire est nul:
Le vecteur \(n\) est unitaire par définition, son produit par lui-même est une constante (qui vaut \(1\) ):
On définitla matrice de la première forme fondamentale ou tenseurmétrique \(\underline{\underline{m}}\) dont les composantes sont:
Si la base était orthonormale, ce tenseur se réduirait au tenseur identité.
L’inverse du tenseur métrique est aussi la première forme fondamentale de l’espace dual, ses composantes valent:
Pour passer des vecteurs tangents de la base covariante à la base contravariante, on utilise le tenseur métrique [1]_ :
Avec:
On définit aussi le tenseur de la seconde forme fondamentale, liée à la courbure, dont les composantes sont les suivantes:
Condition cinématique de non-interpénétration#
Pour traduire la non-interpénétration, on procède comme suit:
On couple les points des surfaces de contact deux à deux: c’est l’ appariement;
On impose entre les deux points appariés une condition de non-interpénétration suivant une direction donnée (généralement la normale à l’une des deux surfaces au point concerné).
Problème d’appariement#
Par convention, on appellera surface esclave la surface de contact définie sur le solide \(1\) et surface maître la surface appariée définie sur le solide \(2\) . On effectue l’appariement en recherchant, pour tout point \(x\) de la frontière \({\gamma}_{c}^{1}\) , le point \(\overline{x}\) de \({\gamma}_{c}^{2}\) le plus proche. Cela revient à résoudre le problème de minimisation suivant:
La solution \(\stackrel{ˉ}{\zeta}\) 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.
|
Figure 3 : Projection du nœud esclave sur la maille maître |
Les détails de l’implémentation de l’algorithme utilisé pour le problème d’appariement sont disponibles dans le document [R5.03.50] (la méthode d’appariement est strictement équivalente entre les méthodes discrètes et la méthode continue).
Définition du jeu#
L’appariement permet de définir le point \(\overline{p}\) de la surface \({\Gamma}_{c}^{2}\) , apparié à \(p\) à l’instant \(t\) dans la configuration d’origine, et le point \(\overline{x}\) apparié à \(x\) à l’instant t dans la configuration actuelle. Le vecteur distance \(d\) entre les points appariés est tout simplement:
C’est le résultat de la minimisation décrite dans le système ():
La première condition d’optimalité s’écrit:
Avec la définition des vecteurs tangents sur la surface déformée \({\gamma}_{c}^{2}\) :
On peut remonter aux vecteurs de base de la surface initiale \({\Gamma}_{c}^{2}\) :
En utilisant la définition du tenseur gradient de la transformation, on exprime explicitement la relation entre \({T}_{\alpha}\) et \({t}_{\alpha}\) :
Nous obtenons la normale par simple produit vectoriel:
La base construite par les trois vecteurs \(({t}_{1},{t}_{2},n)\) n’a aucune raison a priori d’être orthonormale. Par contre, par construction, la normale \(n\) est unitaire et normale au plan défini par les vecteurs \(({t}_{1},{t}_{2})\) qui définissent le plan tangent à la surface déformée \({\gamma}_{c}^{2}\) . En reprenant l’égalité () on a nécessairement:
\({d}_{n}\) est une quantité algébrique signée. Comme \(x={\Phi}_{t}^{1}(\zeta )\) , on obtient:
\({d}_{n}\) est appelé jeu entre les deux solides au point \(x\) . Ce qui nous donne la relation entre la base \(({t}_{1}^{1},{t}_{2}^{1},{n}^{1})\) définie sur le plan tangent au point esclave et la base \(({t}_{1},{t}_{2},n)\) construite par la projection du point \(x\) sur la surface déformée (maître) \({\gamma}_{c}^{2}\) :
Les deux repères sont coïncidents si les points sont en contact (si \({d}_{n}=0\) ). La distance se décompose en une partie normale \({d}_{n}\) et une partie tangentielle \({d}_{\tau}\) :
On obtient pour l’expression de \({d}_{\tau}\) :
Avec \({\underline{\underline{{P}_{}}}}^{\tau}\) la matrice de projection sur le plan tangent:
On peut écrire la condition de non-interpénétration entre \(x\in {\gamma}_{c}^{1}\) et \(\overline{x}\in {\gamma}_{c}^{2}\) suivant la direction de recherche \(n=-{n}^{2}(\overline{x})\) , normale entrante unitaire au point \(\overline{x}\) (dirigée vers l’intérieur de \({\Omega}_{t}^{2}\) ):
Cette quantité est objective.
Remarque:
Compte tenu du choix de la normale, le jeu est «inversé» par rapport au cas du contact en formulation discrète dans Code_Aster : il y a interpénétration si \({d}_{n}\ge 0\) ;
Vitesse relative#
Pour écrire la loi de frottement de Coulomb, nous aurons besoin de la notion de vitesse de glissement relative entre les deux solides. Nous reprenons pour cela la définition du jeu normal ():
On commence par dériver le terme de droite par rapport au temps:
Avec la convention habituelle de sommation sur les indices muets \(\alpha\) . Puis on dérive le terme de gauche:
On reconnaît la définition de \({t}_{\alpha}\) dans le terme \(\frac{\partial {\Phi}_{t}^{2}}{\partial {\zeta}_{\alpha}}\) (voir ()). On égalise les deux termes () et ():
En posant \(v(x,t)\) la vitesse du point \(x\) et \(\overline{v}(\overline{x},t)\) la vitesse du point \(\overline{x}\) apparié:
On obtient finalement en utilisant ():
On note \({v}_{\text{gliss}}\) la vitesse de glissement:
Cette quantité est objective car elle est invariante par changement de repère. La vitesse de glissement relative s’identifie à la différence de vitesse des deux solides, seulement lorsqu’ils sont en contact (\({d}_{n}=0\) ).
Opérateur de discontinuité#
On note \(\mathrm{〚}g\mathrm{〛}(p)\) la discontinuité de la grandeur \(g\) entre le solide \(1\) et le solide \(2\) . Cette discontinuité s’écrit au point \(p\) sur \(1\) qui se projette en \(\stackrel{ˉ}{p}\) sur \(2\) suivant l’application d’appariement:
Par la suite on abandonnera la notation \(\mathrm{〚}g(p)\mathrm{〛}\) au profit de \(\mathrm{〚}g\mathrm{〛}\) en faisant l’hypothèse implicite que la discontinuité s’écrit sur le couple \((p,\overline{p})\) des points appariés. On peut décomposer cette discontinuité en une partie normale \({\mathrm{〚}g\mathrm{〛}}_{n}\) et une partie tangentielle \({\mathrm{〚}g\mathrm{〛}}_{\tau}\) :
Avec:
Où \(n\) est la normale à la surface maître, dirigée vers l’intérieur du solide \(2\) .
La partie tangentielle:
Avec les notations introduites, le jeu \({d}_{n}\) peut s’écrire:
Il convient de remarquer que les quantités projetées sur la normale sont nécessairement scalaires, d’où la notation en \({\mathrm{〚}g\mathrm{〛}}_{n}\) et non en \({\mathrm{〚}g\mathrm{〛}}_{n}\) .
Variations des quantités cinématiques#
Les surfaces de contact sont décrites par un jeu de quatre paramètres \(({\zeta}_{1}^{1},{\zeta}_{2}^{1},{\zeta}_{1}^{2},{\zeta}_{2}^{2})\) , dans l’espace de référence (paramétrique) des éléments finis. L’approche maître/esclave permet de réduire ce nombre en décrivant les variables cinématiques \(\zeta\) par rapport à un corps pris comme référence (en l’occurrence, le corps défini comme «maître»). Le choix adéquat de la surface de référence n’est pas immédiat et reste une hypothèse de modélisation à prendre en considération, surtout dans les cas d’un grand glissement relatif des deux surfaces. Il est nécessaire de distinguer les dérivées lagrangiennes des dérivées eulériennes. Ce paragraphe vise à évaluer les différentes variations. On définit la variation \(\delta f({x}_{0})\) d’une fonction \(f(x)\) à une variable en un point donné \({x}_{0}\) par la quantité:
Dans le cas d’une fonction à plusieurs variables, on distingue la variation totale (par rapport à toutes ses variables), notée \(\tilde{\delta}f\) et la variation partielle par rapport à une variable \(x\) , notée \({\delta}_{x}f\) . La variation totale est la somme de toutes les variations partielles. Notons \(\tilde{\delta}\) la variation totale d’une quantité cinématique. Si on note \(t\) la variable temps, on aura, dans le cas d’une quantité eulérienne \(\stackrel{ˉ}{a}\) :
C’est donc la somme d’une variation \({\delta}_{t}\) par rapport au temps \(t\) et de la variation \({\delta}_{\zeta}\) par rapport aux coordonnées paramétriques:
Si on considère maintenant une quantité lagrangienne \(a\) , sa variation par rapport aux coordonnées paramétriques \(\zeta\) est nulle:
Les quantités liées à la surface esclave sont des quantités lagrangiennes, celles liées à la surface maître sont des quantités eulériennes.
Les formules de Weingarten#
On utilise d’abord le fait que la variation d’une constante est nulle quelque soit le paramètre qui varie. Or d’après () le produit \((n\cdot n)\) est constant, donc:
En développant () on obtient:
Le produit est nul si, et seulement si, \(\delta n\) est orthogonal à \(n\) , et donc \(\delta n\) appartient au plan tangent:
On remarque que la variation de la normale appartient au plan tangent, en utilisant les constantes \({c}_{\beta}\) :
Nous allons maintenant établir des formules très pratiques dont on aura besoin plus tard. La première chose à faire est de donner la formule établissant la valeur des coefficients \({c}_{\beta}\) . On commence par établir le produit scalaire suivant:
Ensuite, on part de la propriété ():
On calcule sa première variation, la variation d’une constante est nulle quelque soit le paramètre qui varie:
On retrouve les composantes de la seconde forme fondamentale, au signe près. Par identification avec (), on a la relation suivante:
La variation de la normale vaut donc:
Par application de (),on a:
En utilisant (),on a:
C’est la formule de Weingarten.
Variation du jeu normal#
A partir de la définition du jeu normal ():
On applique la première variation totale:
Par définition \(\overline{x}\) est une quantité eulérienne, sa variation s’exprime donc par la formule ():
De même pour la variation de la normale:
Comme \(x\) est une quantité lagrangienne, \(\tilde{\delta}x={\delta}_{t}x\) . Finalement:
Si on fait le produit de cette identité avec la normale \(n\) :
On considère les propriétés (), () et (), c’est-à-dire:
Ce qui simplifie l’expression:
La première variation du jeu n’entraîne pas l’apparition de termes de transport (comme dans l’équation ). \({d}_{n}\) est donc une quantité purement lagrangienne.
Variation des coordonnées paramétriques#
On repart de la première variation de l’expression() (en tenant compte du fait que \({d}_{n}\) est une quantité lagrangienne):
Que l’on multiplie cette fois par \(\frac{\partial \overline{x}}{\partial {\zeta}_{\beta}}\) :
On retrouvele tenseurmétrique \(\underline{\underline{m}}\) dont les composantes covariantes valent:
Si on considèrelapropriété (),c’est-à-dire:
L’expression se simplifie:
On repart de l’expression () qu’on va différencier:
En développant:
C’est le tenseur de courbure () au signe près. Pour trouver \(\delta {\zeta}^{\alpha}\) , il faut résoudre le système suivant:
Avec:
Variation des tangentes#
Comme \({t}_{\alpha}\) est une quantité eulérienne, onapplique() sur ():
Seconde variation du jeu normal#
La seconde variation sera notée \(\Delta \delta\) afin de conserver la distinction entre dérivées lagrangienne et eulérienne (et la notation \(\stackrel{}{˜}\) ). Pour calculer la seconde variation du jeu normal, on ne part pas de l’expression () mais directement de ():
La première variation du terme de droite de () vaut:
On va maintenant calculer la seconde variation (totale) en appliquant l’opérateur \(\tilde{\Delta}\) sur (). En premier lieu, on utilise le fait que \({\delta}_{t}x\) est une quantité lagrangienne, on a donc:
Par contre, \({\delta}_{t}\overline{x}\) est une quantité eulérienne:
Seconde variation du terme convectif de ():
Donc, la seconde variation de () vaut:
La première variation du terme de gauche de () vaut:
On va maintenant calculer la seconde variation en appliquant l’opérateur \(\tilde{\Delta}\) sur (). Elle s’écrit:
En égalant les deux équations () et (), on obtient:
On va multiplier () par la normale \(n\) . On considère les propriétés () et (), c’est-à-dire:
Ce qui nous donne, après simplifications:
On sait que \({d}_{n}\) est une quantité purement lagrangienne (§ 2.3.2 ), donc:
On suppose que les variations \({\Delta}_{t}{\delta}_{t}x\) et \({\Delta}_{t}{\delta}_{t}\overline{x}\) sont nulles. De plus, par analogie avec (), on a :
Ce qui nous donne:
On utilise ():
Donc:
On reconnaît le tenseur de courbure ():
Il reste à préciser l’expression de \(\tilde{\Delta}\tilde{\delta}n\cdot n\) . On a:
La formule de Weingarten () nous donne:
La produit scalaire étant commutatif, en utilisant le tenseur métrique contravariant, on obtient:
Il faut développer les variations par application des expressions de ():
Une fois de plus, on retrouve le tenseur de courbure:
Finalement:
Lois de contact et de frottement#
Dans ce chapitre, nous allons écrire les différentes formes équivalentes aux lois de contact-frottement standard (loi de Signorini et loi de Coulomb). La démonstration d’équivalence entre les différentes formes (Signorini, Coulomb, lagrangien augmenté, lagrangien stabilisé) est disponible dans [].
Dans Code_Aster , ce sont les formes «lagrangien stabilisé» qui sont utilisées.
Loi de contact#
Loi de Signorini#
La loi de Signorini pour le contact unilatéral s’écrit:
où \({d}_{n}\) est le jeu algébrique défini par l’équation () et \({\lambda}_{n}\) est la pression de contact. Le graphe de la loi de Signorini est représenté sur la . La première condition retranscrit la condition de non-adhérence et la seconde condition traduit l’impénétrabilité.
|
Figure 4 : Graphe de la loi de contact unilatéral de Signorini |
Le problème de contact ainsi posé introduit une relation non-univoque (\({\lambda}_{n}\) n’est pas une fonction de \({d}_{n}\) ), semi-définie positive et non-différentiable en \({\lambda}_{n}={d}_{n}=0\) . La loi de contact de Signorini est équivalente à une loi de comportement. Elle est analogue au cas de la plasticité rigide parfaite. Il convient donc de noter que le jeu \({d}_{n}\) est équivalent à la mesure d’une déformation et qu’il doit donc être objectif.
Loi de contact – Analyse convexe#
En utilisant les outils de l’analyse convexe, on peut écrire la loi de Signorini sous la forme suivante:
\(\partial ` désignant le sous-différentiel et :math:`{I}_{K}\) la fonction indicatrice du convexe \(K\) :
La loi de contact de Signorini est constituée d’inégalités (recherche de solution dans dans un convexe donné), on utilise cette écriture lorsqu’on résout le problème de contact par des méthodes d’optimisation sous contraintes: c’est la méthodologie privilégiée des formulations discrètes dans Code_Aster . Dans la méthode continue, on utilise une forme de semi-régularisation en utilisant des opérateurs de projection dans un convexe donné.
Loi de contact – Projection#
La méthode du lagrangien augmenté appliquée au contact consiste à utiliser la projection sur le convexe \({\mathrm{ℝ}}^{-}\) d’une quantité appelée multiplicateur de contact augmenté \({g}_{n}\) , défini par:
où \({\rho}_{n}\) est un réel strictement positif, homogène à une force sur une unité de volume:
Le problème de Signorini \({P}^{C}\) se ré-écrit alors:
L’idée de cette formulation est de remplacer la partie non-différentiable énergétique du contact, par le carré de la projection de \({g}_{n}\) dans un convexe approprié, ce qui rend la fonctionnelle énergie du lagrangien augmentée différentiable. L’opérateur de projection vaut explicitement:
Le multiplicateur de contact augmenté \({g}_{n}\) est une droite dans l’espace \(({d}_{n},{\lambda}_{n})\) . Avec l’opérateur de projection \({\mathit{Proj}}_{{\mathrm{ℝ}}^{-}}(x)\) , on se ramène sur le graphe de Signorini. \({\mathit{Proj}}_{{\mathrm{ℝ}}^{-}}(x)\) est continue sur \({\mathrm{ℝ}}^{-}\) et presque partout différentiable.
Loi de contact – Lagrangien stabilisé#
L’idée de la formulation par champ de signe (appelée aussi formulation stabilisée ) est de dissocier le statut d’un point de contact de la valeur du lagrangien augmenté. On introduit le même multiplicateur \({g}_{n}\) que dans le cas du lagrangien augmenté (équation ). On définit le champ de signe («level-set») noté \({S}_{u}^{x}\) par l’expression suivante:
Cet opérateur est la fonction caractéristique de la variable \({g}_{n}\) sur le convexe \({\mathrm{ℝ}}^{-}\) . Il n’est pas différentiable. Le problème de Signorini \({P}^{C}\) se ré-écrit alors:
Le multiplicateur de contact augmenté \({g}_{n}\) est une droite. Le caractère non-régulier de la loi de contact est le caractère non-différentiable provenant du champ de signe. \({S}_{u}^{x}\) est presque partout continue (sauf en \(x=0\) ) et presque partout différentiable. Les deux opérateurs \({S}_{u}^{x}\) et \({\mathit{Proj}}_{{\mathrm{ℝ}}^{-}}(x)\) sont liés par une relation simple:
C’est cette forme qui est utilisée dans Code_Aster .
Loi de frottement#
Loi de Coulomb#
Pour représenter les phénomènes de frottement, on utilise la loi de Coulomb qui s’écrit comme suit:
où \(\mu\) est appelé coefficient de frottement de Coulomb (supposé constant et isotrope), et \({v}_{\tau}\) est la vitesse relative tangente , définie comme la projection dans le plan tangent à la surface au point de contact du saut en vitesse du solide \(1\) par rapport au solide \(2\) . Quand le contact est glissant, comme \({\lambda}_{n}\) est toujours négatif par hypothèse, on peut écrire:
La quantité \(\frac{{v}_{\tau}}{\mathrm{\parallel }{v}_{\tau}\mathrm{\parallel }}\) correspond à la direction de glissement, inconnue, mais dont on sait qu’elle a lieu dans le plan tangent constitué des deux vecteurs \(({t}_{1},{t}_{2})\) . La troisième loi, analogue à celle de l’équation (c) introduite pour le contact, impose la nullité des forces de frottement lorsqu’il n’y a pas contact. Le graphe de la loi de frottement de Coulomb est représenté sur la .
|
Figure 5 : Graphe de la loi de frottement de Coulomb |
Le frottement induit la notion de seuil. La relation introduite par le frottement de Coulomb est non-univoque et non-différentiable . L’expression de \({v}_{\tau}\) est donnée par:
Le saut de vitesse tangentielle \({v}_{\tau}\) n’est pas objectif (voir []) . En toute rigueur, il faudrait utiliser \({v}_{\mathrm{gliss}}\) (voir p. 13 ). Cependant, ces deux définitions sont identiques quand il y a contact (\({d}_{n}=0\) ). La loi de frottement de Coulomb s’écrit en vitesse et non en déplacement , contrairement à la loi de Signorini. Ceci constitue une des difficultés à considérer dans le cas d’une résolution quasi-statique du problème de contact-frottement. On ne peut pas remplacer directement la vitesse par le déplacement sans faire l’hypothèse que le chargement dans la zone de contact varie de manière monotone, ce qui, en pratique, est une hypothèse impossible à satisfaire pour un chargement quelconque de la structure.
Néanmoins, on peut écrire la loi de Signorini sous forme incrémentale, en considérant l” incrément de déplacement , cette manière de faire est correcte sur le plan du respect de la loi de Coulomb. Mais la difficulté est qu’en ne connaissant pas a priori la valeur de la quantité totale de la pression de contact, le problème de Coulomb doit se résoudre en deux temps.
Loi de frottement – Analyse convexe#
En utilisant les outils de l’analyse convexe, on peut écrire la loi de Coulomb sous la forme suivante:
Avec le convexe \(C({\lambda}_{n})\) qui représente le cône de Coulomb:
Loi de frottement – Projection#
La méthode du lagrangien augmenté appliquée au frottement est analogue au cas du contact, sauf que le convexe considéré change. Considérons le multiplicateur de frottement augmenté \({g}_{\tau}\) , défini par:
où \({\rho}_{\tau}\) est un réel strictement positif, homogène à une force sur une unité de volume:
Le problème de Coulomb \({P}^{F}\) se ré-écrit alors:
L’opérateur de projection n’est pas orthogonal:
Loi de frottement – Projection normalisée#
Dans le contact, nous avons utilisé la même définition du multiplicateur augmenté entre la formulation lagrangienne augmentée et la formulation lagrangienne stabilisée. Par contre pour le frottement, nous allons procéder autrement. En premier lieu, on introduit une quantité \(\Lambda\) appelée semi-multiplicateur de Lagrange :
Cette quantité normalisée appartient à la boule unité:
On l’appelle semi -multiplicateur de Lagrange car cette quantité a deux significations:
Dans le cas du frottement adhérent, \(\Lambda\) représente la densité d’effort d’adhérence. C’est donc formellement le multiplicateur de Lagrange associé à la condition d’adhérence (vitesse relative tangentielle nulle);
Dans le cas du glissement, \(\Lambda\) indique la direction de glissement (normalisée) et on a \(\Lambda =\frac{-{\lambda}_{\tau}}{\mathrm{\parallel }{\lambda}_{\tau}\mathrm{\parallel }}=\frac{{v}_{\tau}}{\mathrm{\parallel }{v}_{\tau}\mathrm{\parallel }}\) ; ce qui n’est pas la définition d’un multiplicateur de Lagrange au sens «classique» (voir p. 3.2.3 );
On remarque que \({\lambda}_{n}\) étant toujours négatif ou nul (convention de pression habituelle), \({\lambda}_{\tau}\) est bien de signe opposé à \(\Lambda\) , qui indique la direction de glissement. Ce changement de variable permet de passer de l’inconnue en effort tangent vers la direction de cet effort. Dès lors, on introduit le semi-multiplicateur augmenté de frottement \({h}_{\tau}\) (pour le distinguer du \({g}_{\tau}\) «classique», voir () tel que:
Il convient de noter que le coefficient d’augmentation \({\rho}_{t}\) a une unité différente du \({\rho}_{\tau}\) dans le cas du lagrangien augmenté. En effet, il s’écrit comme le rapport d’un temps sur une longueur:
Dans la formulation lagrangienne augmentée, le coefficient d’augmentation \({\rho}_{\tau}\) est de la même unité que le coefficient d’augmentation de contact (voir ). Le problème de Coulomb \({P}^{F}\) se ré-écrit :
Cet opérateur de projection est orthogonal:
On notera:
On va procéder au calcul de la variation par rapport à la variable \(x\) :
On notera la matrice de projection sur la boule unité:
Loi de frottement – Lagrangien stabilisé#
On définit le champ de signe («level-set») noté \({S}_{f}^{{h}_{\tau}}\) par l’expression suivante:
Cet opérateur est la fonction caractéristique de la variable \({h}_{\tau}\) sur le convexe \(B(0,1)\) . Il n’est pas différentiable. Le problème de Coulomb \({P}^{F}\) se ré-écrit alors:
La forme lagrangienne stabilisée est une formulation exacte du problème de frottement. C’est cette forme qui est utilisée dans Code_Aster .
Loi de contact-frottement#
Lagrangien augmenté#
Nous écrivons le lagrangien augmenté du problème de contact-frottement. Considérons le convexe \(K\) tel que:
Alors le lagrangien augmenté du problème de contact/frottement s’écrit:
La projection sur le convexe \(K\) vaut explicitement:
Le cas (a) correspond au contact adhérent, le cas (b) est le contact glissant et le cas (c) est le cas sans contact (ni frottement).
Lagrangien stabilisé#
Nous écrivons le lagrangien stabilisé du problème de contact/frottement:
Avec les deux champs de signe:
C’est cette forme qui est utilisée dans Code_Aster .
Différences entre Lagrangien augmenté et Lagrangien stabilisé#
Les deux écritures appartiennent à la catégorie des méthodes d’optimisation sans contraintes. Les principales différences viennent de l’usage d’un «champ de signe» (ou level-set), et, pour le frottement, d’un changement de variable réduisant le disque de Coulomb à un disque de rayon unité.
En dissociant l’écriture du lagrangien augmenté de celle du statut (contact/décollement et adhérent/glissement), ceci nous permet d’introduire les inconnues de statuts explicitement dans le cas du lagrangien stabilisé;
Dans le lagrangien stabilisé, on utilise un semi-multiplicateur de frottement, ce qui a pour principal effet d’avoir une unité différente entre le lagrangien augmenté de frottement \({\rho}_{t}\) et le lagrangien augmenté de contact \({\rho}_{n}\) ; C’est cette forme qui est utilise dans Code_Aster .
On trouvera la démonstration de l’équivalence entre ces deux écritures dans [].
Formulation hybride continue du contact/frottement#
Pour exprimer l’équilibre de la structure, en prenant en compte le contact-frottement, nous avons deux manières de procéder:
On exprime la forme classique du principe des travaux virtuels en considérant le contact comme un effort d’interaction de surface (on l’intègre dans l’expression du travail des efforts intérieurs); Pour compléter le système, on y ajoute l’expression faible des lois de contact-frottement;
On écrit directement l’énergie du système sous forme d’inégalité variationnelle mixte (sans contraindre l’espace des solutions par les conditions de contact/frottement). L’expression des critères d’optimalité par l’utilisation du lagrangien augmenté (problème de point-selle) permet de retrouver le système établi à partir du PTV;
Espaces#
L’espace de Hilbert \({H}^{1}\) , espace vectoriel complet, usuellement l’espace de Sobolev muni d’une norme appropriée (le produit scalaire), est l’espace des fonctions de carré mesurables \({L}^{2}\) et dont la dérivée au sens faible appartient aussi à \({L}^{2}\) . L’espace dual de \({H}^{1}\) sera noté \({H}^{-1}\) . On désigne par \({\mathit{CA}}^{i}\) l’espace contenant les fonctions cinématiquement admissibles tels que:
\({H}^{1/2}(\Gamma )\) est l’espace des traces de fonctions (appartenant à \({H}^{1}\) ) sur la frontière \(\Gamma\) et \({H}^{-1/2}(\Gamma )\) est son espace dual. Pour simplifier, on notera \(H={H}^{-1/2}(\Gamma )\) et \(H={({H}^{-1/2}(\Gamma ))}^{2}\) son homologue à deux dimensions (pour le frottement).
Principe des travaux virtuels#
Les équations pour l’équilibre du corps \({B}^{i}\) sont:
Où \({\underline{\underline{P}}}^{i}\) désigne le premier tenseur des contraintes de Piola-Kirchhoff (non-symétrique), \({u}^{i}\) le champ des déplacements et \({\lambda}^{i}\) la densité d’effort due aux interactions de contact frottant entre les deux solides.
Travail virtuel des forces internes#
Le principe des travaux virtuels (PTV) permet d’exprimer la forme bilinéaire des efforts internes:
Où \({\stackrel{̃}{u}}^{i}\) désigne le champ des déplacements virtuels. \({\underline{\underline{S}}}^{i}\) est le second tenseur des contraintes de Piola-Kirchhoff, symétrique et purement lagrangien, lié au premier tenseur \({\underline{\underline{\Pi}}}^{i}\) par la relation:
Travail virtuel des efforts extérieurs#
La densité des forces volumiques \({f}_{v}^{i}\) et surfaciques \({f}_{s}^{i}\) appliquées aux deux solides produisent un travail \({G}_{\text{ext}}^{i}\) :
Travail virtuel des efforts de contact-frottement#
\({G}_{\text{cf}}^{i}\) est le travail des efforts de contact-frottement sur le solide \({B}^{i}\) :
\({\lambda}^{i}\) est la densité des efforts de contact-frottement subis par le solide \({B}^{i}\) . En utilisant la procédure d’appariement, décrite au sous-paragraphe § 2.2.1 , le principe de l’action et de la réaction s’écrit localement, sur les configurations initiales sous la forme suivante:
On pose \({\Gamma}_{c}={\Gamma}_{c}^{1}\) et \({\lambda}^{1}=\lambda\) et on peut donc prendre:
Ce qui nous donne une nouvelle forme du terme de réaction du contact:
\({\tilde{u}}^{i}\) sont les champs de déplacements virtuels, cinématiquement admissibles sur chacun des solides. En outre, on a posé \({\Gamma}_{c}={\Gamma}_{c}^{1}\) et \({\lambda}^{1}=\lambda\) . On signale que la densité d’efforts \({\lambda}^{2}\) est prolongée par zéro aux points de \({\Gamma}_{c}^{2}\) sans vis-à-vis sur \({\Gamma}_{c}^{1}\) . On décompose la densité d’effort de contact \(\lambda\) en une partie normale \({\lambda}_{n}.n\) et une partie tangentielle \({\lambda}_{\tau}\) telles que:
En projetant sur la boule unité (voir () dans le § 3.2.4 ), on a:
La réaction est discontinue et prend deux valeurs (contact ou non). On utilise le champ de signe \({S}_{u}^{{g}_{n}}\) introduit au § 3.1.4 :
Le semi-multiplicateur de frottement est également discontinu et prend deux valeurs (glissant ou adhérent). En utilisant le champ de signe \({S}_{f}^{{h}_{\tau}}\) introduit au § 3.2.5 :
En réinjectant l’équation dans l’équation :
Nous allons ajouter ici le terme dit «de stabilisation», qui permettra de retrouver le comportement d’un lagrangien augmenté standard. L’idée est la suivante, en cas de contact avec adhérence (\({S}_{u}^{{g}_{n}}=1\) et \({S}_{f}^{{h}_{\tau}}=1\) ), la réaction de contact frottement () s’écrit:
On peut modifier l’écriture de la réaction de contact/frottement:
Cet ajout ne change rien dans le cas du contact adhérent, car, dans ce cas, \({d}_{n}=0\) et \({v}_{\tau}=0\) , ce qui permet de retrouver (). Par contre, dans les autres cas (pas de contact ou contact glissant), les paramètres \({\rho}_{n}\) et \({\rho}_{t}\) permettent d’éviter le zéro sur la diagonale de la matrice. Finalement, le terme de réaction de contact-frottement dans le PTV peut se décomposer en une composante pour le contact:
Et une composante pour le frottement:
Formulation faible de la loi de contact#
Nous introduisons la forme faible de la loi de contact. Il suffit pour cela de reprendre l’équation (a) en la multipliant par le champ-test du Lagrange de contact:
Formulation faible de la loi de frottement#
Nous introduisons la forme faible de la loi de frottement. Il suffit pour cela de reprendre l’équation (b) en la multipliant par le champ-test du Lagrange de frottement:
Application du principe des travaux virtuels#
Si on considère le problème quasi-statique, on cherche à résoudre l’équilibre sous forme faible par application du principe des travaux virtuels:
Inégalité variationnelle mixte#
Globalement, il s’agit de trouver la solution du problème suivant (inégalité variationnelle mixte):
\({W}_{int}^{i}\) et \({W}_{\text{ext}}^{i}\) correspondent à l’énergie interne et à l’énergie externe. Nous allons écrire l’expression de l’énergie de contact et de frottement. Soit \({l}_{n}\) la densité continue et différentiable de l’énergie de contact:
Et \({l}_{t}\) la densité continue et différentiable de l’énergie de frottement:
Par rapport à une formulation lagrangienne augmentée classique, on a procédé au changement de variable \({\lambda}_{\tau}=\mu {\lambda}_{n}\Lambda\) (voir § 3.2.4 ). Pour la densité d’énergie de contact:
On rappelle que le jeu normal \({d}_{n}\) s’écrit:
De même, le champ de signe dépend des champs de déplacements et de lagrangien de contact:
Pour la densité d’énergie de frottement:
Pour retrouver l’expression de l’équilibre au sens faible, il suffit d’appliquer la condition d’optimalité () sur () et (). Par variation \({\delta}_{u}\) sur l’énergie de contact \({W}_{\text{c}}\) , on obtient:
On retrouve l’expression du travail virtuel des efforts de contact (). De même, par variation \({\delta}_{{u}^{i}}\) sur l’énergie de contact \({W}_{\text{f}}\) , on obtient l’expression du travail virtuel des efforts de frottement ():
Pour retrouver l’expression du travail virtuel de la loi de contact (équation ), on calcule la variation de \({W}_{c}\) par le lagrangien de contact \({\lambda}_{n}\) et pour le travail virtuel de la loi de frottement (équation ), on exprime la variation de \({W}_{f}\) par le lagrangien de frottement normalisé \(\Lambda\) .
Formulation variationnelle pénalisée#
Pour écrire la formulation variationnelle pénalisée, on régularise les lois de contact et de frottement en écrivant:
Les deux paramètres \({\kappa}_{n}\) et \({\kappa}_{\tau}\) sont des paramètres de pénalisation strictement positifs. \({(.)}^{+}\) désigne la partie positive et \({\mathit{Proj}}_{B(0,1)}\) la projection sur la boule unité (voir définition au § 3.2.4 ). La forme faible de la réaction de contact s’écrit:
Pour la réaction de frottement:
Les champs de signes sont modifiés:
Dans la formulation pénalisée, les multiplicateurs de contact \({\lambda}_{n}\) et les semi-multiplicateurs de frottement \(\Lambda\) n’ont théoriquement plus de raison d’être mais on les conserve pour des raisons d’architecture informatique. On écrit donc une forme faible correspondant à la loi de contact/frottement : on a donc accès aux pressions de contact en formulation pénalisée:
Ainsi qu’à la densité d’effort d’adhérence grâce à la forme faible de la loi de Coulomb:
Dynamique – Écriture en vitesse/impulsion#
Le modèle de Signorini est inadapté aux problèmes de contact-impact, c’est-à-dire aux problèmes de contact en dynamique. Comme les champs de déplacements sont irréguliers (champ de signe), leur intégration temporelle par des schémas de différences finies classiques de type Newmark provoque des oscillations parasites (non-physiques) au moment des chocs. Ces oscillations sont d’autant plus importantes que l’ordre du schéma est élevé.
Loi de Signorini-Moreau pour le contact#
Moreau a ré-écrit le modèle de Signorini pour l’adapter au cas de la dynamique:
Avec la vitesse normale définie par:
Le modèle de Moreau se justifie en écrivant que si le modèle de Signorini est vrai à l’instant initial \({t}_{0}\) , alors le respect des conditions de Moreau revient à résoudre le problème de Signorini pour tout \(t>{t}_{0}\) . Pour \(t={t}_{0}\) , les lois de Signorini-Moreau sont strictement équivalentes à celles introduites dans le § 3.1.1 :
Toujours avec le multiplicateur de contact augmenté \({g}_{n}\) défini par:
Pour \(t>{t}_{0}\) , les lois de Signorini s’écrivent:
avec les deux champs de signes suivants:
Et le multiplicateur de contact augmenté en vitesse \({\dot{g}}_{n}\) , défini par:
On remarquera que la contrainte sur le jeu s’est déplacée en partie sur la vitesse normale \({\mathrm{〚}v\mathrm{〛}}_{n}\) . Il faut aussi noter le changement de nature du paramètre d’augmentation \({\rho}_{n}\) qui a désormais la dimension d’une force sur une vitesse (et non plus d’une force sur un déplacement).
Loi de Coulomb pour le frottement#
Le modèle de Coulomb pour le frottement est déjà naturellement écrit en vitesse:
Application du principe des travaux virtuels#
Pour le problème en dynamique, on ajoute le terme d’inertie \({G}_{\text{dyn}}^{i}\) à l’équation d’équilibre. Le système à résoudre est donc:
\({G}_{\text{dyn}}^{i}\) est le travail des efforts d’inertie:
Il convient d’ajouter au système des conditions initiales:
Pour compléter le système, il faut aussi décrire le schéma d’intégration temporelle, qui relie les déplacements aux vitesses:
En s’appuyant sur les travaux de Jean [], nous utiliserons un schéma thêta d’ordre un:
Si \(\theta =0\) , on retrouve un schéma purement explicite. Si \(\theta =1\) , on retrouve un schéma purement implicite de type Euler (on retrouve en effet l’approximation de la vitesse réalisée en quasi-statique):
Le schéma est stable si:
où \({\omega}_{m}\) est la pulsation maximale du système dynamique. En pratique, on recommande de choisir \(\theta\) entre \(0.5\) et \(1\) .
Travail virtuel des efforts de contact-frottement#
Nous allons maintenant écrire le travail virtuel des efforts de contact-frottement sous forme forte. On a une composante pour le contact:
Et une composante pour le frottement:
On ne voit pas le champ de signe des vitesses \({S}_{v}^{{\dot{g}}_{n}}\) dans l’expression de la forme forte de la réaction due au frottement, mais il est présent implicitement dans l’expression du multiplicateur de contact \({\lambda}_{n}\) .
Formulation faible de la loi de contact#
Nous écrivons maintenant la forme faible de la loi de contact établie par l’équation ():
Formulation faible de la loi de frottement#
Nous écrivons maintenant la forme faible de la loi de frottement établie par l’équation ():
Conservation de l’énergie#
La formulation lagrangienne stabilisée pour le contact conserve l’énergie totale du système, la quantité de mouvement linéaire et la quantité de mouvement angulaire. Dans le cas du frottement, l’énergie est conservée si le frottement est adhérent (\({S}_{u}^{{d}_{n}}{S}_{v}^{{\dot{g}}_{n}}=1\) et \({S}_{f}^{{h}_{\tau}}=1\) ), par contre elle décroît (dissipation) dans le cas du frottement glissant (\({S}_{u}^{{d}_{n}}{S}_{v}^{{\dot{g}}_{n}}=1\) et \({S}_{f}^{{h}_{\tau}}=0\) ).
Algorithmes de résolution du système non-linéaire#
Formulation incrémentale#
Nous allons transformer les formulations variationnelles mixtes continues en une forme incrémentale, car le problème est toujours résolu sous cette forme, que ce soit en quasi-statique ou en dynamique. On se donne une discrétisation de l’intervalle de temps \(I=\left[0,T\right]\) à calculer:
À une étape \(k\) de chargement donnée, correspondant à l’instant fictif \({t}_{k}\) , pour \(k\in \left[1,{n}_{T}\right]\) , on suppose connus les champs à l’instant \({t}_{k-1}\) notés \({W}_{k-1}=({u}_{k-1}^{1},{u}_{k-1}^{2},{\lambda}_{n,k-1},{\Lambda}_{k-1})\) . On cherche les nouveaux champs solutions à l’instant \({t}_{k}\) notés \({W}_{k}=({u}_{k}^{1},{u}_{k}^{2},{\lambda}_{n,k},{\Lambda}_{k})\) . La formulation hybride s’écrit alors sous la forme générique suivante:
Trouver les champs \({W}_{k}\in {\text{CA}}^{1}\times {\text{CA}}^{2}\times H\times H\) tels que:
\(\forall {\tilde{W}}_{k} \in {\text{CA}}^{1} \times {\text{CA}}^{2} \times H \times H\)
La formulation variationnelle donnée par le système () est fortement non-linéaire. En effet, outre les non-linéarités «classiques» dues aux grandes déformations ou aux comportements non-linéaires des matériaux, plusieurs niveaux de non-linéarité inhérente au contact-frottement peuvent être par ailleurs distingués:
Une non-linéarité géométrique due à la méconnaissance a priori des couples de points en contact entre les deux surfaces;
Une non-linéarité due à la méconnaissance de l’état de contact des points appariés;
Une non-linéarité due à la dépendance en \({\lambda}_{n}\) du seuil de frottement \({\lambda}_{s}=\mu {\lambda}_{n}\) .
L’ensembles’écrit sous la forme d’un système non-linéaire:
Avec \({W}_{k}\) qui représentent les inconnues du système et \(p(W)\) les paramètres du problème. L’idée générale est de résoudre ce système non-linéaire en fixant certains paramètres du contact-frottement puis de linéariser le système par une méthode de Newton:
Fixer les paramètres du calcul \({p}_{0}({W}_{0})\) ;
Résoudre le système par la méthode de Newton \(F(W,{p}_{0}({W}_{0}))=0\) ;
Boucler en 1 tant que \({p}_{0}\mathrm{\ne }p(W)\) ;
La méthode de Newton consiste à appliquer un développement de Taylor au premier ordre:
En fonction des paramètres qui seront fixés, la méthode de résolution se décline en plusieurs variantes:
La méthode des points fixes;
La méthode de Newton partielle;
La méthode de Newton généralisée;
Discrétisation temporelle pour la loi de Coulomb#
Il convient de remarquer que la loi de frottement de Coulomb introduit la notion de vitesse relative de glissement \({v}_{\tau}\) entre les deux solides. En quasi-statique, on utilise un schéma d’Euler implicite. Par contre, en dynamique, deux situations se présentent:
On utilise la formulation en vitesse/impulsion proposée dans p. 4.5 ;
On utilise un schéma de Newmark standard en déplacement;
Dans ce dernier cas, on n’utilisera pas la vitesse donnée par le schéma, mais la forme incrémentale de type Euler implicite, comme en statique. Cette vitesse s’écrit sous forme incrémentale:
où \({\mathrm{〚}\Delta u\mathrm{〛}}_{\tau}\) représente l’incrément de la projection du saut de déplacement \(u\) sur le plan tangent à la surface de contact au point de contact considéré. Cette écriture incrémentale en déplacement du problème de Coulomb est équivalente à la formulation originale en vitesse. En pratique, on intègre le \(\Delta t\) dans le paramètre \({\rho}_{\tau}\) , c’est-à-dire qu’on demande à l’utilisateur de les choisir en prenant en compte qu’on doit y intégrer le \(\Delta t\) «typique» du problème. On doit donc prendre garde à ce que le paramètre \({\rho}_{t}\) que l’on introduit via COEF_FROT vaut \({\stackrel{ˉ}{\rho}}_{t}=\frac{{\rho}_{t}}{\Delta t}\) et non \({\rho}_{t}\) .
Paramètres du calcul#
Les paramètres du calcul viennent de la manière de traiter les non-linéarités du problèmes de contact-frottement. Selon les hypothèses de résolution, il existe plusieurs algorithmes de résolution. Dans l’algorithme initial (voir [], [] et []), les trois non-linéarités du contact/frottement sont résolues par une série de points fixes imbriqués. L’algorithme utilise cinq boucles imbriquées:
- Boucle \(\text{<K>}\) sur les pas de temps
- Boucle \(\text{<J>}\) sur la géométrie
- Boucle \(\text{<L>}\) sur le seuil de frottement
- Boucle \(\text{<I>}\) sur le statut de contact
Boucle \(\text{<N>}\) sur Newton
Fin boucle \(\text{<N>}\)
Fin boucle \(\text{<I>}\)
Fin boucle \(\text{<L>}\)
Fin boucle \(\text{<J>}\)
Fin boucle \(\text{<K>}\)
Il ne reste au niveau \(\text{<N>}\) que les termes non-linéaires n’ayant pas pu être traités. Il n’y a donc que les forces internes et la non-linéarité liée au terme non-différentiable de projection sur la boule unité de la loi de Coulomb.
Il y a trois paramètres que l’on peut fixer ou pas par les mots clefs ALGO_RESO_* , correspondant aux boucles J, L et I, donc huit algorithmes. Ils n’ont pas pour vocation d’être tous utilisés. En pratique seul trois sont utiles.
Code |
Géométrie (boucle <J>) ALGO_RESO_GEOM |
Seuil de frottement (boucle <L>) ALGO_RESO_FROT |
Statuts de contact (boucle <I>) ALGO_RESO_CONT |
|
Newton généralisé |
GNM |
NEWTON |
NEWTON |
NEWTON |
Newton partiel en géométrie (Mode par défaut) |
PNM |
POINT_FIXE |
NEWTON |
NEWTON |
Points fixes |
FPM |
POINT_FIXE |
POINT_FIXE |
POINT_FIXE |
Newton partiel |
Mode non recommandé |
POINT_FIXE |
POINT_FIXE |
NEWTON |
Newton partiel |
Mode non recommandé |
NEWTON |
NEWTON |
POINT_FIXE |
Newton partiel |
Mode non recommandé |
NEWTON |
POINT_FIXE |
NEWTON |
Newton partiel |
Mode non recommandé |
POINT_FIXE |
NEWTON |
POINT_FIXE |
Newton partiel |
Mode non recommandé |
POINT_FIXE |
POINT_FIXE |
NEWTON |
Dans l’algorithme de Newton généralisé (voir []),aucun paramètre n’est fixé et l’intégralité des non linéarités est traité au sein de la même itération de Newton.
Boucle \(\text{<K>}\) sur les pas de temps Boucle \(\text{<N>}\) sur Newton Fin boucle \(\text{<N>}\) Fin boucle \(\text{<K>}\) |
En pratique, seules trois versions sont intéressantes:
Le Newton généralisé GNM est l’algorithme qui donne les résultats le plus rapidement, mais il manque parfois de robustesse (et provoque donc des échecs de convergence). De plus, il est très dépendant du choix des paramètres d’augmentation et nécessite un contrôle des cycles pour être performant (voir p. 8 );
L’algorithme FPM est le plus robuste, mais il est très couteux (beaucoup d’itérations de Newton). De plus des pathologies de cyclage peuvent apparaître dans certaines situations. Par contre, il est relativement insensible au choix des paramètres d’augmentation;
L’algorithme intermédiaire PNM propose de garder le paramètre géométrique dans la boucle de point fixe et de résoudre les autres non-linéarités par Newton. C’est un excellent compromis entre les deux algorithmes précédents.
Boucle sur la géométrie#
Le problème de la non-linéarité géométrique vient du fait que l’on ne sait pas a priori quels vont être les couples de points en contact. Cette méconnaissance entraîne également celle du repère local (normal et plan tangent) en chaque point potentiel de contact.
Cette non-linéarité est particulièrement prononcée pour les problèmes avec des surfaces de contact gauches (c’est-à-dire non planes) et/ou potentiellement soumises à de grandes transformations (grands glissements relatifs) . En pratique, une seule correction est souvent suffisante si les surfaces de contact subissent peu de variation géométrique (petits déplacements).
Il est important de comprendre ce qu’implique la réactualisation géométrique: c’est, à chaque itération, un problème tout à fait neuf de contact-frottement qu’on résout, et où il faut repartir de la détermination des couples de points appariés. C’est pourquoi la boucle de géométrie est la plus externe.
La non-linéarité d’appariement entraîne aussi la non-connaissance du repère local. Or le repère local courant (la normale \(n\) et les deux tangentes \({t}_{1}\) et \({t}_{2}\) en 3D, dans la configuration déformée) est présent à plusieurs reprises dans l’écriture des équations de Signorini-Coulomb. Si on prend le cas de la normale, cela signifie qu’elle s’écrit, au cours du déplacement \({u}^{1}\) :
La normale est donc non-linéaire en déplacement. Cette non-linéarité peut être dérivée dans les équations ou l’on peut faire l’hypothèse que la normale reste constante au cours de la boucle <J> sur la géométrie:
Dans notre cas, on fait l’hypothèse que la normale (et les tangentes) reste constante pendant la boucle sur la géométrie. Le nombre d’itérations de la boucle est déterminé automatiquement par le code mais on peut le fixer par NB_ITER_GEOM. Il convient de noter qu’il y a toujours au moins un appariement au début du pas de temps, même si NB_ITER_GEOM=0.
Boucle sur le seuil de frottement#
La non-linéarité de seuil (frottement de Coulomb) vient du fait que le seuil de frottement dépend de la densité normale de l’effort de contact \({\lambda}_{n}\) , qui est lui-même une inconnue du problème. La boucle de point-fixe <L> sur le seuil de frottement \({\lambda}_{s}\) revient à traiter le problème de contact-frottement avec un seuil de frottement fixe , ce qui revient à chaque itération à traiter un problème de Tresca . Cette stratégie est très courante pour résoudre le frottement de Coulomb, et se justifie théoriquement: en fixant le seuil de frottement, on peut dériver un potentiel convexe associé à la loi de Coulomb, qui permet de relier les inconnues de déplacement à leurs variables duales associées, les forces de frottement.
Boucle sur les statuts de contact#
La dernière non-linéarité concerne la méconnaissance de l’état effectif de contact. Nous avons adopté une méthode dérivée de la méthode des contraintes actives implémentée dans Code_Aster pour le contact discret (voir [R5.03.50]). Par défaut, le statut initial des nœuds de contact est non-contactant dans Code_Aster . Par ailleurs \({\lambda}_{s}\) est initialisé à zéro (sauf si l’on utilise SEUIL_INIT). Cela revient à commencer par résoudre un problème sans contact à la première itération, puis avec contact mais sans frottement et enfin à activer le frottement avec un seuil de frottement \({\lambda}_{s}\) fixe (problème de Tresca). On peut néanmoins partir d’un état initialement contactant en utilisant l’option CONTACT_INIT. Dans ce cas, il est possible de résoudre des problèmes hypostatiques avec la formulation continue. C’est un avantage important de la formulation continue du contact par rapport aux formulations discrètes qui ne supportent pas un problème mécanique qui ne soit «tenu» que par le contact.
De plus, la méthode «glissière» permet de simuler du contact bilatéral en supposant que les nœuds entrant en contact le resteront tout au long du transitoire (fonctionnalité activable par le mot-clef GLISSIERE).
Linéarisations des termes#
Le système algébrique non-linéaire () va se résoudre par linéarisation de type Newton. Nous reprenons les expressions de l’équilibre (PTV, voir p. 4.2 ) et de la loi de Signorini et de Coulomb sous forme faible. Les quantités virtuelles sont alors des variations infinitésimales, donc les quantités exprimées en terme de coordonnées actualisées (du genre \(\delta x\) ) se transforment en variation sur les déplacements \(\delta u\) :
Note importante: la linéarisation des termes dans les paragraphes suivants sous-entend l’utilisation d’un schéma de Newmark en dynamique ou le cas quasi-statique, qui utilise l’écriture incrémentale en déplacement de la loi de frottement.
Le PTV comprend en particulier la réaction de contact, qui va s’écrire avec la variation du jeu normal:
Pour la réaction de frottement dans le PTV:
Il est important de comprendre le sens des notations. Dans les deux cas précédents, \(\tilde{\delta}\) s’entend comme la variation totale par rapport aux déplacements (expressions établies à la p. 2.3 ).
Pour la loi de Signorini sous forme faible on considère:
Dans ce cas, \(\delta\) s’entend comme la variation par rapport à la pression de contact \({\lambda}_{n}\) .
Enfin pour la loi de Coulomb sous forme faible:
\(\delta\) s’entend comme la variation par rapport au Lagrange de frottement \(\Lambda\) .
En ajoutant les termes correspondant à l’équilibre de la structure (\({G}_{int}^{i}\) et \({G}_{\text{ext}}^{i}\) ), on a le système non-linéaire suivant à résoudre:
Avec \({W}_{k}\) qui représentent les inconnues du système et \(\delta {W}_{k}\) les quantités virtuelles. La première équation représente l’équilibre de la structure avec contact-frottement, comme la variation s’entend donc par rapport aux déplacements, on a:
L‘équation suivante est l’expression de la loi de contact de Signorini. La variation s’entend par rapport au lagrangien de contact:
Pour la loi de frottement de Coulomb. La variation s’entend par rapport au lagrangien de frottement:
Nous avons déjà fait une hypothèse dans ces deux dernières équations: les multiplicateurs de Lagrange du contact \({\lambda}_{n}\) et les semi-multiplicateurs du frottement \(\Lambda\) n’ont été écrits que sur un des deux solides. On a donc supposé que \({\Gamma}_{c}={\Gamma}_{c}^{1}\) ou \({\Gamma}_{c}={\Gamma}_{c}^{2}\) . La surface privilégiée est appelée surface esclave. On évaluera les intégrales sur cette surface. La linéarisation introduit en particulier des termes correspondant au jacobien généralisé du contact-frottement. Pour l’équilibre:
Pour la loi de contact de Signorini:
Pour la loi de frottement de Coulomb:
Au final, nous obtiendrons le système linéarisé suivant:
Nous devons évaluer tous les termes de la matrice jacobienne \(\left[J\right]\) .
Les variantes pour la linéarisation de la loi de Coulomb#
La loi de Coulomb relie la pression de contact à deux quantités:
L’intensité des contraintes tangentielles (forces de frottement);
La direction de frottement \(\tau\) ;
La dépendance de la direction de frottement par rapport à l’intensité de la pression de contact est un point important et qui montre la difficulté de prise en compte du modèle de Coulomb. Cette dépendance est primordiale en 3D, sur des surfaces courbes.
Si on reprend l’expression de l’équilibre, avec intégration de la réaction de frottement \({G}_{f}\) , on a donc, en toute rigueur, une dépendance à la pression de contact \({\lambda}_{n}\) de \(\mathrm{\parallel }{\lambda}_{\tau}\mathrm{\parallel }\) et de \(\tau\) , mais aussi, par l’expression de la loi faible de frottement, une dépendance de \({\tilde{G}}_{f}\) à la pression de contact \({\lambda}_{n}\) . Trois quantités sont donc variables:
C’est le cas le plus général, toutefois, on peut distinguer en toute rigueur 8 configurations possibles:
Variante |
\({\Delta}_{{\lambda}_{n}}(\mathrm{\parallel }{\lambda}_{\tau}\mathrm{\parallel })=0\) |
\({\Delta}_{{\lambda}_{n}}(\tau )\) |
\({\Delta}_{{\lambda}_{n}}({\tilde{G}}_{f})\) |
Matrice |
1 |
0 |
0 |
0 |
Symétrique |
2 |
0 |
0 |
\(\mathrm{\ne }0\) |
Non symétrique |
3 |
0 |
\(\mathrm{\ne }0\) |
0 |
Non symétrique |
4 |
0 |
\(\mathrm{\ne }0\) |
\(\mathrm{\ne }0\) |
Non symétrique |
5 |
\(\mathrm{\ne }0\) |
0 |
0 |
Non symétrique |
6 |
\(\mathrm{\ne }0\) |
0 |
\(\mathrm{\ne }0\) |
Non symétrique |
7 |
\(\mathrm{\ne }0\) |
\(\mathrm{\ne }0\) |
0 |
Non symétrique |
8 |
\(\mathrm{\ne }0\) |
\(\mathrm{\ne }0\) |
\(\mathrm{\ne }0\) |
Non symétrique |
Toutes ces variantes sont numériquement acceptables (ce sont des variantes de Newton) mais seules quelques unes ont un vrai sens physique. La variante 1 correspond au cas où la pression de contact est constante: c’est le modèle de Tresca ou la variante dite «point fixe sur le seuil de frottement», c’est la seule qui produit une matrice tangente symétrique. La variante 8 est le cas le plus général (méthode de Newton généralisée). À part ces deux cas, seuls deux autres variantes ont un sens physique:
La variante 4 suppose qu’on néglige la dépendance de l’intensité de frottement par rapport à la pression de contact mais qu’on peut avoir des variations de direction des directions tangentielles «brusques». Le cas typique est celui d’une surface courbe en 3D avec de faibles intensités de pression de contact (contact rasant);
La variante 6 suppose qu’on peut négliger la dépendance de la direction de frottement par rapport à la pression de contact (cas 2D ou 3D plan) mais pas la dépendance de la force de frottement. Le cas typique est celui d’une surface plane 3D ou 2D, avec une forte intensité de pression de contact;
En pratique, la variante 8 introduit des termes numériquement mal conditionnés qui ralentissent la convergence. La dépendance à la direction de frottement couple des termes géométriques avec des termes en pression. Nous avons donc introduit la variante 5 qui semble être un bon compromis. Le calcul des termes de la jacobienne est reporté en annexe (voir p. 63 )
Dans les hypothèses de simplification, on a négligé la deuxième variation du jeu tangent (équation ) car les équations qui en découlent sont très complexes et ne semblent pas apporter un gain significatif sur la vitesse de convergence. Bien entendu, on perd ainsi la convergence quadratique théorique de l’algorithme de Newton.
Discrétisation de la formulation variationnelle mixte#
On décrit dans cette section la discrétisation du problème défini par les formulations variationnelles données p. 4 . Les formulations variationnelles utilisées ont les points communs suivants:
L’écriture de PTV pour l’équilibre de la structure, en utilisant une formulation faible de l’équation d’équilibre;
Un champ vectoriel d’inconnues: celui des variables primales représentant soit les déplacements, soit les vitesses (en dynamique, dans le cas du thêta-schéma);
La géométrie est représentée par une triangulation régulière (éléments finis polynomiaux classiques);
Un champ scalaire inconnu représentant le multiplicateur de Lagrange de contact \(\lambda\) ;
Un champ vectoriel inconnu représentant le semi-multiplicateur de Lagrange de frottement \(\Lambda\) ;
Le modèle de contact de Signorini est appliqué de manière faible, avec une formulation Lagrangienne augmentée;
Le modèle de frottement de Coulomb est appliqué de manière faible, avec une formulation Lagrangienne augmentée;
Des champs de signe, discontinus, pour le contact et le frottement ont été introduits dans la formulation.
Dans ce chapitre, nous allons introduire tous les ingrédients nécessaires à la fabrication du modèle discrétisé, puis, dans le chapitre suivant, nous verrons la manière de résoudre le problème.
Important: désormais, nous noterons \(\lambda\) et non plus \({\lambda}_{n}\) la pression de contact pour alléger les notations. La pression de frottement est toujours notée \(\Lambda\) .
Éléments théoriques#
Pavage du domaine#
Les domaines \({\Omega}^{i}\) sont approchés par des domaines \({\Omega}_{h}^{i}\) polygonaux. À chaque domaine \({\Omega}^{i}\) est associé une famille de triangulations \({T}_{{h}^{i}}\) avec le paramètre \({h}^{i}\) qui est le pas du maillage du domaine \({\Omega}^{i}\) . En notant \({K}_{j}^{i}\) les \({\mathit{NE}}^{i}\) éléments de chaque triangulation, nous pouvons écrire:
Pour la discrétisation de la frontière de contact \({\Gamma}_{c}^{i}\) , nous sommes face à une première hypothèse importante. En effet, cette frontière peut être discrétisée de deux manières:
On prend la trace géométrique du maillage du domaine \({\Omega}^{i}\) sous-jacent;
On crée un nouveau maillage indépendant, avec son propre pas de discrétisation.
En pratique, le maillage de la frontière est quasiment toujours construit sur la première hypothèse, car cela impliquerait sinon un travail de maillage difficile à mettre en œuvre, particulièrement en 3D. On utilisera donc toujours la trace géométrique du maillage du domaine \({\Omega}^{i}\) comme maillage pour la frontière de contact. La trace de \({T}_{{h}^{i}}\) forme un maillage \({T}_{{h}^{i}}^{\Gamma}\) de \({\Gamma}_{c}^{i}\) . La frontière discrétisée \({\Gamma}_{c,h}^{i}\) est composée de \({\mathit{NE}}_{c}^{i}\) éléments \({k}_{j}^{i}\) :
Espaces d’approximations discrets#
La version discrète de l’espace \({H}^{1}({\Omega}^{i})\) , sur le domaine \({\Omega}^{i}\) ,que l’on notera \({H}_{h}^{i,q}\) , est telle que:
avec \({\mathit{NE}}^{i}\) le nombre d’éléments provenant de la triangulation (), \(d\) la dimension du problème, \(C({\stackrel{ˉ}{\Omega}}^{i})\) l’espace des fonctions continues sur l’adhérence de \({\Omega}^{i}\) et \({P}_{q}({K}_{j}^{i})\) l’espace de polynômes de degré inférieur ou égal à \(q\) sur le domaine. On a désigné par \({\mathit{CA}}^{i}\) l’espace contenant les fonctions cinématiquement admissibles. La version discrète de cet espace, noté \({\mathit{CA}}_{h}^{i,q}\) est construite comme ():
Nous allons écrire la version discrète des l’espace des traces de fonctions sur la frontière:
Désormais, on oubliera la notation \(q\) pour indiquer l’ordre des polynômes.
La condition Inf-Sup ou condition LBB#
Le caractère mixte de la formulation impose une dépendance entre les espaces éléments finis des multiplicateurs de contact et de frottement et l’espace éléments finis des déplacements. Ce choix relève de la condition de compatibilité (ou condition Inf-Sup ou condition LBB de Ladyzenskaya-Babuska-Brezzi). Un problème mécanique se ramène en général à résoudre une équation variationnelle du type:
où \(U\) est un espace de Hilbert. Il est souvent plus commode de réécrire ce problème sous une forme dite mixte , à deux inconnues \(u\in U\) et \(\lambda \in M\) :
où \(M\) est aussi un espace de Hilbert, dual de \(U\) par la forme bilinéaire \(b\) . Le problème est que la formulation donnée par () n’est pas équivalente à la formulation donnée par (). Ainsi, dans ce cadre très général, l’ équivalence des deux formulations, ainsi qu’une condition nécessaire et suffisante d’existence et d’unicité (en plus, bien entendu, des hypothèses classiques sur l’ ellipticité et la coercivité de la forme bilinéaire \(a\) ) de la solution de (), est donnée par la condition dite de Ladyzenskaya , Babuska et Brezzi (LBB):
Cette condition indique comment il faut choisir l’espace \(M\) , pour que les formulations soient équivalentes. Cependant, il faut avouer que d’un point de vue opératoire, on n’est pas tellement plus avancé, car cette condition, très abstraite, est aussi très technique à mettre en œuvre. D’autant que cette condition est exprimée dans le cadre continu et que son passage dans le cadre discret n’est pas garanti! On se contentera de donner des règles d’application issues de la littérature.
L’application de la condition LBB à notre problème de contact dépend de la nature des interfaces:
Si l’interface est gauche, l’obtention d’une condition Inf-Sup avec une discrétisation discontinue des multiplicateurs de contact nécessite une discrétisation très fine de ces multiplicateurs;
Si l’interface est plane, l’espace \({\text{CA}}_{h}^{\gamma ,i}\) des traces des champs de déplacements sur l’interface de contact comme espace de discrétisation du champ de multiplicateur de contact remplit la condition LBB;
La méthode continue utilise l’espace des traces des champs de déplacements, sachant que la stabilité (au sens LBB) ne sera pas garantie pour les interfaces non planes.
On prend des champs de multiplicateurs dans l’espace \({\text{CA}}_{h}^{\gamma ,i}\) sur toute l’interface de contact [bib2] . Pour des surfaces gauches (), ce choix est commode (puisque nous discrétisons le champ de multiplicateurs de manière identique au champ de déplacement), mais il n’y a aucune preuve de stabilité . En effet, le problème vient de ce que le choix d’une des deux surfaces comme surface d’intégration constitue une approximation plus ou moins grossière, comme l’illustre la . L’expérience montre cependant que cette difficulté peut être contournée par l’utilisation de champs d’approximation aussi lisses que possible (lissage des normales, continuité des champs de multiplicateurs). Nous verrons plus tard que le respect de la condition Inf-Sup est aussi fortement lié au schéma de quadrature numérique utilisé pour évaluer les intégrales.
|
Figure 6 : Surfaces gauches incompatibles |
Pour l’utilisateur du contact dans Code_Aster , ce choix signifie qu’il risque d’avoir des difficultés de convergence dès lors que les interfaces de contact sont gauches, que le glissement relatif est significatif ou que la différence de discrétisation est importante entre les deux surfaces.
Champs de signe#
Le choix des espaces d’approximations (de type éléments finis polynomiaux) pour le champ de déplacement, le champ de multiplicateur de Lagrange de contact et le champ des semi-multiplicateurs de frottement n’est pas suffisant pour construire le problème discret. Il reste le cas des champs de signes (pour le contact, le frottement, le champ de vitesse, etc) qui sont naturellement dans des espaces infinis. On utilise une méthode de collocation pour approximer ces champs de signe. On choisit donc une collection finie de points de la surface de contact discrétisée \({\Gamma}_{c,h}^{i}\) :
Un champ de signe \(S\) sera donc approché par \({S}_{h}\) , qui est la somme discrète suivante:
Attention! Le nombre de points de collocation \({n}_{c}^{i}\) n’est pas nécessairement le même que le nombre de points issus de la discrétisation éléments finis des multiplicateurs de contact/frottement (qui repose, on le rappelle, sur la trace géométrique du volume discrétisé).
Intégration numérique#
Pour évaluer numériquement les intégrales, il convient de choisir un schéma de quadrature numérique idoine. On se fonde pour ça sur l’intégration des termes discontinus (champs de signe) qui «guident» le meilleur schéma à adopter. La question se déplace donc vers le choix des points de collocation optimaux pour obtenir des efforts de contact précis.
Le respect de la condition Inf-Sup a été abordé dans le cadre du choix des fonctions d’approximation des champs de multiplicateurs de contact et de frottement. Toujours dans l’optique de respecter cette condition, il a été dit qu’elle était fortement liée à la précision de l’intégration numérique qu’on est capable d’obtenir. Le phénomène se comprend quand on fait l’analogie avec les méthodes de sous-intégration numérique parfois employées pour résoudre des problèmes de blocages numériques dus à une trop grande richesse de l’approximation, en faisant alors «baisser» le degré d’interpolation des champs en sous-intégrant. Le raisonnement peut donc être renversé facilement: si l’intégration numérique «fausse» (en enrichissant ou en appauvrissant) l’interpolation polynomiale des champs de multiplicateurs, on risque de sortir du domaine de stabilité de la condition LBB.
|
Figure 7 : Intégration numérique et compatibilité |
C’est le cas, notamment en cas d’incompatibilité. Imaginons que l’aire élémentaire d’intégration sur la surface esclave recouvre plusieurs éléments de la surface maîtresse. Comme les fonctions-test ont pour support non-nul un élément fini, le fait d’intégrer numériquement, par une quadrature de Gauss, une fonction non polynomiale, mais identiquement nulle sur une partie de l’aire d’intégration, n’est clairement pas exact, car la quadrature est tout à fait inappropriée (voir ). Or c’est ce qu’on fait quand on essaie d’intégrer un champ fonctionnel d’un élément maître sur un élément esclave beaucoup plus étendu (non compatible). En inversant le choix des surfaces, le résultat est clairement meilleur ( à droite), mais encore sans garantie d’avoir une intégration satisfaisante. Par la suite, nous décrirons la méthode adoptée pour Code_Aster , basée sur une subdivision de l’élément qui ne sert qu’à l’intégration . Cette méthode, qui consiste à raffiner la surface d’intégration, limite l’erreur d’intégration numérique, et a surtout l’avantage d’être extrêmement simple à mettre en œuvre.
Plusieurs voies ont été explorées pour tenter de minimiser ces problèmes d’incompatibilités entre les surfaces. Celle qui a été retenue dans Code_Aster est de proposer différents schémas d’intégration numérique, utilisables suivant les besoins. Il y a plusieurs types de schémas:
Intégration numérique aux nœuds;
Intégration numérique aux points de Gauss ;
Intégration numérique par la méthode de Simpson à 3, 5 ou 9 points d’intégration;
Intégration par la méthode de Newton-Cotes à 4, 5 ou 10 points d’intégration;
L’intégration de Simpson permet d’améliorer les résultats en cas d’incompatibilité de maillage entre les deux surfaces de contact, grâce à la subdivision des éléments d’intégration que nous avons déjà évoquée plus haut.
L’intégration de Newton-Cotes permet d’intégrer exactement des polynômes d’ordre élevé (supérieur à trois). Elle s’avère très utile lorsque les fonctions d’interpolation utilisées sur l’une des deux surfaces de contact sont d’ordre supérieur à un (\(\mathit{P1}-\mathit{P2}\) , \(\mathit{P2}-\mathit{P2}\) …). De plus, les quadratures de Newton-Cotes permettent aussi (comme celles de Simpson) une subdivision de l’élément d’intégration, et s’appliquent donc en cas d’incompatibilité de maillages.
Toutefois, on se gardera bien de croire qu’une utilisation systématique d’une quadrature plus riche (Simpson, Newton-Cotes) est toujours un bon choix. D’abord parce que, dans la plupart des cas, la quadrature aux nœuds s’avère en pratique suffisante (\(\mathit{P1}-\mathit{P1}\) avec des maillages suffisamment fins). Ensuite parce que plus une quadrature est riche (en nombre de points d’intégration), plus les volumes de calcul sont importants, et l’on se rendra bien vite compte que, par exemple, l’utilisation de la quadrature de Newton-Cotes est vite rédhibitoire. Il n’y a donc pas de règle systématique a priori quant au choix d’une quadrature, hormis les quelques indications que nous avons données: il faut tenir compte de la spécificité du modèle et du problème mécanique étudiés, et faire un compromis entre la précision voulue pour la solution et le temps de calcul nécessaire.
Éléments de contact#
Après avoir discrétisé la formulation variationnelle mixte et après avoir donné les schémas d’intégration numérique, il ne nous reste plus qu’à mettre en place des éléments finis pour le contact-frottement.
L’idée d’utiliser des éléments finis n’est pas immédiate. Par exemple, dans les méthodes discrètes de contact frottement [R5.03.50], les conditions cinématiques de contact sont imposées sans utiliser d’éléments finis et en intervenant directement sur le système algébrique à résoudre. Une telle stratégie a l’inconvénient de rendre la programmation et l’architecture complètement dépendante de la formulation, ce qui interdit au contact de bénéficier d’améliorations (architecture, performance) qui seraient réalisées dans le reste du code éléments finis. A contrario, nous allons voir que développer des éléments finis de contact-frottement implique de sortir du schéma standard des éléments finis et que c’est une stratégie sous-optimale du point de vue des performances.
Pour déterminer la matrice de rigidité (et le second membre) du système discrétisé, on doit calculer les termes élémentaires des éléments finis avant d’en effectuer l’assemblage. Ce niveau élémentaire (classique) ne suffit pas pour le traitement de notre problème de contact-frottement. Plusieurs incompatibilités viennent rendre le problème non-standard au sens des éléments finis.
Incompatibilité géométrique#
On constate déjà une incompatibilité géométrique () entre les discrétisations des deux surfaces de contact en vis-à-vis, un élément esclave, pouvant par exemple chevaucher plusieurs éléments maîtres.
|
Figure 8 : Incompatibilité géométrique |
Incompatibilité des statuts#
Le caractère discontinu des états de contact/frottement entraîne une incompatibilité des statuts (). Au sein d’un même élément de contact, peuvent coexister plusieurs états de contact-frottement;
|
Figure 9 : Incompatibilité des statuts |
Cette incompatibilité se règle facilement car on fait l’hypothèse que le champ de signe (indiquant donc le statut de contact et de frottement) est discrétisé par collocation et que les points de collocation ont été choisis pour être ceux de l’intégration numérique. Ce qui veut dire que le statut de contact n’est pas porté par un élément de contact, mais par un point d’intégration d’un élément de contact.
Construction des éléments de contact#
Il reste à construire ces éléments de contact. Nous allons illustrer le principe sur un cas 2D.
|
Figure 10 : Exemple 2D d’appariement entre deux surfaces incompatibles |
Sur la figure , nous représentons un exemple 2D de surfaces géométriquement incompatibles. Les croix représentent les points d’intégration et voici les appariements et les statuts de contact:
Point intégration |
Élément esclave |
Élément maître apparié |
Statut |
\(a\) |
\({e}_{1}-{e}_{2}\) |
\({m}_{2}-{m}_{3}\) |
Adhérent |
\(b\) |
\({e}_{1}-{e}_{2}\) |
\({m}_{2}-{m}_{3}\) |
Glissant |
\(c\) |
\({e}_{2}-{e}_{3}\) |
\({m}_{3}-{m}_{4}\) |
Glissant |
\(d\) |
\({e}_{2}-{e}_{3}\) |
\({m}_{4}-{m}_{5}\) |
Pas de contact |
\(e\) |
\({e}_{3}-{e}_{4}\) |
\({m}_{4}-{m}_{5}\) |
Pas de contact |
\(f\) |
\({e}_{3}-{e}_{4}\) |
aucun |
Pas de contact |
Tableau 6.2.3-1
Il y a donc trois éléments de contact qu’on représente sur la figure . Chaque élément est géométriquement défini par les nœuds esclaves et maîtres. Le statut est représenté par le point de collocation qui agit comme un point d’intégration pour chaque élément de contact.
|
Figure 11 : Création des trois éléments de contact |
Nous avons mis les points de collocation \((a,b,c,d)\) comme des points de Gauss, mais leur emplacement va bien entendu dépendre du schéma d’intégration. Chaque élément de contact représente donc une matrice ou un vecteur élémentaire sur les degrés de liberté portés par les nœuds esclaves et par les nœuds maîtres. En 3D, la plus grande maille de surface est le QUAD9 (quadrangle à neuf nœuds), ce qui fait \(2\times 9\times 3=54\) degrés de liberté pour les déplacements. Dans le cas du contact avec frottement en 3D, on ajoute sur la surface esclave trois degrés de liberté supplémentaires par nœud: le Lagrange du multiplicateur de contact et les deux Lagrange du semi-multiplicateur de frottement, ce qui donne \(9\times 3=27\) degrés de liberté supplémentaires. Donc, en 3D, pour l’élément de contact avec frottement, la taille maximale de la matrice élémentaire est de \(81\) par \(81\) degrés de liberté.
En pratique la méthode est la suivante:
Création de l’élément de contact sur les nœuds maîtres et esclaves;
Pour chaque élément de contact, boucle sur les points de collocation;
Pour chaque point de collocation, calcul de la matrice (ou du vecteur) élémentaire.
C’est à cause du dernier point que la création de l’élément fini est sous-optimal. En effet, on appelle le calcul de la matrice élémentaire pour chaque point de collocation. Or, pour les éléments finis classiques, la boucle sur les points d’intégration est à l’intérieur du calcul élémentaire.
Nous voyons donc que l’usage d’un schéma de collocation permet de faire coexister des états de contact différents au sein du même élément. Toutefois, il subsiste un problème de compatibilité entre l’approximation EF des fonctions de forme et le schéma d’intégration. Prenons l’exemple d’un segment (figure ).
|
Figure 12 : Incompatibilité EF/collocation sur un segment |
Sur ce segment, le déplacement est approché par une fonction de forme linéaire et il vaut 0 au nœud de gauche et il est supérieur à zéro au nœud de droite. Du point de vue du contact, ça signifie qu’un point situé au milieu du segment va nécessairement se décoller. Or le schéma d’intégration par collocation n’indique pas forcément cet état. Il y a un découplage entre l’approximation polynomiale des déplacements et la discrétisation par collocation des champs de signe. Ce qui signifie que la méthode de collocation est parfaite pour des zones avec des états de contact constants sur l’élément mais qu’elle est beaucoup moins pertinente si l’état de contact varie. Une telle situation génère des oscillations dans les zones de transition. Ces oscillations ne peuvent pas être réduites en utilisant plus de points de collocation mais en raffinant le maillage dans cette zone.
Remarques:
Les éléments de contact-frottement n’ayant pas de sens géométrique, ils ne sont pas conformes, sans aucun impact sur la qualité du résultat: ils sont infiniment aplatis lorsqu’il y a contact, ils peuvent être très distordus et avoir des connectivités «étranges» dans le cas du raccordement entre des surfaces de degrés différents;
Le fait d’utiliser le schéma d’intégration aux nœuds ne rend pas les éléments de contact plus performants du point de vue de leur utilisation, puisque le principe décrit ci-dessus est générique pour tous les schémas d’intégration;
Il n’y a pas de procédure dédiée aux maillages parfaitement compatibles. Deux maillages choisis pour être compatibles ne donneront pas des résultats plus rapides que deux maillages incompatibles;
Le fait de devoir construire des éléments finis tardifs se traduit par une refactorisation de la matrice globale à chaque itération de Newton, puisque la topologie de la matrice va être modifiée.
Préparation de la discrétisation#
Discrétisation en temps#
La discrétisation (pseudo)-temporelle sera indiquée par l’indice \(k\) . Par exemple, le déplacement \({u}_{k}\) au pas de temps courant s’écrit comme la somme du déplacement cumulé \({u}_{k}^{\mathit{cum}}\) depuis le début du processus de Newton du pas de temps courant et de l’incrément de déplacement de l’itération courante de Newton \(\Delta u\) que l’on cherche:
Cette écriture itérative est aussi valable pour les autres inconnues nodales (champ des multiplicateurs de Lagrange):
Et:
Dans Code_Aster , on a fait le choix de travailler en quantité totale (et non incrémentale) pour les multiplicateurs de Lagrange, ils sont donc initialisés à zéro à chaque début de pas de temps et l’on a donc:
Et:
On utilisera une notation qui sépare au pas de temps courant \(k\) les quantités connues \({\stackrel{ˆ}{u}}_{k}\) et les quantités inconnues \(\Delta u\) :
Pour les coordonnées actualisées, comme on a \(x=p+u\) :
Pour le champ de pression de contact:
Pour le champ de pression de frottement:
Cette notation reste valable pour d’autres quantités comme le multiplicateur augmenté de contact:
Et le semi-multiplicateur augmenté de frottement:
Pour l’opérateur de projection sur la boule unité:
Discrétisation du repère local#
On ne précise pas la manière dont sont interpolées les normales et les tangentes, car tout dépendra des options précisées dans DEFI_CONTACT (activation du lissage par exemple). Il faut simplement savoir qu’on utilise les mêmes fonctions de forme que la description de la géométrie. On écrira directement que la normale ou les tangentes sont des quantités exprimées au point d’intégration courant:
On aura besoin de transformer des inconnues de l’espace 3D usuel vers le plan tangent. On utilise pour ça la matrice rectangulaire \(\left[{T}_{h}\right]\) telle que:
La discrétisation de la matrice de projection normale se note:
En développant les termes de la normale \(n\) (discrétisée) et ses trois composantes dans le repère global \((x,y,z)\) , on a:
Enfin, la discrétisation de la matrice de projection tangente se note:
Attention à ne pas confondre \(\left[{P}_{h}^{\tau}\right]\) avec la matrice \(\left[{T}_{h}\right]\) () sur le plan tangent. \(\left[{P}_{h}^{\tau}\right]\) vaut:
Les matrices \(\left[{P}_{h}^{\tau}\right]\) et \(\left[{P}_{h}^{n}\right]\) sont symétriques:
Discrétisation des champs solutions#
Considérons le champ vectoriel décrivant la géométrie (coordonnées des nœuds) sur l’élément esclave \({p}^{e}\) . Ce champ a deux ou trois composantes suivant la dimension de l’espace, nous le noterons sous la forme suivante:
La composante \({p}_{x}^{e}\) est approximée par:
Avec \({\mathit{NN}}_{e}\) le nombre de nœuds sur l’élément esclave, \({p}_{x,j}^{e}\) la J-ième valeur nodale du champ suivant la dimension \(x\) et \({N}_{j}^{e}\) les fonctions de forme. En projetant sur les bases éléments finis, les fonctions de forme étant identiques pour les trois dimensions de l’espace:
On écrira sous forme plus compacte les valeurs nodales \(\left\lbrace {p}_{h}^{e}\right\rbrace\) :
De manière analogue pour l’élément maître:
Pour les champs de déplacements:
Et pour la coordonnée courante:
Le champ de multiplicateur de contact est scalaire:
Avec \({\psi}_{j}\) les fonctions de forme utilisées pour les multiplicateurs de Lagrange. On écrira sous forme plus compacte:
La normale et les tangentes doivent également être discrétisées. On part de la définition du repère tangent sur la maille maître, et on interpole les normales et les tangentes grâce aux mêmes fonctions de forme. On notera \(\left\lbrace {n}_{h}\right\rbrace\) , \(\left\lbrace {t}_{h}^{1}\right\rbrace\) et \(\left\lbrace {t}_{h}^{2}\right\rbrace\) ces quantités discrétisées. Lors de l’écriture des différentes relations, nous avons supposé que la densité d’effort normal de frottement était une quantité globale, définie dans le repère global:
La composante \({\Lambda}_{x}\) sur un élément de bord, s’écrit comme suit:
Avec \({\mathit{NN}}_{c}\) le nombre de nœuds sur l’élément de bord, \({\Lambda}_{x,j}\) la J-ième valeur nodale de la composante \(x\) du champ de pression de frottement et \({\phi }_{j}\) les fonctions de forme. En fait, comme il existe naturellement une relation étroite entre le contact et le frottement (celle de la définition du repère local), les fonctions de base \({\psi}_{j}\) et \({\phi }_{j}\) ne sont pas indépendantes, et on décomposera le semi-multiplicateur de Lagrange de frottement global \(\Lambda\) en deux composantes (en 3D) sur le plan tangent défini au point de contact par \(({t}^{1},{t}^{2})\) :
Que l’on va ré-écrire sous forme plus compacte:
Si on considère le cas 3D, le champ du semi-multiplicateur de frottement \(\Lambda\) est donc de dimension trois fois le nombre de nœuds esclaves. Sa version discrète \(\left\lbrace {\Lambda}_{h}\right\rbrace\) , projection sur le plan tangent est de dimension deux fois le nombre de nœuds esclaves. Nous n’avons donc plus que trois bases éléments finis a priori indépendantes:
La base décrivant la géométrie de la structure;
La base décrivant le champ primal d’inconnues (déplacement ou vitesse);
La base décrivant le multiplicateur de Lagrange du contact et la projection de semi-multiplicateur de frottement sur le plan tangent au point de contact.
Intégration numérique#
La transformation de l’intégrale continue d’une quantité scalaire \(f\) en intégrale discrète s’écrit:
Avec \(\mathit{NP}\) le nombre de points d’intégrations, \({\zeta}_{c}\) la coordonnée paramétrique du point d’intégration sur la maille esclave, \({\omega}_{c}\) le poids d’intégration de ce point et \({J}_{c}\) le jacobien de la transformation du domaine d’intégration. Ce point est projeté de la surface esclave vers la surface maître par l’application d’appariement \({A}_{t}\) . L’évaluation sur la surface maître des quantités «appariées» notées \(\stackrel{ˉ}{f}\) , s’écrit donc simplement:
Dans le cas du contact, il convient de remarquer que chaque point de contact est aussi un point d’intégration (méthode de discrétisation du champ de signe par collocation). L’ensemble des intégrales continues seront donc discrétisées par un schéma numérique confondu avec la méthode de collocation des champs de signe.
Discrétisation des sauts#
Les différentes quantités liées au «saut» sont des fonctions d’une quantité sur la surface esclave avec sa projection sur leur surface maître. Considérons d’abord la projection du saut sur la normale:
En rappelant que \(u\) est une quantité définie sur la surface esclave et \(\stackrel{ˉ}{u}\) sa projection sur la surface maître. Nous aurons aussi besoin de développer les termes de type tangentiels:
On rappelle la discrétisation définie sur l’inconnue \(u\) par () utilise les fonctions de forme définies au point d’intégration courant \({\zeta}_{c}\) sur la surface esclave. La projection \(\stackrel{ˉ}{u}\) de cette quantité sur la surface maître utilisera donc le point d’intégration transporté \({A}_{t}({\zeta}_{c})\) sur la surface maître:
Que l’on notera de manière plus compacte par:
On a finalement:
C’est un scalaire. Pour le saut tangentiel:
Si on considère maintenant le saut normal de déplacement \({d}_{n}\) :
Donc:
De même:
D’où:
Les quantités dérivées par rapport aux coordonnées paramétriques:
En notant \(\left[{\stackrel{ˉ}{B}}_{\alpha}^{m}\right]\) la matrice des dérivées des fonctions de forme par rapport à la coordonnée paramétrique \({\zeta}^{\alpha}\) . Soit:
Discrétisation des (semi)-multiplicateurs#
Considérons le multiplicateur augmenté de contact \({g}_{n}\) :
Sa discrétisation donne:
Considérons maintenant le semi-multiplicateur augmenté de frottement \({h}_{\tau}\) :
Nous considérons l’écriture incrémentale en déplacement (voir § 5.2 ):
La discrétisation du semi-multiplicateur donne:
Statuts de contact et de frottement#
On a défini le champ de signe de contact noté \({S}_{u}^{{g}_{n}}\) par l’expression suivante:
Cet opérateur est la fonction caractéristique du multiplicateur augmenté de contact \({g}_{n}\) sur le convexe \({\mathrm{ℝ}}^{-}\) . On a également défini le champ de signe de frottement noté \({S}_{f}^{{h}_{\tau}}\) par l’expression suivante:
Cet opérateur est la fonction caractéristique du semi-multiplicateur augmenté de frottement \({h}_{\tau}\) sur le convexe \(B(0,1)\) . La discrétisation des champs des signes est triviale car on utilise la collocation sur ces champs de nature discontinues (voir § 6.1.4 ). La valeur discrétisée d’un multiplicateur (de contact ou de frottement) est égale à sa valeur au point d’intégration. Ces opérateurs étant discontinus, leur linéarisation n’a pas de sens. On écrit donc directement la valeur de ces champs au temps courant \(k\) mais juste avant l’itération de Newton courante:
Discrétisation des quantités virtuelles#
Nous allons maintenant considérer l’écriture des formes discrétisées des quantités virtuelles. Pour \(\delta {d}_{n}\) :
Pour \(\delta {d}_{\tau}\) :
La variation du semi-multiplicateur augmenté de frottement vaut:
Système discrétisé résultant#
On repart du système (), qui représente l’équilibre du système (PTV) sous forme incrémentale:
Après discrétisation, nous obtenons le système linéarisé suivant:
Ce système est linéarisé (voir § 5.4 ):
Après linéarisation, on se retrouve donc avec un système de la forme suivante:
Avec le vecteur des inconnues nodales (réelles ou virtuelles) :math:`langle delta Wrangle ` :
La matrice \([K]\) :
Le vecteur \(\left\lbrace L\right\rbrace\) :
Les matrices \(\left[{K}_{\mathit{ee}}^{t}\right]\) et \(\left[{K}_{\mathit{mm}}^{t}\right]\) sont les matrices tangentes provenant de la linéarisation des forces internes (discrétisation des matrices jacobiennes \({J}_{uu}^{u}\) , voir p. 11.1.1 ). Les autres matrices proviennent de la linéarisation des termes de contact et de frottement (voir § 5.4 ). Le système à résoudre sera donc:
Dans les expressions développées, on considérera toujours le cas du Newton généralisé, en se basant sur la discrétisation «incrémentalo-itérative» précisée dans le § 6.3.1 . Les quantités discrétisées seront écrites en un seul point d’intégration de coordonnée \({\zeta}_{c}\) , de telle manière qu’on introduira le poids d’intégration \({\omega}_{c}\) de ce point et le jacobien \({J}_{c}\) de la transformation du domaine d’intégration:
Remarque importante: Il convient de remarquer qu’avec les notations ci-dessus c’est en fait \(-\left[L\right]\) que l’on trouve à droite du signe égal. Cependant comme les efforts de contact-frottement sont ajoutés dans le code aux forces internes qui sont elles-mêmes retranchées au second membre, on n’a pas besoin de tenir compte de ce changement de signe. Les expressions qui sont écrites dans les termes élémentaires du code sont donc bien identiques à celles présentées.
L’expression du vecteur second membre se retrouve en annexe à la p. 69 et celle de la matrice tangente à la p. 72 .
Gestion des incompatibilités dans le traitement du contact#
La présence de contact et de conditions aux limites de Dirichlet imposées sur des points d’une surface commune peut générer des problèmes de redondance aboutissant à l’échec du calcul (en général, la matrice de rigidité tangente devient singulière). Dans ce qui suit, nous présentons quelques cas typiques de redondance que nous avons identifiés, et les traitements spécifiques proposés dans [].
Redondance entre les conditions aux limites et les conditions de contact-frottement#
Exemple#
Le cas le plus courant est une redondance des conditions de contact avec les conditions aux limites. Prenons l’exemple d’un bloc en contact frottant sur un plan rigide (cas-test ssnv128, voir ): on modélise la moitié de ce bloc, à l’extrémité de laquelle on applique les conditions de symétrie: \({u}_{x}=0\) . Le chargement est une pression imposée \(F\) sur les faces de coté et de dessus. La face du bas est, comme on l’a dit, en contact frottant avec son support. Sur le segment \([\mathit{AC}]\) , on va avoir redondance entre la condition aux limites et la condition de glissement.
Lorsque l’on a discrétisé la formulation variationnelle mixte, on a écrit qu’il y avait couplage entre les déplacements et les semi-multiplicateurs de frottement (matrices \(\left[{A}_{u}^{\Lambda}\right]\) et \(\left[{A}_{\Lambda}^{u}\right]\) ), ce qui veut dire que la condition de non-glissement va donner \({u}_{x}=0\) , alors que cette dernière condition est déjà imposée dans le système global (la matrice \([B]\) des conditions aux limites). On impose donc deux fois la relation \({u}_{x}=0\) . Cette situation aboutit à une erreur de pivot nul lors de l’inversion de la matrice de rigidité tangente.
|
Figure 13 : Plaque en contact frottant sur un plan rigide |
Choix de l’espace d’approximation#
La solution retenue est de toujours donner la préséance aux conditions aux limites sur les conditions d’adhérence-glissement. Pour cela, il suffit de découpler les multiplicateurs de frottement et les déplacements pour la composante du nœud concerné par la redondance (pour reprendre notre exemple, on découplera \({u}_{x}\) et \({\Lambda}_{x}\) aux nœuds \(A\) et \(C\) ). Cela revient à modifier l’espace d’approximation du champ du multiplicateur de frottement en utilisant une fonction qui vaut zéro pour le degré de liberté idoine. Cette opération est illustrée sur la figure dans le cas 2D où la condition de redondance est au point \(A\) .
On prend \({H}_{h}\) pour espace d’approximation du semi-multiplicateur de frottement \({\mu}_{h}\) :
Cet espace est conforme du point de vue de la condition LBB. En pratique, on met simplement des zéros aux emplacements idoines des matrices \(\left[{A}_{e}^{\Lambda}\right]\) et \(\left[{A}_{\Lambda}^{e}\right]\) [bib3] .
L’option SANS_GROUP_NO_FR/SANS_NOEUD_FR de Code_Aster permet de supprimer la contribution à la matrice de frottement d’un ensemble de nœuds, tout en gardant les conditions aux limites et de contact unilatéral en ces nœuds.
|
Figure 14 : Modification des espaces d’approximation pour la redondance frottement/condition limite |
En 2D, une et une seule direction de glissement existe. La suppression de la condition d’adhérence pour les points problématiques est donc suffisante pour éliminer toute redondance. Elle se fait en découplant le semi-multiplicateur de frottement associé à une liaison des inconnues de déplacement.
En 3D par contre, le glissement ou l’adhérence se produisant dans le plan tangent, on peut vouloir privilégier une direction de glissement particulière (par exemple perpendiculaire à un blocage). Pour cela, l’utilisateur fournit sous le mot-clé DIRE_EXCL_FROT la direction à exclure. Cette direction permet la construction d’un repère local qui facilite le découplage d’un seul semi-multiplicateur de frottement. Si aucune direction n’est fourni alors cela revient à ne pas traiter le frottement sur les nœuds concernés (découplage des deux multiplicateurs).
Cyclage pendant le processus de Newton#
Problématique#
L’utilisation d’un algorithme de Newton sur des quantités non différentiables fait parfois apparaître des phénomènes de cyclage numérique. De manière rigoureuse, on peut utiliser l’algorithme de Newton sur des quantités quasi-différentiables (au sens classique), en l’occurrence, on parlera ici de B-différentiabilité [], c’est-à-dire que les points de non-différentiabilité appartiennent à un espace de mesure nulle. En pratique, avec les approximations numériques introduites par la discrétisation et le passage à une algèbre en virgule flottante fait qu’il est «fort peu probable» que l’on rencontre ces points de non-différentiabilité, mais des problèmes de cyclage ne sont pas exclus, particulièrement dans le cas où le contact/frottement se couple avec d’autres non-linéarités.
Un cyclage correspond au passage d’un point de non-différentiabilité lors du processus de Newton (c’est-à-dire d’une itération de Newton à l’autre).
Dans [], on en distingue de trois sortes (voir figure ) :
Le cyclage sur le statut de contact: un point du maillage passe alternativement d’un cas de contact à un cas de non-contact;
Le cyclage sur le statut de frottement: un point du maillage passe alternativement d’un cas de glissement à un cas d’adhérence;
Le cyclage sur le seuil de glissement: un point du maillage passe alternativement d’un cas de glissement avant à un cas de glissement arrière;
D’un point de vue plus heuristique, on montre également que les cycles d’ordre 2 causent souvent la non-convergence.
Figure 15: Les différents cyclages de contact/frottement
Détection#
Les trois sortes de cycle sont détectés automatiquement.
Traitement#
Traitement automatique du cyclage#
Il y a la possibilité de traiter automatiquement les cycles de type contact/décollement, adhérence/glissement et glissement avant/glissement arrière (mot-clef ADAPTATION dans DEFI_CONTACT).
Pour le mot-clef ADAPTATION=”CYCLAGE”, on construit une base du jacobien généralisé sur deux itérations successives: on se sert de cette base pour modifier les matrices de contact et de frottement courantes. Cette manipulation sert pour le cyclage contact/décollement. Les cas de cyclage glissement avant/glissement arrière sont traités de sorte que le point est pris adhérent dès que la norme de son glissement est presque nul. Dans le cas où l’utilisateur demande un traitement type cyclage, on fait un basculement automatique des points non convergentsen pénalisation avec un contrôle automatique du coefficient de pénalisation (si l’algorithme initial était du Lagrangien augmenté). La formulation variationnelle en ALGO_CONT=”STANDARD’est ainsi modifiée. En effet, dans l’expression de la formulation variationnelle de départ. En notant, \({G}_{\text{c},k}^{\mathit{STD}}\) la force de contact faible liée ALGO_CONT=”STANDARD” pour l’ensemble des points qui ne posent pas de problèmes, \({G}_{\text{c},k}^{\mathit{PENA}}\) pour les points qui ont basculé, on a:
Enfin, en mode cyclage, pour renforcer la robustesse, on n’attend pas la convergence de Newton qui assure que les pressions et les jeux sont bien dans le cône des valeurs admissibles \({\mathrm{ℝ}}^{-}\) .
Pour le mot-clef ADAPTATION=”ADAPT_COEF”, il y a deux cas:
En lagrangien augmenté (ALGO_CONT=”STANDARD” et ALGO_FROT=”STANDARD”), on va modifier le coefficient d’augmentation pour «casser» le cycle entre le statut adhérent et le statut glissant au même point d’intégration. Cette modification du coefficient est sûre car elle ne change pas les résultats, mais juste la vitesse de convergence. Ce mode sert pour gérer le cyclage adhérent/glissant. Pour cela, l’algorithme suivant est employé:
Pour un coefficient d’augmentation donné, on calcule le statut adhérent/glissant au cycle précédent et au cycle courant du point détecté comme «cyclique»: est-on proche du point de discontinuité (tolérance fixée à 5%)?
Si oui, on modifie itérativement le coefficient d’augmentation pour sortir de cette zone critique tout en ne s’en rapprochant pas de l’autre coté;
Si non, on ne change rien;
On propage le coefficient adapté sur les autres points selon un algorithme spécifique (voir § 8.3.2 );
En lagrangien stabilisé (ALGO_CONT=”PENALISATION”), on va modifier les coefficients de régularisation en pénalisation pour détecter automatiquement les cas de matrices singulières et traiter automatiquement la pénétration maximum calculée au cours des itérations: voir § 8.3.3 .
Remarques sur le cyclage adhérent/glissant :
Du fait de la mise à l’échelle de la projection du lagrangien de frottement (voir § 3.2.4 ), la tolérance donnée (ici 5% de part et d’autre du point de discontinuité) ne dépend pas des unités du problème;
La procédure de détection du cyclage a exhibé un point qui varie entre le statut glissant et le statut adhérent. La stratégie retenue ici consiste à considérer les cas où l’un des statuts est proche du point de discontinuité. Si ce n’est pas le cas (cyclage entre deux zones éloignées de la discontinuité ou entre deux zones proches de la discontinuité), l’algorithme ne fait rien.
L’algorithme d’adaptation du coefficient d’augmentation est une simple boucle qui, dans un premier temps, augmente le coefficient d’augmentation en le multipliant par deux, et si on échoue, repart dans l’autre sens en divisant le coefficient par deux. Le coefficient est choisi dans l’intervalle \([{10}^{-8}{,10}^{+8}]\) .
Algorithme de propagation du coefficient d’augmentation cas adhérent/glissant#
L’algorithme précédent n’est efficace que si on met en place une procédure qui va propager (ou pas) le coefficient optimal trouvé sur un point cyclique. En effet, en pratique, un tel coefficient optimal a de fortes chances de l’être aussi pour d’autres points. Cette méthode permet ainsi d’anticiper les risques de cyclage des points qui ne sont pas encore en contact et rend l’algorithme beaucoup plus performant.
Pour cela, lors de la procédure d’adaptation du coefficient (algorithme du § 8.3.1 ), on stocke les valeurs maximales et minimales trouvées sur toute la zone de contact et on applique l’algorithme suivant:
Boucle sur les points de contact
Avec le coefficient maximal sur la zone de contact, on évalue sur tous les points si on est proche ou pas du point de discontinuité, on a alors quatre situations:
Ce nouveau coefficient ne modifie pas la situation du point sur le graphe de frottement par rapport à l’ancien coefficient→ on propage le nouveau coefficient;
Ce coefficient modifie la situation du point sur le graphe de frottement mais l’ancienne situation était loin de la discontinuité→ on ne propage pas le nouveau coefficient;
Ce coefficient modifie la situation du point sur le graphe de frottement, l’ancienne situation était proche de la discontinuité et la nouvelle position s’éloigne de la discontinuité→ on propage le nouveau coefficient;
Ce coefficient modifie la situation du point sur le graphe de frottement, l’ancienne situation était proche de la discontinuité et la nouvelle position change de statut (franchissement de la discontinuité)→ on ne propage pas le nouveau coefficient;
Avec le coefficient minimal sur la zone de contact, on évalue sur tous les points si on est proche ou pas du point de discontinuité, on évalue les mêmes situations que dans le cas précédent;
Fin de la boucle sur les points de contact
Algorithme d’initialisation du coefficient de pénalisation#
Nous allons détailler si l’utilisation du mot-clef ADAPTATION=”ADAPT_COEF” dans le cas du contact pénalisé. Dans certains calculs avec la méthode de pénalisation sur le contact, la matrice est singulière dès la première itération du premier pas de temps. Cette singularité est due à un mauvais choix des coefficients de pénalisation de départ. Pour détecter numériquement cette situation et la traiter, on s’appuie sur la matrice assemblée avant la première factorisation. En effet, on constate que cette situation apparaît si la plus grande valeur sur la diagonale (hormis les termes de Lagrange dus aux conditions aux limites) est très éloignée de la plus petite valeur sur cette même diagonale.
Le profil de la matrice assemblée en regardant la diagonale est:
Avec \(\stackrel{~}{K}\) la matrice de rigidité globale réduite aux termes diagonaux sans les termes de contact et \({K}_{\mathit{pena}}\) la matrice de pénalisation du contact telle que:
On omet volontairement les termes de frottement par souci de simplicité. Dans cette dernière expression, on rappelle que:
\({S}_{c}\) désigne le statut de contact;
\({g}_{n}\) désigne le gap courant de contact;
\(\delta {g}_{n}\) désigne le gap virtuel de contact;
L’algorithme d’initialisation traite deux cas de figures actuellement :
Après avoir assemblée la matrice \(K\) , on calcule \({r}_{min}=min(\stackrel{~}{K}\cup {K}_{\mathit{pena}})\) et \({r}_{max}=max(\stackrel{~}{K}\cup {K}_{\mathit{pena}})\)
Premier cas de figure: \({r}_{min}\) est très petit mais non nul et \(\left|(\log({r}_{max}))\right|/\left|(\log({r}_{min}))\right|<4\) .
Ce cas de figure correspond au cas le plus souvent rencontré lorsqu’il y a une matrice singulière détectée. Elle traduit deux hypothèses observées en pratique.
La première hypothèse est qu’on est en présence des matériaux non mous (rigidité globale de structure non négligeable donc \({r}_{min}\) très petit vient du contact) et qu’il y au moins un point en contact initialement (\({K}_{{\mathit{pena}}_{11}}\ne 0\) ). Si \({r}_{n}\) a été choisi trop grand (ex: \({r}_{n}\gg E\) module de Young du corps dur) alors soit \({K}_{{\mathit{pena}}_{22}}\) est nulle soit \({K}_{{\mathit{pena}}_{11}}\) est très supérieure à \(max(\stackrel{~}{K})\) : il faut diminuer \({r}_{n}\) . Une valeur inférieure à \({r}_{max}\) suffit pour bien initialiser le coefficient de pénalisation. Sinon si \({r}_{n}\) a été choisi trop petit alors soit \({K}_{{\mathit{pena}}_{11}}\) est nulle soit \({K}_{{\mathit{pena}}_{22}}\) est très supérieure à \(max(\stackrel{~}{K})\) : il faut augmenter \({r}_{n}\) . Une valeur inférieure à \({r}_{max}\) suffit pour bien initialiser le coefficient de pénalisation.
La deuxièmehypothèse est qu’on est en présencedes matériaux non mous (rigidité globale de structure non négligeable donc \({r}_{min}\) très petit vient du contact) et qu’il n’y a aucun nœud en contact initialement. Dans ce cas seul \({K}_{{\mathit{pena}}_{22}}\) est à l’origine de la singularité. Une valeur inférieure à \({r}_{max}\) suffit pour bien initialiser le coefficient de pénalisation.
En pratique une seule formule suffit pour ce cas de figure :
:math:`{r}_{{n}_{mathit{corrige}}}={10}^{(left |
(log({r}_{max}))right |
/3.0)}` |
Deuxième cas de figure: \({r}_{min}\) est nul ou \({r}_{max}\) très grand et \({r}_{min}\) non nul
Ce cas de figure est symptomatique du fait qu’il y a au moins un matériau souple dans le système mécanique. En pratique, il n’y a pas singularité de matrice observée mais pour améliorer la convergence on propose les formules suivantes:
Si \({r}_{min}\) nul:
\({a}_{min}\) est la plus petite arête dans la zone de contact.
Si \({r}_{max}\) très grand et \({r}_{min}\) n’est pas nul
:math:`{r}_{{n}_{mathit{corrige}}}={10}^{(left |
(log({r}_{min}))right |
)}` |
Conseils si le traitement automatique échoue ou n’est pas possible#
Pour l’utilisateur, l’apparition de ces cyclages est le signe d’une certaine difficulté à converger et on peut traiter ce problème de trois manières:
Utiliser un algorithme plus robuste comme le Newton partiel ou le Newton incomplet (voir mots-clefs ALGO_RESO_* ). Dans ce dernier cas, un autre type de cyclage peut apparaître: c’est le flip-flop, un cyclage de type contact mais d’ordre 15;
Modifier les paramètres d’augmentation (mots-clefs COEF_FROT et COEF_CONT). Ici, il n’y a pas de règles pré-établies sur le sens de la modification (diminuer ou augmenter la valeur ) et sur son amplitude. Néanmoins, il convient d’observer que le paramètre d’augmentation \({\rho}_{t}\) est très sensible.
Modifier son modèle. Soit en raffinant/déraffinant le maillage (attention! le raffinement peut amplifier les phénomènes de cyclage), soit en dégradant le modèle (passage à un matériau élastique, voire rigide).
Fonctionnalités et vérification#
Cette documentation décrit la modélisation du contact-frottement par la méthode dite «continue», pour son utilisation, on peut se référer à la documentation d’usage de la commande DEFI_CONTACT et au manuel de modélisation U2.04.04.
Il existe une soixantaine de cas-tests validant les différentes fonctionnalités liés au contact/ frottant , parmi lesquels:
Le cas-test du contact de Hertz: voir ssnv104;
Le cas-test du patin frottant: voir ssnv128;
Les cas-tests du benchmark NAFEMS sur le contact/frottement: voir ssna122, ssnp154, ssnp155, ssnp156 et ssnp157;
Le cas test «shallow ironing», voir ssnp153;
Bibliographie#
Alart P., Curnier A. A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput. Methods. Appl. Mech. Engrg., vol.9, pp.353-375, 1991.
Hallquist J.O., Goudreau G.L., Benson D.J. Sliding interfaces with contact-impact in large-scale langrangian computations. Comput. Methods. Appl. Mech. Engrg., vol.51, pp.107-137, 1985.
Laursen T.A., SimoJ.C. A continuum-based finite element formulation for the implicit solution of multibodies, large deformation frictional contact problems. Int. J. Num. Meth. Engng., vol.36, pp.3451-3485, 1993.
Ben-Dhia H., Zarroug M. Hybrid frictional contact particles-in elements . Revue Européenne des Éléments Finis, n°9, pp.417-430, 2002.
Zarroug M. Élémentsmixtes de contact frottant en grandes transformations et applications . Thèse de Doctorat de l’École Centrale Paris, 2002.
Zammali C. Contribution à la modélisation mécanique et numérique des problèmes de contact-impact. Thèse de Doctorat de l’École Centrale Paris, 2005.
Torkhani M. Contribution au développement numérique d’éléments de contact et modélisation de l’usure des structures minces. Thèse de Doctorat de l’École Centrale Paris, 2008.
Yastrebov V.A. Computational contact mechanics - Geometry, detection and numerical techniques . Thèse de Doctorat de l’École Supérieure des Mines de Paris, 2011.
Kudawoo D. Problèmes industriels de grande dimension en mécanique numérique du contact: performance, fiabilité et robustesse. Thèse de Doctorat de l’Université d’Aix-Marseille, 2012.
JeanM. The non-smooth contact dynamics method . Comput. Methods. Appl. Mech. Engrg., vol.177, pp.235-257, 1999.
Alart P. Méthode de Newton généralisée en mécanique du contact . J. Math. Pures Appl., vol.76, pp.83-108, 1997.
Klarbring. Large displacement frictional contact: a continuum framework for finite element discretization . Eur. J. Mech. A/Solids, vol.14, pp.237-253, 1995.
Annexe A : calcul des termes de la jacobienne#
Équation d’équilibre#
Pour l’équation d’équilibre les quantités virtuelles sont des variations totales par rapport aux déplacements. La notation \(\tilde{\delta}\) sera donc employée. La première ligne de la matrice jacobienne \(\left[J\right]\) est constituée de trois termes:
On sait déjà que les forces internes et externes des corps ne dépendront pas du lagrangien de contact, ni du lagrangien de frottement:
Mais également que la pression de contact ne dépend pas de la pression de frottement:
Alors que la pression de frottement dépend de la pression de contact (voir § 5.4.1 ):
Il ne s’agit pas d’hypothèses de simplification, mais de la conséquence logique de l’écriture du modèle physique.
Seconde variation par rapport aux déplacements#
Dans ce paragraphe, la seconde variation s’entendra comme variation totale par rapport aux déplacements \({u}^{i}\) . Elle sera notée \({\Delta}_{{u}^{i}}=\tilde{\Delta}\) . Nous n’allons pas revenir en détail sur le calcul des termes non linéaires \({J}_{uu}^{int}\) et \({J}_{uu}^{\text{ext}}\) correspondant aux forces internes \({G}_{int}^{i}\) et externes \({G}_{\text{ext}}^{i}\) puisqu’ils concernent le comportement et la cinématique des solides sans contact-frottement. Par exemple, la variation du travail des efforts internes introduit les matrices tangentes \({\underline{\underline{K}}}_{t}^{i}\) telles que:
Nous commençons par évaluer \({J}_{uu}^{\text{c}}\) écrite avec la variation de jeu normal:
Trivialement:
En développant l’expression:
Reprenons les expressions établies dans le § 2.3 . La première variation totale du jeu normal vaut (voir ()):
Et donc, par application de ():
La seconde variationest beaucoup plus complexe (voir ()):
On utilise le tenseur métrique \(\underline{\underline{m}}\) et le tenseur de la seconde forme fondamentale \(\underline{\underline{\kappa}}\) . Toutes ces expressions n’interviennent qu’en cas de contact (\({S}_{u}^{{g}_{n}}=1\) ) et donc quand \({d}_{n}=0\) . Cette hypothèse simplifie beaucoup l’expression. Pour (), en supprimant le terme en \({d}_{n}\) et avec ():
Pour trouver \(\delta {\zeta}_{\alpha}\) et \(\Delta {\zeta}_{\alpha}\) , il faut résoudre le système suivant (version contravariante de ()):
On a pris \({d}_{n}=0\) , ce qui simplifie les expressions de (). Finalement, on a:
Avec:
Et:
Nous allons maintenant considérer le terme d’équilibre \({J}_{uu}^{\text{f}}\) correspond à la réaction pour le frottement:
L’équation () utilise la variation totale du jeu tangent \({\tilde{\delta}d}_{\tau}\) par rapport aux déplacements. Nous allons faire des hypothèses de simplification. On va négliger la deuxième variation du jeu tangent, c’est-à-dire:
Ce qui nous donne:
On écrit la variation du semi-multiplicateur de Lagrange de frottement \({h}_{\tau}\) . En utilisant l’hypothèse cinématique ():
Et donc:
Ainsi:
De manière analogue:
Pour la variation du semi-multiplicateur normé \(\tilde{\Delta}\tau\) , qui correspond à la projection sur la boule unité, on ne considère que la variation en glissement (\({S}_{f}^{{h}_{\tau}}=0\) ) puisque le terme s’annule dans le cas adhérent. On a (voir ()):
Par application de la variation d’une fonction composée:
Ce qui nous donne:
Avec:
Et:
Seconde variation par rapport à la pression de contact#
Maintenant, la seconde variation s’entend comme variation par rapport à \({\lambda}_{n}\) . Nous commençons par évaluer \({J}_{uc}^{\text{c}}\) tel que:
On rappelle que la variation du jeu normal vaut ():
Et:
Il reste:
Nous considérons maintenantévaluer \({J}_{uc}^{\text{f}}\) tel que:
Nous avons supposé (variante 5, voir § 5.4.1 ):
Et nous avons:
Il reste:
Avec:
Et:
Seconde variation par rapport à la pression de frottement#
Finalement, la seconde variation s’entendra comme variation par rapport à \(\Lambda\) . Le seul terme non-nul est \({J}_{uf}^{\text{f}}\) tel que:
Les variations de \({h}_{\tau}\) par rapport à \(\Lambda\) valent :
On a:
De même, trivialement:
Ce qui nous donne:
Avec:
Et:
Loi de contact#
Pour la loi de contact (sous forme faible) les quantités virtuelles sont des variations par rapport à la pression de contact \({\lambda}_{n}\) . La notation \(\delta\) sera donc employée dans ce sens. La seconde ligne de la matrice jacobienne \(\left[J\right]\) est constituée de trois termes:
Comme la pression de contact ne dépend pas de la pression de frottement, on sait déjà que:
Seconde variation par rapport aux déplacements#
Dans ce paragraphe, la seconde variation s’entendra toujours comme variation totale par rapport aux déplacements \({u}^{i}\) . Elle sera notée \({\Delta}_{{u}^{i}}=\tilde{\Delta}\) . Nous commençons par évaluer \({J}_{cu}\) telle que:
On a \(\tilde{\Delta}{\delta \lambda }_{n}=0\) .Il ne reste donc que le terme en \(\tilde{\Delta}{d}_{n}\) dont l’expression est donnée par ():
On observe immédiatement que \({J}_{cu}\mathrm{\ne }{J}_{uc}\) à cause de (). En frottement, la matrice n’est pas symétrique.
Seconde variation par rapport à la pression de contact#
La seconde variation s’entendra comme variation par rapport à \({\lambda}_{n}\) . Nous calculons \({J}_{cc}\) telle que:
La seconde variation \({\Delta}_{{\lambda}_{n}}{\delta \lambda }_{n}=0\) est nulle et \({\Delta}_{{\lambda}_{n}}({\lambda}_{n})=\Delta {\lambda}_{n}\) , donc:
On remarqueraque ce terme n’intervientque lorsqu’il n’y a pas de contact (\({S}_{u}^{{g}_{n}}=0\) ).
Seconde variation par rapport à la pression de frottement#
La seconde variation s’entendra comme variation par rapport à \(\Lambda\) . Le terme \({J}_{cf}\) est nul:
Loi de frottement#
Pour la loi de frottement (sous forme faible) les quantités virtuelles sont des variations par rapport à la pression de frottement \(\Lambda\) . La notation \(\delta\) sera donc employée dans ce sens. La troisième ligne de la matrice jacobienne \(\left[J\right]\) est constituée de trois termes:
Seconde variation par rapport aux déplacements#
Dans ce paragraphe, la seconde variation s’entendra toujours comme variation totale par rapport aux déplacements \({u}^{i}\) . Elle sera notée \({\Delta}_{{u}^{i}}=\tilde{\Delta}\) . Nous commençons par évaluer \({J}_{fu}\) telle que:
On réutilise ():
On a aussi \(\tilde{\Delta}\delta \Lambda =0\) et \(\tilde{\Delta}\Lambda =0\) . Ainsi:
Avec:
Et:
Seconde variation par rapport à la pression de contact#
La seconde variation s’entendra comme variation par rapport à \({\lambda}_{n}\) . Nous commençons par évaluer \({J}_{fc}\) telle que:
En utilisant la variante 5 (voir § 5.4.1 ), on a \({\Delta}_{{\lambda}_{n}}({\tilde{G}}_{f})=0\) , ce qui implique \({J}_{fc}=0\) .
Seconde variation par rapport à la pression de frottement#
La seconde variation s’entendra comme variation par rapport à \(\Lambda\) . Nous commençons par évaluer \({J}_{ff}\) telle que:
Les variations de \({h}_{\tau}\) par rapport à \(\Lambda\) valent ():
La pression de contact ne dépend pas de la pression de frottement, donc:
On a aussi():
Finalement:
Il n’y a une contribution qu’en contact glissant et quand le contact est inactif.
Annexe B : vecteur second membre#
Les expressions considérées ici ne concernent que les termes liés au contact frottement. Il s’agit donc de discrétiser l’expression des travaux virtuels donnés par les équations () à (). On rappelle que le système à résoudre est:
Termes d’équilibre#
On commence par l’expression de la réaction de contact dans l’équation d’équilibre, à partir de \({G}_{c}\) () que l’on écrit sur la configuration connue au pas de temps \(k\) , depuis le début du processus de Newton:
Après discrétisation:
Le signe \(–\) provient de l’expression de l’équilibre dans (). On obtient, en utilisant () et en séparant les contributions esclave et maître:
Ces termes ne sont actifs que lorsqu’il y a contact (\({S}_{u,k}^{{g}_{n}}=1\) ).
On considère ensuite la réaction de frottement ():
En utilisant l’opérateur de projection (). Après discrétisation:
On obtient:
On observe un terme actif en contact glissant (\({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Et un terme actif en contact adhérent (\({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) :
Termes de la loi de Signorini#
Nous allons maintenant considérer le terme correspondant à l’inconnue en pression de contact, soit à partir de ():
Après discrétisation:
On obtient:
Cette expression peut se décomposer en deux parties. Si le contact est inactif (\({S}_{u,k}^{{g}_{n}}=0\) ), on a:
Si le contact est actif (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) ), on a des termes se simplifiant entre \({\stackrel{ˆ}{\lambda}}_{k}\) et \({\stackrel{ˆ}{g}}_{n,k}\) , ce qui donne:
Termes de la loi de Coulomb#
Enfin, le terme correspond à l’inconnue en pression de frottement, à partir de l’expression faible de la loi de Coulomb ():
Après discrétisation:
On obtient:
Cette expression complexe peut se décomposer en trois parties. Si le contact est inactif (\({S}_{u,k}^{{g}_{n}}=0\) ), on a:
Si on est en contact glissant (\({S}_{u,k}^{{g}_{n}}=1\) et \({S}_{f,k}^{{h}_{\tau}}=0\) ), on a:
Si on est en contact adhérent (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) ), on a, certains termes se simplifient entre la première et la deuxième ligne de (), en particulier, on a:
Avec:
Donc:
En résumé:
Sans contact |
\({\left\lbrace {L}_{c}^{c}\right\rbrace }_{\text{sans contact}}=-\frac{1}{{\rho}_{n}}{\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}\left\lbrace \psi \right\rbrace\) \({\left\lbrace {L}_{f}^{f}\right\rbrace }_{\text{sans contact}}={\omega}_{c}{J}_{c}{\left[\psi \right]}^{T}{\left[T\right]}^{T}\left\lbrace {\stackrel{ˆ}{\Lambda}}_{k}\right\rbrace\) |
Contact glissant |
\(\left\lbrace {L}_{c}^{e}\right\rbrace =-{\left[{N}^{e}\right]}^{T}\lbrace {n}_{h}\rbrace {\stackrel{ˆ}{g}}_{n,k}{\omega}_{c}{J}_{c}\) \({\left\lbrace {L}_{c}^{m}\right\rbrace }_{\text{glissant}}=+{\left[{N}^{m}\right]}^{T}\lbrace {n}_{h}\rbrace {\stackrel{ˆ}{g}}_{n,k}{\omega}_{c}{J}_{c}\) \({\left\lbrace {L}_{f}^{e}\right\rbrace }_{\text{glissant}}=-\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[{N}^{e}\right]}^{T}\left[{P}^{\tau}\right]\left\lbrace {\stackrel{ˆ}{\tau}}_{k}\right\rbrace\) \({\left\lbrace {L}_{f}^{m}\right\rbrace }_{\text{glissant}}=+\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[{\stackrel{ˉ}{N}}^{m}\right]}^{T}\left[{P}^{\tau}\right]\left\lbrace {\stackrel{ˆ}{\tau}}_{k}\right\rbrace\) \({\left\lbrace {L}_{c}^{c}\right\rbrace }_{\text{glissant}}=-{\stackrel{ˆ}{d}}_{n,k}{\omega}_{c}{J}_{c}\left\lbrace \psi \right\rbrace\) \({\left\lbrace {L}_{f}^{f}\right\rbrace }_{\text{glissant}}=\frac{1}{{\stackrel{ˉ}{\rho}}_{t}}\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[\psi \right]}^{T}{\left[T\right]}^{T}(\left\lbrace {\stackrel{ˆ}{\Lambda}}_{k}\right\rbrace -\left\lbrace {\stackrel{ˆ}{\tau}}_{k}\right\rbrace )\) |
Contact adhérent |
\({\left\lbrace {L}_{c}^{e}\right\rbrace }_{\text{adhérent}}=-{\left[{N}^{e}\right]}^{T}\lbrace {n}_{h}\rbrace {\stackrel{ˆ}{g}}_{n,k}{\omega}_{c}{J}_{c}\) \({\left\lbrace {L}_{c}^{m}\right\rbrace }_{\text{adhérent}}=+{\left[{N}^{m}\right]}^{T}\lbrace {n}_{h}\rbrace {\stackrel{ˆ}{g}}_{n,k}{\omega}_{c}{J}_{c}\) \({\left\lbrace {L}_{f}^{e}\right\rbrace }_{\text{adhérent}}=-\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[{N}^{e}\right]}^{T}\left[{P}^{\tau}\right]\left\lbrace {\stackrel{ˆ}{h}}_{\tau ,k}\right\rbrace\) \({\left\lbrace {L}_{f}^{m}\right\rbrace }_{\text{adhérent}}=+\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[{\stackrel{ˉ}{N}}^{m}\right]}^{T}\left[{P}^{\tau}\right]\left\lbrace {\stackrel{ˆ}{h}}_{\tau ,k}\right\rbrace\) \({\left\lbrace {L}_{c}^{c}\right\rbrace }_{\text{adhérent}}=-{\stackrel{ˆ}{d}}_{n,k}{\omega}_{c}{J}_{c}\left\lbrace \psi \right\rbrace\) \({\left\lbrace {L}_{f}^{f}\right\rbrace }_{\text{adhérent}}=-\mu {\stackrel{ˆ}{\lambda}}_{k}{\omega}_{c}{J}_{c}{\left[\psi \right]}^{T}{\left[T\right]}^{T}\left[{P}^{\tau}\right]\left\lbrace {\mathrm{〚}\stackrel{ˆ}{x}\mathrm{〛}}_{k}\right\rbrace\) |
Annexe C : matrice tangente#
Les expressions considérées ne concernent toujours que les termes liés au contact-frottement. Il s’agit donc de discrétiser les expressions linéarisées. Le processus de linéarisation fait apparaître des termes matriciels que nous allons discrétiser ici.
Quantités matricielles préliminaires#
On commence par décrire la forme discrétisée d’un certain nombre de quantités. La discrétisation de la direction \({\stackrel{ˆ}{\tau}}_{k}\) nous donnera:
Discrétisation de l’opérateur de projection sur la boule unité:
Le processus de linéarisation de la quantité relative à la boule unité (§ 3.2.4 , voir ()) nous donne:
On notera au passage l’identité:
La discrétisation utilise la définition (), avec () on peut donc écrire:
Quantités pour l’équilibre#
On commence par considérer les termes matriciels issus de la linéarisation de l’équation d’équilibre. La discrétisation des jacobiennes () et ():
On commence par ():
Donc:
On utilise :math:`left[{P}^{n}right]=lbrace {n}_{h}rbrace langle {n}_{h}rangle ` . Finalement:
À partir de ():
Considérons la première intégrale, en discrétisant, on obtient:
Avec les expressions:
Avec la matrice \(\left[G\right]\) telle que:
Considérons la seconde intégrale, en discrétisant, on obtient:
Avec les expressions:
Avec la matrice \(\left[H\right]\) telle que:
Enfin, la dernière contribution:
Avec les expressions:
Avec la matrice \(\left[L\right]\) telle que:
On considère maintenant la discrétisation des jacobiennes () et ():
La première partie concerne le cas adhérent (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) ):
La discrétisation de \({\underline{\underline{P}}}^{\tau}\) et des déplacements nous donne:
La seconde partie concerne le cas glissant (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Ce qui s’écrit, en utilisant en particulier la discrétisation de la projection sur la boule unité ():
Nous allons maintenant discrétiser les quantités concernant la seconde colonne, relatives à la pression de contact. À partir de la réaction de contact et de la jacobienne ():
On doit discrétiser:
Ce qui nous donne:
À partir de la réaction de frottement, c’est-à-dire des jacobiennes () et ():
La première partie concerne le cas adhérent (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) ):
La discrétisation de \({\underline{\underline{P}}}^{\tau}\) , des déplacements \(u\) et de la pression de contact \(\lambda\) nous donne:
La seconde partie concerne le cas glissant (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Ce qui s’écrit, en utilisant en particulier la discrétisation de la projection sur la boule unité ():
Nous allons maintenant discrétiser les quantités concernant la troisième colonne, relatives à la pression de frottement. À partir de la réaction de frottement et des jacobiennes () et ():
La première partie concerne le cas adhérent (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) ):
La discrétisation de \({\underline{\underline{P}}}^{\tau}\) , des déplacements \(u\) et de la pression de frottement \(\Lambda\) nous donne:
La seconde partie concerne le cas glissant (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Ce qui s’écrit, en utilisant en particulier la discrétisation de la projection sur la boule unité ():
Quantités pour la loi de contact#
On considère les termes matriciels issus de la linéarisation de la loi de contact et donc à la seconde ligne de la matrice globale. La discrétisation de la jacobienne () correspondant à la première colonne de la matrice globale:
À partir de la jacobienne:
On obtient deux quantités:
Par la discrétisation de la jacobienne (), on obtient l’expression des quantités correspondant à la seconde colonne:
À partir de la jacobienne:
On obtient:
Quantités pour la loi de frottement#
On considère les termes matriciels issus de la linéarisation de la loi de frottement et donc à la troisième ligne de la matrice globale. La discrétisation des jacobiennes () et () correspondant à la première colonne de la matrice globale:
La première partie concerne le cas adhérent (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=1\) ):
La discrétisation de \({\underline{\underline{P}}}^{\tau}\) , des déplacements \(u\) et de la pression de frottement \(\Lambda\) nous donne:
La seconde partie concerne le cas glissant (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Ce qui s’écrit, en utilisant en particulier la discrétisation de la projection sur la boule unité ():
La discrétisation de la jacobienne () correspondant à la troisième colonne de la matrice globale:
La jacobienne s’écrivant:
On préfère la diviser en deux. Une contribution pour le cas de contact glissant (\({\stackrel{ˆ}{S}}_{u,k}^{{g}_{n}}=1\) et \({\stackrel{ˆ}{S}}_{f,k}^{{h}_{\tau}}=0\) ):
Et une contribution pour le cas sans contact:
Écriture du glissement infinitésimal : GRAND_GLIS#
Dans le cas où le contact est établi et que le glissement est actif et que ce dernier est contraint par la loi de Coulomb, il faut linéariser la forme faible de la force induisant le glissement, la loi de Coulomb restant inchangée.
\({G}_{f}={\int}_{{\Gamma}_{c}}\mu {\widehat{\lambda}}_{k}{\widehat{S}}_{u,k}^{{g}_{n}}\left(1-{\widehat{S}}_{f,k}^{{h}_{\tau}}\right){\widehat{\tau}}_{k}\cdot {\delta d}_{\tau}d{\Gamma}_{c}\)
Jusqu’à présent, le glissement relatif a été approximé de la façon suivante:
Note: dans ce qui suit, on peut si besoin, remplacer \(\Delta ` par :math:\)delta ` .
En toute rigueur, si on adopte une écriture cinématique propre à la mécanique des milieux continus,la vitesse relative en un point de contact est :
Par la suite, on a:
\(\frac{\partial [{y}_{\mathit{escl}}-{y}_{\mathit{mait}}]}{\partial t}\Delta t-\frac{\partial {y}_{\mathit{mait}}}{\partial {\zeta}^{\alpha}}\Delta {\zeta}^{\alpha}=\frac{\partial {g}_{n}n}{\partial t}\Delta t-\frac{\partial {y}_{\mathit{mait}}}{\partial {\zeta}^{\alpha}}\Delta {\zeta}^{\alpha}=\frac{\partial {g}_{n}}{\partial t}\delta tn+{g}_{n}\frac{\partial n}{\partial t}\Delta t-\frac{\partial {y}_{\mathit{mait}}}{\partial {\zeta}^{\alpha}}\Delta {\zeta}^{\alpha}\)
On va noter:
\({v}_{n}=\frac{\partial {g}_{n}}{\partial t}\Delta tn\) et \({v}_{\tau}={g}_{n}\frac{\partial n}{\partial t}\Delta t-\frac{\partial {y}_{\mathit{mait}}}{\partial {\zeta}^{\alpha}}\Delta {\zeta}^{\alpha}\)
L’écriture continue de la vitesse relative est donc:
On s’intéresse ici à l’écriture de \({v}_{\tau}\) . En première approximation, et en supposant que le point est constamment en contact, on a:
On peut se référer à l’équation (66) pour prendre en compte la variation dûe à la normale \({g}_{n}\frac{\partial n}{\partial t}\) . D’un point de vue vocabulaire en mécanique des milieux continus, on dit que le point de contact esclave est décrite suivant une cinématique Lagrangienne tandis que son point projeté maître est décrit suivant une cinématique Eulerienne. L’équation précédente, montre qu’en négligeant le terme \({g}_{n}\frac{\partial n}{\partial t}\) , toute l’information nécessaire au calcul du glissement dans le plan tangent est contenue dans la variation du point projeté maître à condition que le statut du point contact soit égal à 1 et le reste durant le processus de Newton dans le pas de temps. On va expliciter le terme \(\Delta {\zeta}^{\alpha}\) .
Pour la linéarisation du frottement, on aura également besoin de la deuxième variation de \(\Delta \delta {\zeta}^{\alpha}\) .
On propose à l’utilisateur les deux écritures du glissement:
soit par projection de l’incrément de déplacement sur le plan tangent équation (4339) via le mot-clef GRAND_GLIS=”NON”,
soit en adoptant une écriture continue du glissement tangent relatif au point maître: GRAND_GLIS=”OUI”.
En théorie, en petits glissements, les deux écritures donnent le même résultat. Pour l’implémentation, la technique de discrétisation reste la même pour les inconnues cinématiques telles que décrite dans la section 6.