Implantation du modèle de plasticité cristalline

Le cadre de modélisation viscoplastique a été retenu dans ces travaux pour deux raisons : d’une part, comme révélé par la campagne expérimentale de caractérisation du chapitre 4, les matériaux étudiés présentent une sensibilité visqueuse qui, tout en étant faible n’est pas négligeable. D’autre part, la plasticité cristalline non visqueuse conduit à une indétermination concernant les systèmes de glissements activés au cours de la déformation. Les procédures mentionnées dans la littérature pour lever cette ambiguïté s’accompagnent d’un coût numérique parfois considérable qui se rajoute au temps de   calcul (déjà significatif) nécessaire pour déterminer l’état de contraintes aux points d’intégration. Dans l’optique d’une utilisation des outils développés sur des cas industriels, la viscoplasticité représente alors un bon compromis entre l’efficacité numérique et la précision des calculs. Des algorithmes explicites sont retenus pour l’intégration de la loi de comportement. En effet, les algorithmes implicites, basés sur une variable primaire (taux de glissement, gradient de déformation plastique…) dont l’estimation initiale est corrigée par des méthodes itératives de type Newton-Raphson, présentent un coût numérique plus élevé que celui des algorithmes explicites pour une précision comparable ((Ling et al., 2005), (Dumoulin et al., 2009)). De plus, dans le cadre d’applications à des procédés de mise en forme, les fortes non-linéarités géométriques et les conditions de contact fluctuantes entre les outils et la tôle nécéessitent forcément, indépendemment de la loi de comportement employée, de petits incréments de chargement. Cette caractéristique confère de facto un avantage certain aux algorithmes explicites en termes de compromis performance – précision – robustesse. 

Implantation dans ABAQUS/Explicit 

ABAQUS/Explicit et VUMAT

Le solveur dynamique explicite ABAQUS/Explicit a été conçu pour résoudre des problèmes dynamiques. Cependant, il est appliqué avec succès pour des problèmes quasistatiques et l’expérience montre qu’il se révéle souvent plus efficace que le solveur implicite ABAQUS/Standard. En effet, ABAQUS/Standard doit s’assurer que l’incrément de déplacement calculé vérifie au mieux l’équilibre mécanique global ; pour cela il effectue des corrections itératives sur la solution par un algorithme de type Newton-Raphson jusqu’à atteindre une précision satisfaisante. Il passe en outre par la formulation et l’inversion d’une matrice tangente, opération également très coûteuse en temps de calcul lorsque le modèle contient un grand nombre de noeuds. De plus, la convergence de l’algorithme de correction n’est pas assurée pour des conditions de contact complexes. Par contre, ABAQUS/Explicit détermine l’incrément de déplacement après intégration temporelle explicite de la relation d’équilibre dynamique : u¨ = M−1 · (F − I) (3.1)  avec u le vecteur des déplacements nodaux, M la matrice de masse, F le vecteur des efforts externes et I celui des efforts internes. L’équilibre mécanique n’est pas vérifié et aucune itération n’est faite ; les erreurs s’accumulant à chaque intégration, la stabilité est donc conditionnelle et il faut utiliser des pas de temps suffisamment réduits. Le pas de temps doit être inférieur au temps nécessaire à une onde de compression pour se propager à travers le plus petit élément du maillage, soit : ∆t ≤  L e cd  avec cd = s λ + 2µ ρ (3.2) où L e est une longueur caractéristique d’élément, cd est la vitesse de propagation d’une onde de compression, λ et µ sont les coefficients de Lamé du matériau et ρ sa densité. La possibilité fournie aux utilisateurs par ABAQUS/Explicit pour définir une loi de comportement est la routine VUMAT (Vectorized User MATerial). A partir des variables fournies en entrée par ABAQUS/Explicit à chaque incrément, elle a pour objectif de calculer le tenseur de contraintes de Cauchy au point d’intégration. Aussi, bien que la VUMAT oeuvre au point d’intégration, les réponses d’un certain nombre de points (NBLOCK) sont évaluées simultanément lors d’un seul appel de la routine. Cette structure vectorisée permet d’accroitre significativement l’efficacité numérique en répartissant le nombre de points évalués sur différents processeurs via les procédures de calcul parallèle. Le formalisme de la VUMAT comporte toutefois certaines particularités dont il faut tenir compte lors de l’implantation d’un modèle. En effet, dans le souci de garantir l’objectivité incrémentale des lois de comportement, le repère de travail employé par ABAQUS/Explicit est le référentiel tournant associé à la dérivée de Green-Naghdi dénommé référentiel en rotation propre. Comme détaillé au chapitre 2, ce référentiel est généré par la rotation pure R issue de la décomposition polaire du gradient de la transformation F = R · U où U est un tenseur d’élongation pure. Pour l’implantation d’une loi de plasticité cristalline au travers d’une VUMAT, on peut donc distinguer trois repères : – le reférentiel global qui est fixe dans l’espace ; – le référentiel en rotation propre qui suit la rotation moyenne de la matière au point matériel (au sens de Green-Naghdi), dans lequel le code de calcul travaille et exprime les tenseurs. Nous le désignerons par l’anglicisme MECS (Material Element Coordinate System) et un tenseur {•} exprimé dans ce référentiel s’écrira {•ˆ} ; Chapitre 3. Outils numériques 55 – le référentiel cristallographique défini par trois axes cristallographiques et qui suit la rotation du réseau cristallin. Il sera dénoté CACS (Crystal Associated Coordinate System). Un tenseur {•} exprimé dans le CACS sera noté {•e}. 

