National audienceL'implantation numérique efficace des méthodes de traitement du contact reste encore un sujet de recherche à part entière. On détaille ici la mise en oeuvre d'une méthode de type Lagrangien augmenté. Afin de gagner en performance, l'algorithme comporte une procédure particulière de reconnaissance (séquentiel et parallèle) de noeuds venant en contact sur la surface de contact, la reconstruction à bon escient de la matrice pour les éléments de contact et une formulation en calcul parallèle

Cailletaud, Georges

Feyel, Frédéric

Yastrebov, Vladislav

HAL-MINES ParisTech

Traitement des proble`mes de contact en calcul paralle`leVladislav Yastrebov, Georges Cailletaud, Fre´de´ric FeyelTo cite this version:Vladislav Yastrebov, Georges Cailletaud, Fre´de´ric Feyel. Traitement des proble`mes de contacten calcul paralle`le. dir. M. Raous, P. Pasquet, C. Rey. Neuvie`me colloque en calcul desstructures, May 2009, Giens, France. pp.647-652, 2009. <hal-00399755>HAL Id: hal-00399755https://hal.archives-ouvertes.fr/hal-00399755Submitted on 29 Jun 2009HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Traitement des problèmes de contact en calcul parallèleV. A. Yastrebov1, G. Cailletaud1, F. Feyel21 MINES ParisTech, MAT, Centre des matériaux (CNRS UMR 7633)BP 87, 10 rue Henri Desbruères, 91003 Evry cedex{Vladislav.Yastrebov, Georges.Cailletaud}@mines-paristech.fr2 ONERA29, Avenue de la Division Leclerc, BP72, F92322 Chatillon Cedex - FranceFrederic.Feyel@onera.frRésumé — L’implantation numérique efficace des méthodes de traitement du contact reste encoreun sujet de recherche à part entière. On détaille ici la mise en œuvre d’une méthode de type Lagran-gien augmenté. Afin de gagner en performance, l’algorithme comporte une procédure particulièrede reconnaissance (séquentiel et parallèl) de nœuds venant en contact sur la surface de contact,la reconstruction à bon escient de la matrice pour les éléments de contact et une formulation encalcul parallèle.Mots clés — Contact-frottement, méthode de Lagrangien augmenté, reconnaissance de contact,parallélisation.1 IntroductionS’il y a maintenant une offre importante de méthodes numériques capables de traiter le pro-blème de contact en mécanique, l’implantation numérique de celles-ci reste encore une questiondélicate. La raison est que, au-delà de la méthode elle-même, il faut traiter les questions de rapiditéet de robustesse de la reconnaissance de contact, de stabilisation de la solution et de parallélisation.Les difficultés sont particulièrement évidentes lorsqu’on passe à des problèmes de grande taille,impossibles à traiter séquentiellement. L’approche cohérente qui est exposée dans cet article meten avant les traits principaux et les points délicats d’implantation de la méthode du Lagrangienaugmenté et d’une méthode de reconnaissance basée sur un maillage auxiliaire.2 Méthode du Lagrangien augmentéLa méthode de Lagrangien augmenté (LA) a été introduite par Hestenes [5] et Powell [7] pourles problèmes de programmation non linéaire sous contrainte. Alart et Curnier [1] et Pietrzak [6]ont bien adapté cette méthode pour le problème de contact avec friction dans le cadre des grandesdéformations. Le lecteur est renvoyé à l’article de Pietrzak pour un exposé détaillé de la méthode.On ne cite ici que les principales équations, appliquées à un milieu discrétisé par des élémentslinéaires en 2D.Pour chaque zone de contact entre deux corps, on définit un impacteur et une cible, l’impacteurcorrespondant aux nœuds susceptibles d’entrer en contact sur un corps, et la cible correspond àla surface de contact sur le deuxième corps ; on introduit donc une certaine asymétrie. Il y a deuxvaleurs importantes d’un point de vue cinématique : la distance normale entre les surfaces decontact et la vitesse de glissement tangentielle. Dans un milieu discrétisé, on définit la distanceentre chaque nœud de l’impacteur et sa projection sur la surface-cible (ou sur la prolongation decette surface dans certains cas) : gn = (x1−x2) ·n, où x1 est la coordonnée du nœud de l’impacteur,x2 la coordonnée d’un des nœuds de la cible et n une normale de la ligne reliant deux nœuds de lasurface-cible x2 et x3. La vitesse de glissement tangentielle est définie par Curnier [2] parg˙t(ξ,τ) = F(ξ,τ) ·∂Xp(ξ,τ)∂τ,où le τ est le temps, le F est le gradient de la transformation F = ∂x/∂X, Xp la coordonnée dela projection du nœud de l’impacteur sur la surface-cible dans la configuration de référence. Ilest d’usage de désigner par X les coordonnées des nœuds dans la configuration initiale, et x lescoordonnées dans la configuration actuelle.Après discrétisation en temps (on considère deux instants τ et τ + dτ), et en espace (on selimite à la considération d’un élément de contact composé d’un nœud de l’impacteur x1(τ) et deuxnœuds de la cible x2(τ), x3(τ), fig. 1), on a l’expression suivante :g˙t(τ,τ+dτ) =[γ(τ+dτ)γ(0)t(τ+dτ)⊗ t(0)]·[(β−α)t(0)γ(0)dτ]= γ(τ+dτ)t(τ+dτ)(β−α)dτ,où α =(x1(τ)−x2(τ)) · (x3(τ)−x2(τ))γ2(τ)et β =(x1(τ+dτ)−x2(τ+dτ)) · (x3(τ+dτ)−x2(τ+dτ))γ2(τ+dτ),γ(τ) = ‖x3(τ)−x2(τ)‖ et t = (x3−x2)/γ est le vecteur tangentiel normalisé.Figure 1 – Géométrie d’un élément de contactL’approche par Lagrangien augmenté est une procédure qui associe les méthodes des multi-plicateurs de Lagrange et de pénalisation, et qui s’exprime sous la forme :La(λ,x, pˆn) =Ncorps∑i=1U i(x)+Ncontact∑j=1Uˆ jc (λ,x, pˆn),où λ est un vecteur de multiplicateurs de Lagrange, pˆn = pn + rgn la pression normale de contactaugmentée par la valeur cinématique gn multipliée par r qui caractérise la mesure de la pénalité,U i l’énergie potentielle du ie corps du système et Uˆjc est l’énergie de contact qui correspond à laje zone de contact.Uˆ jc (x,λ, pˆn) =ZΓjcℓn(x,λn)dΓ+ZΓjcℓt(x,λt , pˆn)dΓ,Γjc désigne la surface de contact, ℓn et ℓt sont les fonctions régularisées suivantes :ℓn(x,λn) = gn(x)λn +r2‖gn(x)‖2−12rdist2R− [σn]ℓt(x,λn, pˆn) = g˙t(x) ·λt +r2‖g˙t(x)‖2−12rdist2Cˆ[pˆn][σt ] ,où σn = λn + rgn(x), σt = λt + rg˙t(x) et Cˆ[pˆn] est le cône de frottement de Coulomb, définissantla valeur maximale de la force tangentielle, qui dépend de la force normale augmentée pˆn.Une condition nécessaire d’existence de la solution est∇x,λLa =Ncorps∑i=1∇xUi(x)+Ncontact∑j=1∇x,λUˆjc (λ,x, pˆn) = G(x)+F(x,λ) = 0où G(x) correspond à la part sur la structure, et F(x,λ) correspond à la part dans le contact. Laforme de F, qui correspond aux forces nodales dans les éléments de contact estF(x,λ)=000−λr , si σn > 0 (non contact) ;σnn+ t⊗ t ·σt−(1−β)σnn−σt ·(gnγ t⊗n+(1−α)t⊗ t)−βσnn+σt ·(gnγ t⊗n−αt⊗ t)−1r(λ−σnn−σt), si ‖σt‖ ≤ −µσn; σn < 0 (adhérence)(n− εµt)σn(−(1−β)n+ εµ(gnγ n+(1−α)t))σn(−βn− εµ(gnγ n−αt))σn−1r(λ−σn(n− εµt)), si ‖σt‖>−µσn; σn < 0 (glissement)Il s’ensuit que la matrice tangente élémentaire est définie comme Ac = ∇F. Dans le cadred’un logiciel de calcul par éléments finis, il est nécessaire de définir les nœuds de l’impacteurqui entrent en contact et de les lier avec les nœuds de la cible. On forme pour cela des nouveauxéléments, les éléments de contact, en groupant un nœud de l’impacteur et quelques nœuds de lacible (au minimum deux nœuds dans le cas 2D pour les éléments linéaires). En dehors des nœudsqui forment cet élément, on ajoute des degrés de liberté correspondant aux multiplicateurs deLagrange, qui sont la force de contact normale λn et la force de contact tangentielle λt . Il est ensuitenécessaire de définir des forces nodales – un vecteur de vecteurs F(x,λ) –, et la matrice tangentede rigidité (matrice de tenseurs) Ac. Après chaque incrément de chargement, il faut reconstruireles éléments de contact, c’est-à-dire répéter la procédure de reconnaissance sur la zone de contact,qui s’est déformée. Cette partie de l’algorithme global de contact, utilisée pendant tout le cycle derésolution, est donc très importante : la vitesse de résolution dépend fortement de l’efficacité de laprocédure de reconnaissance, exposée dans la section suivante.3 Reconnaissance du contactLa question de reconnaissance du contact n’est pas triviale, et bon nombre de méthodes ontété proposées pour résoudre le problème. Dans cette section, on décrit une méthode inspirée de[3], qui s’est montrée efficace pour une large gamme de problèmes de contact. La méthode estbasée sur la division de l’espace de contact par un maillage régulier. Ce maillage est limité parun parallélépipède qui contient entièrement la surface-cible avec un jeu de valeur w défini parl’utilisateur. Ce paramètre définit la distance maximale sur laquelle on cherche les nœuds quientrent en contact. Pendant la procédure de reconnaissance, on cherche d’abord pour chaque nœudde l’impacteur le nœud le plus proche de la cible qui se trouve dans la même maille que le nœudde la cible ou dans les mailles voisines, à la suite de quoi on définit s’il y a une projection de cenœud sur les surfaces-cibles contenant le nœud trouvé.Il faut vérifier que le choix du paramètre w, défini par l’utilisateur, est tel que la surface-ciblese trouve entièrement dans l’union des cercles, de centres les nœuds de la cible, et de rayons ri =w(fig. 2a) : lcible ∈SNciblei=1 Ji, Ji = {r | ‖r− ri‖< w}⇔ w> L/2, où L = maxi=2,Ncible{‖ri− ri−1‖}.La procédure précise de reconnaissance comporte quatre étapes :1. Définition des limites du maillage de reconnaissance (de complexité Ncible)Il est nécessaire de déterminer les coordonnées du nœud de la cible le plus éloigné ; cetteprocédure se passe dans le repère global, c’est-à-dire qu’on cherche Rimin = min j=1,Ncible[xij]etRimax = max j=1,Ncible[xij], où i = 1,2 en 2D et i = 1,2,3 en 3D, le Ncible est le nombre de nœudsde la cible. Les limites de la boîte de recherche vont de {R1min−w,R2min−w,R3min−w} à {R1max +w,R2max+w,R3max+w} en 3D. L’introduction du jeu w est obligatoire, faute de quoi certains nœudsde l’impacteur peuvent être oubliés.La boîte se divise en N1×N2(×N3) mailles où Ni = max{[Rimax−Rimin2w], 1}, la notation [x] défi-nissant, dans tout cet article, la partie entière, si bien que les dimensions approximatives de chaquemaille sont 2w× 2w(×2w). Cette heuristique fournit une approche optimale en temps de recon-naissance.2. Remplissage des vecteurs (de complexité < Ncible +Nimp)Lorsque le maillage est formé, pour chaque nœud de la cible et chaque nœud de l’impacteurqui entrent dans la boîte de recherche, le numéro de maille où il se trouve doit être déterminé :N2D =[r1− r1min∆x]+nx[r2− r2min∆y], N3D =[r1− r1min∆x]+nx[r2− r2min∆y]+nxny[r3− r3min∆z]A chaque nouvelle opération, l’identificateur du nœud s’ajoute dans un vecteur de vecteurs Ai j, oùi est le numéro de la maille dans laquelle se trouve ce nœud, et j est le numéro d’ordre du nœudtrouvé dans cette maille. Cette procédure s’utilise séparément pour les nœuds de l’impacteur et dela cible, si bien qu’on obtient deux vecteurs de vecteurs I et T , qui correspondent respectivementà l’impacteur et à la cible.3. Définition des nœuds les plus proches (de complexité ∼4NimpNcible(nxny)2en 2D, ∼8NimpNcible(nxnynz)2en 3D)On calcule le minimum de la distance entre le nœud impacteur et les nœuds de la cible dansla maille courante et un sous-ensemble des mailles voisines. Ce sous-ensemble est construit en negardant que les mailles situées du «même côté» que le nœud impacteur (fig. 2b, 3 mailles sur 8en 2D et 7 sur 26 en 3D). Le choix de 2w pour la taille de maille assure qu’une recherche sur cesous-ensemble est suffisante.4. Définition d’une surface de projection (de complexité < 2Ncible en 2D)Chaque nœud est attaché à une surface (en 3D) ou une ligne (en 2D), il ne reste plus qu’àdéterminer sa projection sur la surface. Si aucune projection n’existe, il est nécessaire de vérifiersi le nœud de l’impacteur se trouve dans «l’angle mort» de la surface-cible discrétisée. C’estnotamment important lorsqu’il y a des conditions symétriques (fig. 2b).La complexité totale de la méthode de reconnaissance est ∼ 4Ncible +Nimp +4NimpNcible(nxny)2en 2Det ∼ (2+ k)Ncible +Nimp +8NimpNcible(nxnynz)2en 3D, où le k est le nombre moyen de surfaces attachéesà chaque nœud de la cible. La comparaison de la méthode proposée et de la méthode «brutale»sur le cas réel est représentée sur la figure 3a. Les avantages de la nouvelle méthode sont larobustesse et la rapidité. Cependant, l’efficacité de la méthode diminue si le maillage de la cibleou de l’impacteur est très irrégulier ou si la cible est plate.Figure 2 – a. Condition sur w ; b. Etapes de reconnaissance et exemples d’angle mortFigure 3 – a. Gain en temps en 2D de la méthode de reconnaissance par rapport à la méthode«brutale» ; b. variation de la distance d’adhérence en fonction du nombre de nœuds sur la surfacede contact4 ParallélisationLa parallélisation de la méthode de contact comporte deux parties : la parallélisation de la re-connaissance, et celle de la résolution. Elle est basée sur une méthode générique de décompositionde domaines [4] (contenant notamment la méthode FETI – Finite Element Tearing and Intercon-necting).On considère donc N processeurs, et une partition du modèle total en N sous-domaines.Dès que les nœuds qui entrent en contact sont trouvés, il est nécessaire d’ajouter au vol les nou-veaux éléments de contact aux sous-domaines qui contiennent les zones de contact. La versionparallèle de l’algorithme comprend 5 étapes :1. Chaque sous-domaine dont une partie du maillage contient une surface de contact l’envoieau premier sous-domaine.2. Le sous-domaine 0 définit les paramètres du maillage de reconnaissance : w, Rimin, Rimax, nx,etc. Il définit ensuite le numéro de maille pour chacun des nœuds de la surface de contact. Lemaillage de reconnaissance est alors divisé en N parties de même volume (fig. 4, zone A) chacuneconfiée à un sous-domaine. A chacune des parties envoyées, on ajoute en plus les nœuds de la ciblequi se trouvent dans les zones attenantes à la frontière de la division (fig. 4, zone B) ; l’épaisseurde cette zone est égale à une maille. Cela nous permet de ne plus échanger d’informations entreles sous-domaines pendant la suite de la reconnaissance.3. Chacun sous-domaine effectue indépendamment des autres la reconnaissance dans la partieproposée et envoie ses résultats au premier sous-domaine.4. Le premier sous-domaine définit alors les domaines où les éléments de contact sont à créer,et envoie ces informations aux sous-domaines concernés.5. On effectue ensuite la résolution parallèle avec la méthode générique de décomposition dedomaine.Figure 4 – Partition du maillage de reconnaissance pour la procédure parallèle5 ConclusionL’approche cohérente pour traiter le problème de contact en séquentiel et en parallèl est expo-sée. La nouvelle méthode de reconnaissance de contact est proposée en deux version : séquentielleet parallèle. La dernière permet d’effectuer la recherche par les sous domaines indépendamment.L’algorithme proposé dans cet article a été implanté dans Zébulon, un logiciel de calcul par élé-ments finis développé conjointement par l’ONERA, Northwest Numerics (Seattle, USA) et leCentre des Matériaux (CNRS UMR 7633). La méthode est validée sur le problème axisymétriquede contact entre une estampe rigide et un demi-espace en prenant en compte le frottement [8].L’erreur maximale pour les coefficients de frottement différents (µ = 0.2; 0.3; 0.4) ne dépasse pas10% pour une surface de contact contenant 400 nœuds (fig. 3b).Références[1] P. Alart, A. Curnier. A mixed formulation for frictional contact problems prone to Newtonlike solution method. Computer Methods in Applied Mechanics and Engineering, Vol 92,issue 3, pp. 353-375, 1991[2] A. Curnier, Q.-Ch. He and A. Klarbring. Continuum Mechanics Modeling of Large Defor-mation Contact with Friction. in : Contact Mechanics, M. Raous, M. Jean and J.-J. Moreau.(Eds.), Plenum, New York, 1995[3] W. Fujun, C. Jiangang, Y. Zhenhan. A contact searching algorithm for contact-impactor pro-blems. Acta Mechanica SINICA (English series), Vol 16, No 4, pp. 374-382, 2000[4] P. Gosselet, C. Rey. Non-overlapping domain decomposition methods in structural mecha-nics. Archives of Computational Methods in Engineering, Vol 13, No 4, pp 515-572, 2006[5] M. Hestenes, Multiplier and gradient method, J. Opt. Th. Appl., 4, pp. 303-320, 1969[6] G. Pietrzak, P. Alart, Large deformation frictional contact mechanics : continuum formula-tion and augmented Lagrangian treatment. Comput. Methods Appl. Mech. Engrg., 177, pp.351-381, 1999[7] B. T. Powell, A method for nonlinear constraints in minimisation problems, in : R. Fletcher,ed., Optimisation, Academic Press, London, 1969[8] D. A. Spence, The Hertz contact problem with finite friction. Journal of Elasticity, 5, No 3-4,pp. 298-319, 1975

Traitement des problèmes de contact en calcul parallèle

https://hal.archives-ouvertes.fr/hal-00399755

Traitement des problèmes de contact en calcul parallèle

Abstract

Similar works

Full text

Available Versions

HAL-MINES ParisTech