Illustration 2: Marche en temps avec adaptation de maillage
u2.08.09 Adaptation de maillage en non-linéaire#
Résumé :
Ce document présente la méthodologie de réalisation d’une étude non-linéaire avec adaptation de maillage. Cette adaptation du maillage au cours du transitoire est obtenue à l’aide du logiciel HOMARD appelé via la commande MACR_ADAP_MAIL [U7.03.01].
Un exemple illustre les possibilités et la mise en œuvre de stratégies de remaillage.
Adaptation en non-linéaire#
Dans cette partie, on décrit le problème type que l’on cherche à résoudre ainsi que les difficultés essentielles qu’il présente.
Problème considéré#
Le problème type auquel on s’intéresse est un transitoire impliquant un modèle avec une loi de comportement non-linéaire. On souhaite adapter le maillage support au cours de la résolution suivant le schéma de la figure . Le principe est le suivant :
on résout le transitoire sur un ou plusieurs pas de temps
on adapte le maillage sur la base du calcul d’un indicateur d’erreur
on reprend le transitoire
En fait, ce schéma dissimule une grande part de la complexité du problème. Prenons l’exemple de l’utilisation d’une loi de comportement élasto-plastique de von-Mises avec écrouissage isotrope. À chaque pas de temps, on stocke l’état du domaine de calcul \(\Omega\) , donné par le champ de déplacement, le champ de contraintes et le champ de variables internes \((u,\sigma ,\alpha )\) . C’est ce triplet porté par le maillage \(M\) qui détermine complètement l’équilibre de \(\Omega\) . Or dans le cas de l’utilisation de l’adaptation de maillage, le maillage devient un élément évolutif. Si on suppose que le maillage est adapté à chaque pas de temps, il faut aussi transférer le triplet \((u,\sigma ,\alpha )\) d’un maillage à l’autre. C’est la phase de projection. Malheureusement, une fois la projection réalisée, rien n’assure que la nouveau triplet assure l’équilibre de du domaine sur le nouveau maillage. Il faut pour cela corriger le nouveau triplet. C’est la phase de ré-équilibrage.
Phase de rééquilibrage#
Dans la suite, on utilisera les hypothèses suivantes :
utilisation d’une loi de comportement élasto-plastique de von-Mises avec écrouissage isotrope
adaptation du maillage à chaque pas de chargement
On pose les notations suivantes :
\({M}_{i}\) désigne le maillage au pas de chargement \(i\)
\(({u}_{i},{\sigma}_{i},{\alpha}_{i})\) désignent le triplet champ de déplacement, champ de contraintes et champ de variables internes au pas de chargement \(i\)
\({(x)}_{{M}_{i}}\) désigne la grandeur \(x\) portée par le maillage \({M}_{i}\)
On détaille maintenant le processus de calcul. On part d’un état équilibré \({({u}_{i},{\sigma}_{i},{\alpha}_{i})}_{{M}_{i}}\) porté par le maillage \({M}_{i}\) au pas de chargement \(i\) . On procède à la résolution du problème au pas de chargement \(i+1\) . On obtient l’état équilibré \({({u}_{i+1},{\sigma}_{i+1},{\alpha}_{i+1})}_{{M}_{i}}\) . Sur la base du calcul d’un indicateur d’erreur, on produit un nouveau maillage \({M}_{i+1}\) .
Pour continuer le calcul, il faut projeter l’état équilibré \({({u}_{i+1},{\sigma}_{i+1},{\alpha}_{i+1})}_{{M}_{i}}\) sur le nouveau maillage \({M}_{i+1}\) . Nous détaillerons cette projection dans la suite. On obtient alors l’état \({({u}_{i+1}^{\mathrm{proj}},{\sigma}_{i+1}^{\mathrm{proj}},{\alpha}_{i+1}^{\mathrm{proj}})}_{{M}_{i+1}}\) mais qui n’est pas forcément équilibré. En effet, la phase de projection peut induire des perturbations sur chacun des 3 champs et supprimer leur cohérence. Il s’agit donc de rééquilibrer l’état mécanique à l’instant \(i+1\) du domaine \(\Omega\) discrétisé par \({M}_{i+1}\) .
Pour ce faire, on se reporte à la documentation [R5.03.01] «Algorithme non-linéaire quasi-statique» et en ne prenant pas en compte les variables de commandes, on trouve que le système que l’on résout à la première itération de chaque pas de temps est le suivant :
Les seconds membres \({L}_{i+1}^{\text{méca}}-{Q}_{i}^{T}.{\sigma}_{i}\) et \({\mathrm{u}}_{i+1}^{d}-\mathrm{B}.{\mathrm{u}}_{i}\) mesurent respectivement :
l’écart à l’équilibre entre le chargement et un état mécanique de contraintes
l’écart à la vérification des conditions aux limites
La résolution du système non-linéaire sous-jacent amène à l’annulation de ces écarts et la production d’un état mécanique équilibré. C’est exactement ce que l’on cherche à faire non pas en partant de l’état mécanique de l’instant précédent \({({u}_{i},{\sigma}_{i},{\alpha}_{i})}_{{M}_{i}}\) comme dans () mais à partir de l’état mécanique projeté \({({u}_{i+1}^{\mathit{proj}},{\sigma}_{i+1}^{\mathit{proj}},{\alpha}_{i+1}^{\mathit{proj}})}_{{M}_{i+1}}\) .
Si on le transpose au cas qui nous intéresse, le système à résoudre est donc le suivant :
Du point de vue de la modélisation, cela revient donc à relancer la résolution du pas de chargement \(i+1\) en utilisant comme état initial l’état mécanique projeté \({({u}_{i+1}^{\mathrm{proj}},{\sigma}_{i+1}^{\mathrm{proj}},{\alpha}_{i+1}^{\mathrm{proj}})}_{{M}_{i+1}}\) .
L’utilisation du l’état mécanique projeté \({({u}_{i+1}^{\mathrm{proj}},{\sigma}_{i+1}^{\mathrm{proj}},{\alpha}_{i+1}^{\mathrm{proj}})}_{{M}_{i+1}}\) comme état initial est capitale. Elle permet d’abord de partir d’un point de départ le plus proche possible du point d’arrivée et donc d’augmenter les chances de voir converger la méthode de Newton. Mais il y a une raison bien plus importante. Dans le cadre d’utilisation de lois de comportement irréversible, les variables internes, telle que la déformation plastique cumulée, portent l’histoire mécanique de chaque point de la structure. Il faut partir d’un état de variables internes le plus fidèle possible à l’état réel de la structure pour bien traduire son histoire mécanique.
Cela nous amène à une nouvelle constatation: il est essentiel de produire un état mécanique projeté \({({u}_{i+1}^{\mathit{proj}},{\sigma}_{i+1}^{\mathit{proj}},{\alpha}_{i+1}^{\mathit{proj}})}_{{M}_{i+1}}\) le plus fidèle possible. C’est ce que nous allons aborder maintenant.
Projection#
La projection des champs se fait avec l’opérateur PROJ_CHAMP; pour toute information sur sa syntaxe et ses fonctionnalités complètes, le lecteur se reportera à [U4.72.05].
Nous nous intéressons à la projection de l’état mécanique qui est composé de 3 champs. Nous allons examiner la projection de chacun de ces champs.
Projection du champ de déplacement#
Le champ de déplacement DEPL est un champ aux nœuds continu. Sa projection ne pose pas de problème particulier. La méthode de collocation, qui est utilisée par défaut pour ce type de champ, est parfaitement adaptée.
Projection du champ de contraintes#
Le champ de contraintes (et/ou d’efforts) SIEF_ELGA est nativement calculé aux points d’intégration. De manière à minimiser la perturbation induite par la projection, il est préférable d’exploiter directement le champ au point d’intégration. C’est la méthode que l’opérateur de projection choisira automatiquement pour ce champ.
Il est possible de faire la projection en utilisant un champ par éléments aux nœuds ou un champ aux nœuds. Cela implique des manipulations supplémentaires assez techniques. Cette approche est déconseillée.
Projection du champ de variables internes#
Le champ de variables internes VARI_ELGA est nativement calculé aux points d’intégration. De manière à minimiser la perturbation induite par la projection, il est préférable d’exploiter directement le champ au point d’intégration. C’est la méthode que l’opérateur de projection choisira automatiquement pour ce champ.
Comme pour le champ de contraintes, il est possible d’utiliser un champ par éléments aux nœuds ou un champ aux nœuds mais cela pose plus de problème. En effet, si l’on considère le cas de la plasticité de Von Mises, la déformation plastique cumulée \(p\) doit toujours être positive, \(p\ge 0\) . Considérons la figure . Elle représente un élément SEG2 à 2 nœuds \(\mathit{N1}\) et \(\mathit{N2}\) et à 2 points d’intégration \(\mathit{PG1}\) et \(\mathit{PG2}\) . Cet élément porte un champ à ses points d’intégration de valeur \(\mathit{VG1}\) et \(\mathit{VG2}\) , qui sont toutes deux positives. Or, suite à l’extrapolation aux sommets, on constate que la valeur \(\mathit{VN1}\) au sommet \(\mathit{N1}\) reste bien positive tandis que la valeur \(\mathit{VN2}\) au sommet \(\mathit{N2}\) devient négative. De manière à être utilisée comme valeur initiale, la déformation plastique cumulée \(p\) doit être positive ou nulle; la valeur \(\mathit{VN2}\) sera donc mise à zéro.
L’utilisation d’un champ de variables internes aux nœuds est fortement déconseillée car, en plus des perturbations inhérentes à la projection comme au paragraphe précédent, elle impose un traitement supplémentaire pour rendre le champ admissible qui le perturbe notablement.
Illustration 3: Perturbation liée à l’extrapolation
Mise en œuvre#
Schéma général#
Le schéma général de l’adaptation est présenté sur la figure . Il comprend les 4 phases décrites précédemment:
résolution sur un ou plusieurs pas de chargement
adaptation du maillage sur la base du calcul d’un indicateur d’erreur
projection des champs à la fin du dernier pas de temps sur le nouveau maillage
équilibrage des champs à la fin du dernier pas de temps sur le nouveau maillage
Illustration 4: Processus complet d’adaptation
Exemple#
Un exemple de mise en œuvre pratique est présenté dans le cas-test Code_Aster SSNP158 [V6.03.158].
Références#
Les documents [bib1] à [bib4] traitent de l’outil d’adaptation de maillage HOMARD.
Les documents [bib5] à [bib9] traitent des estimateurs d’erreur.
NICOLAS & al. https://www.code-aster.org/outils/homard
NICOLAS, T. FOUQUET, “Logiciel HOMARD - Volume 1 – Présentation générale”, rapport EDF H-I23-2008-04107
«Macro-commande MACR_ADAP_MAIL», [U7.03.01]
«Macro-commande MACR_INFO_MAIL», [U7.03.02]
«Estimateur d’erreur de Zhu-Zienkiewicz en élasticité 2D», [R4.10.01]
«Estimateur d’erreur en résidu», [R4.10.02]
«Indicateur d’erreur spatiale en résidu pour la thermique transitoire», [R4.10.03]
«Détection des singularités et calcul d’une carte de taille», [R4.10.04]
DELMAS: «Estimateurs d’erreur en quantités d’intérêt», Compte-rendu AMA 06-66