Algorithme en repère global

La loi de comportement telle que formulée en 2.3.1.3 est exprimée dans le repère global. Le tenseur des contraintes ainsi déterminé dans le repère global doit être retourné, au code, exprimé dans le MECS. Pour ce faire, le tenseur de rotation R qui transforme une base orthonormée directe du repère fixe en une base orthonormée directe du MECS est employé. ABAQUS/Explicit 6.11 fournit comme entrée de la routine VUMAT le gradient de la transformation en début (à t) et en fin d’incrément (à t + ∆t), notés respectivement Fn et Fn+1. De même, le tenseur d’élongation pure à droite U est fourni en début et en fin d’incrément, soit Un et Un+1. Soulignons ici que d’après la documentation (Hibbitt et al., 1992), le tenseur U fourni est bien exprimé dans le MECS. En effet, F et R sont des tenseurs mixtes, qui ont un pied dans les deux bases à savoir celle du repère global et celle du MECS. On écrit en terme de composantes : FiJ = RiKUKJ (3.3) où les indices en majuscule sont relatifs à la base de Green-Naghdi et ceux en minuscule à la base globale. Ainsi, les tenseurs de rotation en début d’incrément Rn et en fin d’incrément Rn+1 sont obtenus par la relation : R = F · U−1 (3.4) Les variables d’état sont stockées pour chaque point d’intégration en fin d’incrément dans un vecteur et récupérées au début de l’incrément suivant. Pour l’algorithme en repère global, on compte comme variables d’états : – F e le gradient élastique ; – F p le gradient plastique ; – τ s=1..12 c les cissions critiques d’activation des systèmes de glissement ; Chapitre 3. Outils numériques 56 – ρ s=1..12 les densités sur les systèmes de glissement si le modèle d’écrouissage de (Tabourot, 2001) est employé ; – γ le glissement total sur tous les systèmes ; – Q la matrice d’orientation. Pour garder une certaine lisibilité, les variables d’état intégrées numériquement (exception faite de Q) seront regroupées dans un pseudo-vecteur H = (F p , τ s c , ρs , γ). De manière similaire, nous introduisons l’ensemble Gn = ( ˙γ s )n des taux de glissement calculés à l’instant n. Précisons que dans le schéma ici adopté, les taux de glissement à l’instant n sont estimés explicitement à partir de données à l’instant n − 1. On peut donc écrire sous une forme globale : H˙ = h (t, H, G) (3.5) où h représente les différentes équations d’évolution des variables internes. L’intégration numérique des relations est effectuée par des schémas explicites de RungeKutta d’ordre 1 (méthode d’Euler explicite), 2 (méthode du point milieu) ou 4. On considère que l’indice n correspond à l’instant tn et que l’indice n + 1 renvoie à l’instant tn + ∆t où ∆t est l’intervalle de temps entre deux instants consécutifs. En employant la forme unifiée des trois méthodes de Runge-Kutta évoquées, on obtient pour un schéma d’ordre N : ∆Hi = h (tn + ai∆t, Hn + ai∆Hi−1, Gi) ∆t Hn+1 = Hn + PN i=1 bi∆Hi.