r5.03.01 Algorithme non linéaire quasi-statique (STAT_NON_LINE)#
Résumé :
L’opérateur STAT_NON_LINE [U4.51.03] permet dans le cas d’une sollicitation quasi-statique d’intégrer divers types de non-linéarités provenant du comportement du matériau, de grandes transformations géométriques ou des conditions de contact/frottement. On décrit ici l’algorithme de résolution global employé.
L’intégration des relations de comportement proprement dite est décrite dans les documents [R5.03….] et [R7.01….] , (par exemple [R5.03.02] pour l’élasto-plasticité), auxquels on pourra se reporter pour plus de détails.
Pour les calculs en grandes transformations géométriques, on pourra consulter par exemple le document [R5.03.20] sur l’élasticité non linéaire en grands déplacements, ou les documents [R5.03.21], [R5.03.22] sur la thermoélastoplasticité à écrouissage isotrope.
Pour le contact frottement, il existe trois documents: [R5.03.50] sur la formulation discrète du contact/frottement, [R5.03.52] pour la formulation hybride par des éléments de contact/frottement.
Pour tout ce qui concerne le pilotage, il faut se référer au document [R5.03.80].
Méthode de Newton#
Principe de la méthode#
La méthode de Newton est une méthode classique de résolution des équations du type recherche de zéro. Considérons une fonction vectorielle \(F\) non-linéaire du vecteur \(x\) . On cherche le zéro de cette fonction, c’est à dire:
La méthode de Newton consiste à construire une suite de vecteurs \({x}^{n}\) convergeant vers la solution \(x\) . Pour trouver le nouvel itéré \({x}^{n+1}\) , on approche \(F({x}^{n+1})\) par un développement limité à l’ordre un autour de \({x}^{n}\) et on exprime que \(F({x}^{n+1})\) doit être nul:
Soit encore:
Finalement:
\({F}^{'}(x)\) est l’application linéaire tangente associée à la fonction \(F\) . La dérivée au point \(x\) dans la direction \(h\) est définie comme la dérivée directionnelle suivante:
La matrice de \({F}^{'}(x)\) dans les bases choisies pour les espaces vectoriels concernés s’appelle la matrice jacobienne de \(F\) au point \(x\) . Lorsque \(F\) est une fonction d’un espace vectoriel euclidien à valeurs réelles, \({F}^{'}(x)\) est une forme linéaire, et on peut montrer qu’il existe un vecteur (unique), noté \(\nabla F(x)\) et appelé le gradient de \(F\) , tel que :
C’est-à-dire le produit scalaire de \(h\) et du gradient de \(F\) .
Lorsque l’on est proche de la solution, la convergence de la méthode de Newton est quadratique c’est-à-dire que le nombre de zéros après la virgule dans l’erreur double à chaque itération (0,19 – 0,036 – 0,0013 – 0,0000017 par exemple). Mais cette méthode (utilisant la vraie tangente) a plusieurs inconvénients :
Elle nécessite le calcul de la tangente à chaque itération, ce qui est d’autant plus coûteux que la taille du problème est grande (surtout dans le cas où on utilise un solveur direct),
Si l’incrément est grand, la tangente (dite cohérente ou consistante) peut conduire à des divergences de l’algorithme,
Elle peut ne pas être symétrique, ce qui oblige à utiliser des solveurs particuliers.
C’est pour cette raison que l’on peut utiliser d’autres matrices à la place de la matrice tangente: la matrice élastique, une matrice tangente obtenue antérieurement, la matrice tangente symétrisée,…
Adaptation de la méthode de Newton au problème posé#
Dans un premier temps, on ne prend pas en compte les conditions aux limites de Dirichlet. On doit résoudre un système (non-linéaire car dépendant de \({u}_{i}\) ) de la forme:
Les forces internes \({L}_{i}^{int}\) peuvent symboliquement être notées \({Q}_{i}^{T}.{\sigma}_{i}\) , où \({Q}_{i}^{T}\) est la matrice associée à l’opérateur divergence (expression du travail du champ de déformations virtuelles ). Les forces internes s’expriment alors:
Et les forces du chargement mécanique:
Où:
\({w}_{i}\) désigne le champ des déplacements virtuels;
\({f}_{i}\) désigne les forces volumiques s’appliquant à l’instant \({t}_{i}\) sur \(\Omega\) ;
\({t}_{i}\) désigne les forces surfaciques s’appliquant à l’instant \({t}_{i}\) sur la frontière \(\Gamma\) de \(\Omega\) .
L’application de la méthode de Newton conduit à résoudre une suite de problèmes linéaire du type (\(n\) est le numéro de l’itération de Newton, \(i\) celui de la variable d’instant):
On note \(\delta {u}_{i}^{n+1}={u}_{i}^{n+1}-{u}_{i}^{n}\) l’incrément de déplacement entre deux itérations de Newton successives. La matrice \({K}_{i}^{n}\) est la matrice de rigidité tangente en \({u}_{i}^{n}\) et le vecteur \({L}_{i}^{int,n}\) représente les forces internes à la \({n}^{\text{ième}}\) itération de Newton du \({i}^{\text{ème}}\) pas de charge. La quantité \({R}_{i}^{n}=({L}_{i}^{\text{méca},n}-{L}_{i}^{int,n})\) représente les forces non équilibrées, que l’on appelle le «résidu d’équilibre». La matrice \({K}_{i}^{n}\) est la matrice de l’application linéaire tangente de la fonction \({R}_{i}^{n}\) , elle vaut donc:
En l’absence de forces suiveuses [§ 3.2 ], le second terme est nul. Il ne reste donc de la matrice \({K}_{i}^{n}\) que la dérivée au point \({u}_{i}^{n}\) des forces internes par rapport aux déplacements :
Une petite erreur dans l’évaluation des forces internes peut avoir des conséquences graves, car c’est leur calcul exact qui garantit, si l’on converge, que ce sera vers la solution cherchée. Par contre, il n’est pas toujours nécessaire d’utiliser la vraie matrice tangente, dont le calcul et la factorisation sont coûteux. Par exemple, une variante de la méthode utilise la matrice élastique \({K}_{\text{élas}}\) . La méthode utilisant la vraie matrice tangente \({K}_{i}^{n}\) (dite aussi matrice cohérente ou consistante) s’appelle la méthode de Newton ; les méthodes utilisant d’autres matrices (comme par exemple la matrice élastique \({K}_{\text{élas}}\) ) sont appelées méthodes de Newton modifiées ou méthodes de quasi-Newton. Le choix entre une matrice tangente (la dernière obtenue ou une matrice précédente) et une matrice élastique est effectué dans Code_Aster par l’intermédiaire du mot-clé MATRICE= “TANGENTE” ou MATRICE=”ELASTIQUE” du mot-clé facteur NEWTON. De plus, il est possible d’utiliser une matrice de décharge, c’est-à-dire d’une matrice à variables internes constantes (l’évolution des non linéarités n’est donc pas prise en compte dans cette matrice), en dessous d’un certain pas de temps, pour certaines lois de comportement. On se reportera à la documentation [U4.51.03] pour l’utilisation de cette fonctionnalité.
La méthode de Newton à matrice tangente consistante peut s’illustrer simplement à l’aide du schéma de la Figure 2.2.1-a.
Fig. 192 Figure 2.2.1-a#
On enchaîne avec la boucle des itérations de Newton qui permet, à la convergence, d’obtenir les valeurs de \(\Delta {\mathrm{u}}_{i}\) , et donc celles de \({\mathrm{u}}_{i}\) par application de l’équation ():
Phase de prédiction d’Euler#
L’expérience montre que la convergence de la méthode de Newton est fortement dépendante d’un choix judicieux de l’estimation initiale : «plus l’estimation initiale est proche de la solution, plus l’algorithme converge vite». Pour amorcer le processus itératif de la méthode, il est donc utile de déterminer un «bon» incrément initial \(\Delta {\mathrm{u}}_{i}^{0}\) . Pour cela, on linéarise par rapport au temps le problème continu: c’est ce qu’on appelle la phase de prédiction .
Linéarisation#
On va donc linéariser le système () par rapport au temps autour de \({\mathrm{u}}_{i-1}\) . On commence par linéariser les forces internes \({L}_{i}^{int}\) :
On suppose que le chargement mécanique ne dépend pas du temps (les charges suiveuses sont exclues), donc :
En ré-injectant et dans la première équation du système , on obtient pour l’équation d’équilibre:
On a équilibre à l’instant \({t}_{i-1}\) , c’est-à-dire:
Et il reste donc:
On obtient le système d’équations permettant de calculer des valeurs prédictives \(\Delta {\mathrm{u}}_{i}^{0}\) :
Matrice tangente de prédiction#
La quantité \({\frac{\partial {\mathrm{L}}^{int}}{\partial \mathrm{u}}\mid }_{{\mathrm{u}}_{i-1}}\) désigne la dérivée partielle à temps constant de \({L}_{i-1}^{int}\) , elle peut se développer:
La matrice \({K}_{i-1}\) est appelée matrice tangente de prédiction. La dépendance de la matrice \(Q\) par rapport aux déplacements est négligée dans l’hypothèse des petits déplacements : le terme \({\frac{\partial {Q}^{T}}{\partial u}\mid }_{{u}_{i-1}}\) , dit terme de rigidité géométrique, disparaît donc de . Ce terme est pris en compte pour les grandes transformations (voir § 1.2.2 ). Pour les développeurs, précisons que le calcul de la matrice tangente lors de la phase de prédiction est effectué via l’option de calcul RIGI_MECA_TANG.
Variantes de la prédiction#
Il existe d’autres options de prédiction disponibles dans STAT_NON_LINE.
On peut utiliser une matrice élastique \({K}_{\text{élas}}\) à la place de la matrice tangente en vitesse \({K}_{i-1}\) , c’est l’option PREDICTION=”ELASTIQUE” (option RIGI_MECA).
On peut utiliser une matrice sécante \({K}_{\text{sécante}}\) à la place de la matrice tangente en vitesse \({K}_{i-1}\) , c’est l’option PREDICTION=”SECANTE”(option RIGI_MECA_ELAS). La matrice sécante est une matrice élastique dont le module de Young est utilisé en appliquant l’endommagement (voir par exemple [R5.03.18] pour plus de détails)
On peut utiliser un incrément de déplacement précédemment calculé à la place de l’estimation, c’est l’option PREDICTION=”DEPL_CALCULE”. Dans ce cas on ne fait aucune inversion de système et le \(\Delta {u}_{i}^{0}\) est directement donné. Voir documentation [U4.51.03] pour son usage.
On peut utiliser un incrément de déplacement extrapolé par rapport au pas précédent. On calcule l’estimation de l’incrément de déplacement à partir de l’incrément total obtenu comme solution au pas de temps précédent (pondéré par le rapport des pas de temps). C’est l’option PREDICTION=”EXTRAPOL”.
Dans ces deux derniers cas, afin d’assurer que le déplacement initial est cinématiquement admissible, on projette l’estimation sur l’ensemble des champs cinématiquement admissibles (i.e. satisfaisant les conditions aux limites de Dirichlet) selon la norme donnée par la matrice élastique , qui doit donc être calculée.
Vecteur second membre des variables de commande#
Une variable de commande \(\beta (t)\) est une quantité scalaire, fonction du temps et de l’espace [1]_ , donnée a priori par l’utilisateur via le mot-clef AFFE_VARC dans l’opérateur AFFE_MATERIAU. C’est un paramètre du problème et non une inconnue . La quantité \({\frac{\delta {L}^{int}}{\delta t}\mid }_{{t}_{i-1}}\) désigne la différentielle partielle, par rapport à \(t\) et à \(u\) constant, de \({L}^{int}=Q.\sigma (t,\beta (t))\) . Cette notation particulière a pour but d’attirer l’attention sur le fait que pour \({n}_{\mathrm{varc}}\) variables de commandes, on écrit la différentielle totale:
Si on prend comme exemple la variable de commande décrivant la température \(\theta\) :
On suppose que la température varie linéairement entre les deux instants:
La dépendance de \(\sigma\) par rapport au temps et par rapport à la température permet d’écrire:
Car \(Q\) ne dépend pas du temps et donc \(\frac{\delta}{\delta t}({Q}^{T})=0\) . Le vecteur \(\Delta {L}_{i}^{\text{varc}}\) , dont l’expression est donnée par , est l’incrément de chargement de température (attention au changement de signe!) que l’on a généralisé à toutes les variables de commandes : température, irradiation, phases métallurgiques (voir [R4.04.02] ),…
On ne tient pas compte actuellement de la dépendance explicite des contraintes par rapport au temps et donc le premier terme de vaut zéro. Et donc finalement :
L’incrément de chargement des variables de commande \(\Delta {L}_{i}^{\text{varc}}\) , issu de la dérivation des forces internes par rapport aux variables de commande est une estimation de l’effet d’une variation des variables de commandes.
Dans le cas de la température, si l’on note \(K\) le module de compression hydrostatique et \(\alpha\) le coefficient de dilatation thermique, la contrainte thermique s’écrit:
Où \(I\) est la matrice identité. Et donc, en appliquant :
Dans le cas élastique, ce sont les forces internes associées à une dilatation thermique (ce n’est pas à proprement parler un chargement, cela s’assimile plutôt à l’effet d’une déformation initiale). Cette estimation est utilisée dans la phase de prédiction et dans le critère d’arrêt. Si les dilatations thermiques font sortir la structure du domaine élastique (plasticité par exemple), cette estimation sera corrigée lors des itérations de Newton.
Vecteur second membre du chargement mécanique#
L’incrément de chargement mécanique \(\Delta {\mathrm{L}}_{i}^{\text{méca}}\) est composé des charges mortes (indépendantes de la géométrie, comme la pesanteur) et des charges suiveuses. En réalité, il existe des cas (le premier incrément de charge, par exemple) où \({L}_{i-1}^{\text{méca}}\) est inconnu. On rappelle que l’incrément de chargement () s’écrit:
On a équilibre à l’instant \({t}_{i-1}\) , donc en appliquant :
L’expression des forces internes au pas de temps précédent \({L}_{i-1}^{int}\) implique soit de sauvegarder ce vecteur du calcul précédent s’il existe (reprise d’un calcul antérieur), soit d’intégrer la loi de comportement à partir de l’état initial donné par l’utilisateur (ce qui peut être coûteux). Par souci de simplicité et d’efficacité, on choisit de ne pas réintégrer la loi de comportement et on exprime simplement les forces internes comme les forces nodales en prenant les contraintes connues à cet instant, soit:
D’où la nouvelle expression:
Le calcul direct à partir de () demande à l’utilisateur de prendre garde à la cohérence entre le champ des contraintes, les champs de déplacements et de variables internes ( DEPL, SIGMet VARIdans ETAT_INIT). vis-à-vis de l’intégration de la loi de comportement dans le cas d’une reprise de calcul.
Système linéaire#
En ré-injectant l’expression de \({\frac{\partial {L}^{int}}{\partial u}\mid }_{{u}_{i-1}}\) (équation ) , de \({\frac{\delta {L}^{int}}{\delta t}\mid }_{{t}_{i-1}}\) (équation ) et \(\Delta {L}_{i}^{\text{méca}}\) (équation ) dans , le système d’équations en prédiction s’écrit:
Phase de correction de Newton#
A l’issue de la phase de prédiction, nous nous retrouvons avec une estimation de l’incrément des déplacements \(\Delta {u}_{i}^{0}\) . Si cette estimation est exacte (modulo l’application des critères de convergence décrits au § 2.5 ), alors on obtient la solution convergée pour le pas de temps \({t}_{i}\) :
Mais si ce n’est pas le cas, on doit trouver la valeur \(\Delta {\mathrm{u}}_{i}\) de l’incrément de déplacement à partir de la valeur \({\mathrm{u}}_{i-1}\) obtenues à l’équilibre précédent (instant \({t}_{i-1}\) ):
Si la phase de prédiction a convergé, on a donc trivialement:
Sinon, on prend comme valeur initiale \(\Delta {\mathrm{u}}_{i}^{0}\) obtenue à l’issue de la phase de prédiction, avant de corriger par les itérations \(\delta {\mathrm{u}}_{i}^{n}\) de la méthode de Newton. Avec un nombre \({n}_{\mathrm{CV}}\) suffisant d’itérations de Newton (toujours avec l’arbitrage du critère de convergence):
Tant qu’on a pas convergé (si le nombre d’itérations de Newton n’est pas suffisant), on note:
Le déplacement total, pour le pas de temps \(i\) et l’itération de Newton \(n\) s’écrira donc:
À chaque itération, on doit résoudre un système permettant de déterminer \(\delta {\mathrm{u}}_{i}^{n}\) , les incréments des déplacements depuis le résultat \({\mathrm{u}}_{i}^{n-1}\) de l’itération précédente:
On note aussi:
Soit encore:
Linéarisation#
On doit linéariser le système () par rapport aux inconnues en \({\mathrm{u}}_{i}^{n}\) à \({t}_{i}\) constant. On commence par linéariser les forces internes \({L}_{i}^{int,n}\) :
On suppose que le chargement mécanique ne dépend pas du temps (les charges suiveuses sont exclues).
Système linéaire#
En ré-injectant () dans l’équation (), on obtient pour l’équation d’équilibre:
La quantité \({\frac{\partial {L}^{int}}{\partial u}\mid }_{{u}_{i}^{n-1}}\) est appelée matrice tangente cohérente et elle est notée \({K}_{i}^{n-1}\) :
Le système à résoudre s’écrit finalement:
Le vecteur des forces internes \({L}_{i}^{int,n-1}\) est calculé à partir des contraintes \({\sigma}_{i}^{n-1}\) . Celles-ci étant calculées à partir des déplacements \({\mathrm{u}}_{i}^{n-1}\) par l’intermédiaire de la relation de comportement du matériau [§ 1.1 ]. En fait, dans le cas des comportements incrémentaux , \({\sigma}_{i}^{n-1}\) est calculé à partir de \(({\sigma}_{i-1},{\alpha}_{i-1})\) et de l’incrément de déformation \(\epsilon (\Delta {\mathrm{u}}_{i}^{n-1})\) induit par l’incrément de déplacement depuis le début du processus itératif (y compris la phase de prédiction) ou par le gradient de la transformation \(\mathrm{F}\) dans le cas des grandes transformations.
Différence des matrices en prédiction et correction#
Il est important de souligner que la matrice tangente issue de l’option RIGI_MECA_TANG (phase de prédiction) et la matrice tangente issue de l’option FULL_MECA (phase de correction) ne sont en général pas identiques. Supposons que l’on a atteint la convergence pour l’instant \({t}_{i-1}\) et que l’on cherche maintenant à obtenir l’équilibre pour l’instant suivant \({t}_{i}\) . La matrice issue de RIGI_MECA_TANG provient d’une linéarisation des équations d’équilibre par rapport au temps autour de \(({u}_{i-1},{\lambda}_{i-1})\) i.e. autour de l’équilibre à l’instant \({t}_{i-1}\) . C’est donc la matrice tangente du système convergé à l’instant \({t}_{i-1}\) .
Par contre, la matrice issue de FULL_MECA provient d’une linéarisation des équations d’équilibre par rapport au déplacement autour de \({\mathrm{u}}_{i}^{n}\) i.e. autour de l’équilibre à l’instant \({t}_{i}\) .
On peut interpréter les différences entre RIGI_MECA_TANG et FULL_MECA en d’autres termes. On peut ainsi montrer que la matrice issue de RIGI_MECA_TANG correspond à l’opérateur tangent du problème continu en temps, dit aussi problème en vitesse (et relie la vitesse de contrainte à la vitesse de déformation), alors que la matrice issue de FULL_MECA correspond à l’opérateur tangent du problème discrétisé en temps. Le document [R5.03.02] donne l’expression dans chacun des deux cas pour la relation d’élasto-plasticité de Von Mises à écrouissage isotrope ou cinématique linéaire.
On rappelle que le traitement d’une relation de comportement [R5.03.02 § 5] consiste à:
Calculer les contraintes \({\sigma}_{i}^{n}\) et les variables internes \({\alpha}_{i}^{n}\) à partir de l’état initial \(({\sigma}_{i-1},{\alpha}_{i-1})\) et de l’incrément de déformation \(\varepsilon (\Delta {u}_{i}^{n-1})\) induit par l’incrément de déplacement depuis le début du processus itératif (y compris la phase de prédiction).
Calculer les forces internes \({L}_{i}^{int,n}={Q}_{i}^{T}.{\sigma}_{i}^{n}\) .
Calculer (éventuellement) la matrice tangente (option RIGI_MECA_TANG pour la phase de prédiction, option FULL_MECA pour les itérations de Newton).
Critères de convergence#
À la fin de la prédiction et à chaque itération de Newton, on doit estimer si le processus itératif a convergé (l’équilibre de la structure est atteint). On se place au pas de temps courant \({t}_{i}\) et à l’itération de Newton \(n\) (étant entendu que la valeur \(n=0\) correspond à la prédiction). Il existe quatre critères de convergence.
Le critère RESI_GLOB_MAXI consiste à vérifier que la norme infinie [2] du résidu est inférieure à la valeur \(\gamma\) spécifiée par l’utilisateur.
Il n’est pas conseillé d’utiliser ce critère seul, car on ne peut pas facilement avoir une idée des ordres de grandeurs absolus admissibles.
Le critère RESI_GLOB_RELA (choisi par défaut) revient à vérifier que le résidu est suffisamment petit, comme précédemment, et ceci relativement à une quantité représentative du chargement.
\(\eta\) étant la précision relative souhaitée donnée par l’utilisateur sous le mot-clef RESI_GLOB_RELA(ou la valeur par défaut de \({10}^{\text{-6}}\) ).
On peut remarquer que, dans le cas d’utilisation de RESI_GLOB_RELA, le critère peut devenir singulier si le chargement extérieur \({\mathrm{L}}_{i}^{\text{méca}}+{\mathrm{L}}_{i}^{\text{varc}}\) devient nul. Ceci peut arriver en cas de décharge totale de la structure. Si un tel cas de figure se présente (i.e. chargement \({10}^{\text{‑6}}\) fois plus petit que le plus petit chargement observé jusqu’au présent incrément), le code utilise alors le critère RESI_GLOB_MAXI avec comme valeur celle observée à la convergence de l’incrément précédent. Lorsque le chargement redevient non nul, on revient au critère initial.
Le troisième critère est le critère RESI_REFE_RELA: l’idée de ce critère est de construire une force nodale de référence, qui servira à estimer terme à terme, la nullité (approchée) du résidu:
Plus précisément, la force nodale de référence \({\mathrm{F}}_{j}^{\text{ref}}\) est construite à partir de la donnée d’une amplitude de référence \({A}^{\text{ref}}\) qui peut être:
Une contrainte;
Une pression, une température dans le cas de la THM;
Une force généralisée dans le cas des poutres ou des coques;
autres …
La liste est accessible dans la documentation d’utilisation de la commande STAT_NON_LINE [U4.51.03], description de l’opérande RESI_REFE_RELA.
Si on prend comme exemple une contrainte, l’amplitude de référence \({\sigma}^{\text{ref}}\) étant donnée par l’utilisateur via le mot-clef SIGM_REFE . A partir de cette amplitude de contrainte de référence, on définit le tenseur \({\sigma}^{\text{test}}\) : il est nul pour toutes ces composantes, sauf la j -ième qui vaut \({\sigma}^{\text{ref}}\) . On définit alors, pour chaque nœud de chaque élément la force nodale \({\tilde{R}}_{i}^{e}\) (le but étant de donner une idée de l’importance d’une composante en un point de Gauss de la contrainte sur la force nodale) :
Avec \(N\) le nombre de points de Gauss de l’élément, \(M\) le nombre de composantes du tenseur des contraintes ; l’exposant \(\alpha\) servant à noter l’évaluation de quantité au point de Gauss. Par exemple \({\omega}^{\alpha}\) sont les poids des points de Gauss.
Remarque: Pour certains éléments, comme les barres, les grilles ou les membranes, cette définition aboutit à des résidus de référence nuls sur certains axes. Pour y remédier, on détermine les forces de référence nodales via un calcul d’ordre de grandeur basé sur la taille de l’élément.
La force nodale de référence est alors définie par:
où \({\Gamma}_{i}\) est l’ensemble des éléments connectés au nœud \(i\) .
Le quatrième critère est le critère RESI_COMP_RELA: l’idée de ce critère est de séparer les différentes contributions du résidu composantes par composantes (au sens DX, DY, DZ, PRE1, PRE2, TEMP). Chacun des vecteurs obtenus sera ensuite normé par la force interne correspondant à ce résidu. Ce choix de critère de convergence n’a de sens que pour les problèmes fortement évolutifs, typiquement les problèmes THM. Qui plus est, ce choix sera efficace dés lors que l’on aura des problèmes à forts contrastes. Ce sont en effet les zones à «fort gradient» qui vont contrôler la convergence. On définit \({F}^{n}(c)\) la partie du résidu \({\mathrm{Q}}^{T}.{\sigma}_{i}^{n}-{\mathrm{L}}_{i}^{\text{méca}}\) correspondant à la composante \(c\) et \({\mathrm{L}}^{int,0}(c)\) le vecteur des forces internes en prédiction correspondant à ce même composant \(c\) . Le critère RESI_COMP_RELA revient alors à vérifier que ce résidu est suffisamment petit, c’est-à-dire :
La convergence est décrétée réalisée lorsque tous les critères spécifiés par l’utilisateur sont vérifiés simultanément. Par défaut, on fait un test sur le résidu global relatif (RESI_GLOB_RELA) et le nombre maximum d’itérations de Newton (ITER_GLOB_MAXI).
Pour les résidus RESI_GLOB_RELA et RESI_GLOB_MAXI , toutes les composantes du champ de déplacement sont utilisés dans l’évaluation de la norme \({\mathrm{\parallel }.\mathrm{\parallel }}_{\infty}\) , sauf dans deux cas où un traitement particulier est fait au niveau du choix des composantes:
Pour les chargements de type AFFE_CHAR_CINE , le degré de liberté concerné est ignoré dans l’évaluation de la norme du résidu car la procédure d’élimination des inconnues ne permet pas d’accéder aux réactions d’appuis;
Pour le contact continu, les composantes LAGR_C et LAGR_F1/LAGR_F2 sont ignorés dans l’évaluation de la norme car la loi de Signorini-Coulomb est déjà vérifié dans l’algorithme (voir [R5.03.52]) et que ces termes sont dimensionnellement incohérents avec ceux relatifs aux déplacements; Par contre, pour le cas du contact dans XFEM, ces composantes sont conservées car elles servent à vérifier la condition LBB;
Utilisation d’une matrice évolutive tangente-sécante#
La méthode décrite dans ce paragraphe s’applique exclusivement aux problèmes fortement non-linéaires, où une méthode de Newton classique échoue pour tout type de choix de matrice, pour la phase de prédiction ou de correction. Typiquement, la méthode de Newton habituelle est mise en défaut pour les problèmes mal-posés provenant de l’utilisation des lois de comportement adoucissantes (voir par exemple R5.04.01) .
Dans ces situations, une non-convergence se manifeste lorsque l’algorithme ne parvient pas à choisir entre plusieurs solutions admissibles, dans un incrément de (pseudo)-temps donné. Ce défaut de convergence au niveau global se traduit le plus souvent au niveau local (i.e. au point d’intégration) par une alternance répétée entre un état non-dissipatif (élastique) et un état dissipatif (plasticité, endommagement, …) au cours des itérations de Newton consécutives.
Une des stratégies consiste à utiliser la notion de pilotage (voir R5.03.80) pour pallier les insuffisances de Newton. L’autre stratégie consiste à modifier la matrice tangente. C’est cette dernière stratégie que l’on détaille ici.
En suivant l’état de chaque point d’intégration d’une itération à l’autre, on peut repérer les points responsables d’une non-convergence globale. Une fois ces points repérés, on décide de modifier la matrice en espérant que cette modification permettra une convergence globale.
Figure 2‑2.6-a: Loi de comportement adoucissante, utilisation d’une matrice mixte tangente-sécante
L’approche baptisée tangente-sécante, activée sous le mot-clef COMPORTEMENT avec TANGENTE_SECANTE=’OUI’, se justifie par le raisonnement suivant: si la méthode de Newton utilisant la matrice tangente cohérente ne converge pas, c’est souvent parce qu’un certain nombre (variable) de points d’intégration change d’état (non-dissipatif/dissipatif) entre deux itérations globales de Newton. Au niveau local (voir Figure 2-2.6-a ), cela veut dire qu’on continue à alterner entre le domaine 1 (\(\varepsilon <{\varepsilon}_{\text{pic}}\) ) et le domaine 2 (\(\varepsilon >{\varepsilon}_{\text{pic}}\) ). A cause du changement de pente entre 1 et 2, l’utilisation d’une matrice sécante ou tangente cohérente conduit à une approximation très pauvre, d’où l’intérêt de la modifier. Le choix que l’on présente consiste à construire la matrice tangente comme une combinaison linéaire entre la matrice tangente cohérente et la matrice sécante, les deux étant déterminées par les lois de comportement. Actuellement, l’approche n’est disponible que pour la loi de comportement adoucissante ENDO_ISOT_BETON.
Pour gérer l’amorçage de l’option tangente-sécante, on introduit une variable interne supplémentaire par rapport aux variables internes existantes, \(\alpha =({\alpha}_{1},\dots {\alpha}_{n},{\alpha}_{n+1})\) (avec \(n\) le nombre de variables internes du modèle utilisé). Cette variable traduit l’alternance éventuelle entre l’état élastique et l’état adoucissant d’un point de Gauss. On l’initialise à la première itération de Newton de chaque nouveau pas de temps, puis on la fait évoluer pour savoir, en chaque point de Gauss, le nombre d’alternances successives entre les deux états. En ayant cette information, on peut estimer qu’à partir d’un certain seuil (par exemple \({S}_{0}=3\) alternances), l’algorithme de Newton ne pourra plus converger pour l’incrément de temps courant et qu’il est nécessaire de modifier la matrice tangente. Pour modifier la matrice, on s’appuie directement sur la manière dont la matrice tangente cohérente dans ENDO_ISOT_BETON est construite (voir [R7.01.04]). Il s’agit de faire la somme entre la partie dissipative et la partie non dissipative.
où \({K}_{\text{tg}}\) est la matrice tangente, \({K}_{\text{sc}}\) la matrice sécante (partie non-dissipative) et \({K}_{\text{end}}\) la correction endommageante (partie dissipative). Pour la matrice modifiée \({K}_{\text{t-s}}\) , on remplace l’expression \({K}_{\text{tg}}\) dans () par:
où \(\eta\) est une fonction de \({\alpha}_{n+1}\) avec des valeurs entre 0 et 1. La fonction \(\eta\) retenue dans la suite est la suivante:
où \(A\) est une constante et \({S}_{0}\) le seuil sur la valeur du nombre d’alternances successives subies à partir duquel la matrice tangente est modifiée. Pour \({\alpha}_{n+1}={S}_{0}\) la matrice reste tangente cohérente et pour \({\alpha}_{n+1}\gg {S}_{0}\) , elle devient sécante. Des valeurs jugées satisfaisantes pour ces paramètres sont \(A=1,5\) et \({S}_{0}=3\) (valeurs par défaut). Le choix sur l’évolution de la valeur de \({\alpha}_{n+1}\) est primordial pour la performance de l’algorithme. On choisit une augmentation de \({\alpha}_{n+1}\) de \(G=1,0\) , \({\alpha}_{n+1}\to {\alpha}_{n+1}+G\) pour chaque nouvelle alternance entre un état élastique et un état endommageant, puis une diminution de \({\alpha}_{n+1}\) de \(P\) , \({\alpha}_{n+1}\to {\alpha}_{n+1}-P\) , lorsque l’état reste endommageant deux fois de suite. L’objectif de l’utilisation de \(P\) est de permettre le retour à la matrice tangente cohérente lorsqu’un point de Gauss reste endommageant sur plusieurs itérations, puisque la matrice tangente cohérente rend la convergence quadratique, à condition qu’on soit près de la solution. La valeur du taux de diminution \(P\) par rapport au taux d’augmentation \(G\) , est cruciale pour le comportement de l‘algorithme évolutif. L’idée globale consiste à augmenter \({\alpha}_{n+1}\) , lorsqu’on est loin de la solution pour avoir un opérateur plus proche du sécant que du tangent cohérent, puis une fois «proche» de la solution, de basculer en matrice tangente cohérente (qui est la meilleure au sens de Newton).Le rapport \(P/G\) (mot-clef TAUX_RETOUR – 0.05 par défaut) représente le troisième paramètre de l’algorithme, en plus de \(A\) (mot-clef AMPLITUDE) et \({S}_{0}\) (mot-clef SEUIL).
Méthode de Newton-Krylov#
Principe général#
La méthode de Newton-Krylov fait partie des méthodes de Newton inexactes. Elle est uniquement utilisable quand le solveur du système linéaire () est itératif (par opposition à un solveur direct). Cette approche ne modifie pas le choix de la matrice tangente du système comme les méthodes précédentes. Elle joue sur le critère de convergence du solveur itératif utilisé pour le système linéarisé. En adaptant au mieux le critère de convergence de la méthode itérative à la convergence de Newton, on évite de faire des itérations inutiles (dans le solveur linéaire) et on gagne ainsi en temps de calcul.
Mise en œuvre#
Le principe général des méthodes de Newton inexactes est de remplacer la condition que l’incrément de solution \(\delta {\mathrm{u}}_{i}^{n}\) soit la solution exacte du système () par une condition plus faible. On demande que \(\delta {\mathrm{u}}_{i}^{n}\) vérifie la condition:
où \({\eta}_{n}\) est appelé le forcing term .
On peut montrer que la méthode proposée est convergente et que quand la suite \({\eta}_{n}\) tend vers \(0\) , la convergence est super-linéaire (voir [ 8 ] p.96). Plus cette valeur est petite, plus la solution sera proche de celle obtenue par une résolution exacte, mais moins on gagnera de temps à la résolution du système linéaire. Il faut donc trouver un bon compromis entre résoudre rapidement les systèmes linéaires et ne pas trop dégrader la convergence des itérations de Newton.
En examinant la condition (), on constate qu’elle est identique au critère de convergence relatif des solveurs itératifs utilisés pour résoudre le système linéarisé de Newton. Pour vérifier cette condition, il s’agit donc d’utiliser le forcing term comme critère de convergence d’un solveur linéaire itératif.
Comme on l’a vu précédemment, il faut aussi s’assurer que la suite \({\eta}_{n}\) tende vers \(0\) pour conserver la convergence super-linéaire de la méthode de Newton. Pour ce faire, on va asservir \({\eta}_{n}\) à la décroissance du résidu de Newton par la loi d’évolution ( cf. [ 8 ] p.105):
où la constante est choisie telle que \(\gamma =0.1\) .
Cette simple formule n’est pas suffisante dans la pratique car il faut garantir une décroissance adéquate de \({\eta}_{n}\) . Pour cela, on détermine complètement \({\eta}_{n}\) par l’expression suivante:
où la constante vaut \({\eta}_{0}=0.9\) et correspond au critère de convergence utilisé pour la première résolution linéaire.
\({\eta}_{min}\) , quant à elle, est la valeur du critère de convergence du solveur linéaire itératif fournie par l’utilisateur (mot-clé RESI_RELA).
Conditions limites et chargements#
Classification#
Dans les paragraphes suivants, nous allons décrire l’ensemble des chargements possibles dans l’opérateur et leur classification. En premier lieu, il convient de distinguer trois grandes catégories de chargements:
Les chargements de Dirichlet (conditions limites) s’appliquent sur les inconnues nodales au sens large: ce peut être un déplacement, une température ou une pression (pour la THM). Ils peuvent être appliqués par élimination (AFFE_CHAR_CINE) ou par dualisation (AFFE_CHAR_MECA). Quand les conditions limites sont traitées par élimination, leur prise en compte se fait au moment de la résolution du système linéaire, de manière interne et ne modifie donc aucunement les équations déjà décrites précdemment;
Les chargements de Neumann s’appliquent en général aux flux. Ce sont des conditions en efforts, en général, définis sur les mailles (sauf pour le chargement FORCE_NODALE).
Il y a une autre distinction à réaliser entre les chargements suiveurs et ceux qui ne le sont pas. Un chargement est suiveur lorsqu’il dépend de la géométrie de la structure. Typiquement, la pression s’exerce dans la direction opposée à la normale. Donc, lorsque la structure se déforme avec l’évolution de la charge, le chargement, exprimé dans un repère absolu, est transformé. Il existe également des chargements de Dirichlet suiveurs. Ce sont des conditions limites qui dépendent de la géométrie parce qu’ils sont non-linéaires. Les charges qui ne dépendent pas de la géométrie de la structure sont appelées des charges mortes ou fixes (par exemple, la pesanteur).
Chargements suiveurs#
Pour indiquer qu’une charge doit être traitée comme une charge suiveuse dans STAT_NON_LINE, on indique TYPE_CHARGE=”SUIV” sous le mot-clé EXCIT. Un chargement mécanique \({L}^{\text{méca}}({u}_{i})\) comportant des charges suiveuses se décompose donc en deux parties:
L’exposant \({}^{\mathrm{fixe}}\) désigne ici les charges mortes, et \({}^{\text{suiv}}\) les charges suiveuses. Le système d’équations à résoudre devient alors:
Les opérations de dérivation permettant d’écrire la phase de prédiction et les itérations de la méthode de Newton font donc intervenir les dérivées de \({L}_{i}^{\text{méca}}({u}_{i})={L}_{i}^{\text{fixe}}+{L}_{i}^{\text{suiv}}({u}_{i})\) par rapport aux déplacements \(({u}_{i})\) .La phase de prédiction devient:
Et les itérations de Newton consistent à résoudre le système:
Ainsi, au début de chaque pas de charge (prédiction) et à chaque itération de Newton, on doit calculer une matrice de rigidité \({\frac{\partial {\mathrm{L}}^{\text{suiv}}}{\partial \mathrm{u}}\mid }_{\mathrm{u}}\) et un vecteur \({L}^{\text{suiv}}\) liés aux chargements suiveurs.
Les seules charges de Neumann qui peuvent être traitées comme des charges suiveuses dans l’état actuel de l’opérateur STAT_NON_LINE sont:
La pression pour les modélisations 3D, 3D_SI, D_PLAN, D_PLAN_SI, AXIS, AXIS_SI, C_PLAN, C_PLAN_SI [R3.03.04] et COQUE_3D [R3.03.07];
Le chargement de pesanteur pour les éléments CABLE_POULIE [R3.08.05], éléments à trois nœuds comportant une poulie et deux brins de câbles : la force de pesanteur s’exerçant sur l’élément dépend des longueurs respectives des deux brins;
La force centrifuge en grands déplacements, qui pour une vitesse de rotation \(\omega\) est donnée par: \({\int}_{\Omega}\rho .\omega \wedge \left[\omega \wedge \text{OM}\right].d\Omega ={\int}_{\Omega}\rho .\omega \wedge \left[\omega \wedge ({\text{OM}}_{0}+u)\right].d\Omega\) . Disponible pour les modélisations 3D et AXIS_FOURIER;
Le chargement de pesanteur pour toutes les modélisations THM des milieux poreux non saturés [R7.01.10] : en effet, la masse volumique dépend des variables nodales pour tenir compte des relations de comportement des géomatériaux.
Conditions limites dualisées – Cas linéaire#
Nous allons maintenant prendre en compte les conditions limites dualisées dans l’écriture du système linéaire à résoudre. Dans un premier temps, nous allons seulement considérer des relations linéaires . On suppose l’existence d’un potentiel \(J(v)\) tel que la minimisation de ce dernier permet d’exprimer l’équilibre mécanique (au sens faible) de la structure. Si \(J(u)\) est ce minimum alors on a:
\({L}^{int}(u,v)\) désigne la somme des forces internes discrétisées. Elle dépend de la solution recherchée via les lois de comportement, les cinématiques de grandes transformations, etc.
\({L}^{\text{ext}}(u,v)\) désigne la somme des forces externes discrétisées. Elle peut dépendre de la solution comme dans le cas des forces suiveuses, du contact, etc.
\(V\) désigne l’ensemble des solutions tests admissibles;
\(U\) désigne l’ensemble des solutions recherchées admissibles;
\(B\) est un opérateur linéaire de l’espace des champs de déplacements sur un espace de fonctions définies sur une partie du bord de la structure;
Remarque: en mécanique, on peut considérer la présence de forces dissipatives liées au caractère visqueux de la loi de comportement ou aux forces de frottement. Dans ces cas le potentiel est remplacé par l’existence d’un pseudo-potentiel de dissipation. On ne rentre pas dans les détails mais on montre qu’on peut se ramener dans le cadre d’un problème d’optimisation grâce à la transformation de Legendre-Fenchel. Par souci de simplicité, on garde le cadre général de l’équation précédente car la suite des développements n’est pas impactée par l’existence d’un pseudo-potentiel.
Sous forme discrète, la dualisation des conditions aux limites de Dirichlet \(\mathrm{B}.\mathrm{u}={\mathrm{u}}^{d}(t)\) conduit au problème suivant[R3.03.01] :
Les inconnues sont à présent, à tout instant \(t\) , le couple \((u,\lambda )\) , où \(\lambda\) représente les multiplicateurs de Lagrange des conditions aux limites de Dirichlet [R3.03.01]. Le vecteur \({B}^{T}.\lambda\) s’interprète comme l’opposé des réactions d’appui aux nœuds correspondants. En paramétrant par l’instant \({t}_{i}\) :
Trouver l’équilibre consiste donc à annuler en \(({u}_{i},{\lambda}_{i},{t}_{i})\) le vecteur résidu d’équilibre qui sera défini par :
Les inconnues sont calculées de façon incrémentale par l’algorithme de résolution global y compris les multiplicateurs de Lagrange des conditions limites dualisées. À partir de \(({u}_{i-1},{\lambda}_{i-1})\) , solution satisfaisant l’équilibre en \({t}_{i-1}\) , on détermine \(\Delta {u}_{i}\) et \(\Delta {\lambda}_{i}\) qui permettront d’obtenir la solution en \({t}_{i}\) :
Phase de prédiction d’Euler#
Pour la prédiction, il faut toujours linéariser le système () par rapport au temps mais, dans ce cas, autour des deux inconnues (déplacements et multiplicateurs de Lagrange) \(({u}_{i-1},{\lambda}_{i-1})\) . On note d’abord que la linéarisation des réactions d’appui \({B}^{T}.{\lambda}_{i}\) est immédiate car on considère des conditions limites linéaires et donc la matrice \(\mathrm{B}\) est constante (elle ne dépend pas des déplacements ou du temps). Comme \({\lambda}_{i}=\Delta {\lambda}_{i}^{0}+{\lambda}_{i-1}\) , il vient immédiatement:
On suppose que le chargement mécanique ne dépend pas du temps (les charges suiveuses sont exclues) et que les conditions limites de Dirichlet sont aussi linéaires, donc:
On obtient pour l’équation d’équilibre:
On a équilibre à l’instant \({t}_{i-1}\) , c’est-à-dire:
Et il reste donc:
Si on regarde maintenant la deuxième équation du système (), on obtient pour les conditions limites de Dirichlet:
On a équilibre à l’instant \({t}_{i-1}\) , c’est-à-dire:
Il reste finalement:
On obtient le système d’équations permettant de calculer des valeurs prédictives \((\Delta {\mathrm{u}}_{i}^{0},\Delta {\lambda}_{i}^{0})\) :
On remarquera alors que cette expression fait désormais intervenir les multiplicateurs de Lagrange à l’instant précédent , qui sont parfois inconnus (au premier incrément de charge, par exemple). Ce qui veut dire qu’avec cette nouvelle expression, on a déplacé le problème de la connaissance des forces internes à l’instant \({t}_{i-1}\) vers l’ignorance des multiplicateurs de Lagrange à ce même instant! Mais comme les conditions limites sont linéaires, le problème est simplifié. Considérons que la solution de () en ce qui concerne les multiplicateurs de Lagrange \({\stackrel{ˆ}{\lambda}}_{i}\) s’écrive sous forme incrémentale:
Cette solution résout la première équation du système:
L’idée est de rechercher ce \({\stackrel{ˆ}{\lambda}}_{i}\) . Comme l’opérateur \({B}^{T}\) est constant, en appliquant , on a:
On suppose que les conditions limites sont vérifiées, donc:
Les déplacements imposés s’écrivent aussi sous forme incrémentale:
La matrice \(B\) est constante, nous avions donc à l’incrément précédent (le problème a été résolu):
En utilisant et dans :
A l’équilibre, on a donc \({\stackrel{ˆ}{\lambda}}_{i}\) qui satisfait aussi les conditions limites que l’on ré-écrit:
Le vecteur des multiplicateurs de Lagrange \({\stackrel{ˆ}{\lambda}}_{i}\) peut donc être trouvé lors de la phase de prédiction en modifiant l’équation d’imposition des conditions limites par l’expression . Par analogie avec l’incrément des déplacements trouvés en prédiction \(\Delta {u}_{i}^{0}\) on notera \(\Delta {\lambda}_{i}^{0}={\stackrel{ˆ}{\lambda}}_{i}\) :
Un cas particulier concerne l’utilisation d’une excitation de type TYPE_CHARGE= “DIDI” signifiant Dirichlet différentiel, c’est-à-dire par rapport à l’état initial. Cela consiste, pour les conditions aux limites de type blocages, à imposer, non pas \(B.u={u}^{d}\) , mais \(B.(u-{u}_{\text{didi}})={u}^{d}\) Dans ce cas, le système à résoudre dans la phase de prédiction pour le nouvel incrément de charge est:
Phase de correction de Newton#
On doit linéariser le système par rapport aux inconnues en \(({u}_{i}^{n},{\lambda}_{i}^{n})\) à \({t}_{i}\) constant. On commence par linéariser les forces internes \({L}_{i}^{int,n}\) :
La linéarisation des réactions d’appui \({B}^{T}.{\lambda}_{i}^{n}\) est immédiate (conditions linéaires):
On suppose que le chargement mécanique ne dépend pas du temps (les charges suiveuses sont exclues) et que les conditions limites de Dirichlet sont aussi linéaires. On obtient pour l’équation d’équilibre:
Pour les conditions limites, la linéarisation du système () nous donne de manière analogue à :
Le système à résoudre s’écrit finalement:
Adaptation des critères de convergence#
On explicite ici l’adaptation de RESI_GLOB_MAXI et RESI_GLOB_RELA aux cas des conditions limites dualisées. Le critère RESI_GLOB_MAXI consiste à vérifier que la norme infinie du résidu est inférieure à la valeur \(\gamma\) spécifiée par l’utilisateur. Il est nécessaire d’introduire la contribution des conditions limites dualisées:
Le critère RESI_GLOB_RELA (choisi par défaut) revient à vérifier que le résidu est suffisamment petit, comme précédemment, et ceci relativement à une quantité représentative du chargement.
Conditions limites dualisées – Cas non-linéaire#
Considérons maintenant le cas de conditions limites dualisées non-linéaires. Par exemple, LIAISON_SOLIDE en grandes transformations. On suppose l’existence d’un potentiel \(J(v)\) tel que la minimisation de ce dernier permet d’exprimer l’équilibre mécanique (au sens faible) de la structure. Si \(J(u)\) est ce minimum alors on a:
\({L}^{int}(u,v)\) désigne la somme des forces internes discrétisées. Elle dépend de la solution recherchée via les lois de comportement, les cinématiques de grandes transformations, etc.
\({L}^{\text{ext}}(u,v)\) désigne la somme des forces externes discrétisées. Elle peut dépendre de la solution comme dans le cas des forces suiveuses, du contact, etc.
\(V\) désigne l’ensemble des solutions tests admissibles;
\(U\) désigne l’ensemble des solutions recherchées admissibles;
\(B\) est un opérateur linéaire de l’espace des champs de déplacements sur un espace de fonctions définies sur une partie du bord de la structure;
\(C(v)\) est un opérateur non-linéaire de l’espace des champs de déplacements sur un espace de fonctions définies sur une partie du bord de la structure;
Comme conditions limites non-linéaires, on considère l’exemple illustré sur la figure (), on veut que le solide reste rigide, c’est-à-dire que la distance entre le nœud \(M\) est lié aux autres nœuds \({M}_{i}\) reste constante. Soit:
Figure 3.4-1: Exemple de relations non-linéaires
Cette relation est non-linéaire car elle dépend des déplacements des différents nœuds \({M}_{i}\) .
Problème d’équilibre sous conditions d’égalités non linéaires#
On suppose l’existence d’une égalité non-linéaire discrète \(C(u)\) de \({ℝ}^{n}\) dans \(ℝ\) dépendant de la solution \(\textcolor[rgb]{0,0,0}{u\in {R}^{n}}\) tel que l’on ait le problème suivant:
On suppose que \(\textcolor[rgb]{0,0,0}{C(.)}\) est une fonction différentiable [bib3] de classe au moins un. Le problème précédent () et le problème () sont tels que l’on cherche le minimum dans une région de solutions admissibles \(\textcolor[rgb]{0,0,0}{V\subset {R}^{n}}\) . La difficulté réside dans la construction à priori de l’ensemble des solutions admissibles, ce qui n’est pas direct dans le cas du problème () à cause du caractère inconnu (dépendant de la solution) de cet espace admissible. On va utiliser la technique de dualisation par les multiplicateurs de Lagrange, on écrit le Lagrangien augmenté de ce système:
D’où le problème suivant équivalent au problème ():
Pour cela, il suffit d’établir les conditions de stationnarité du lagrangien, soit:
D’où le système non-linéaire suivant à résoudre:
Qu’on écrira sous forme plus compacte:
Linéarisation du problème#
Pour résoudre ce problème on utilise également la méthode de Newton. Mais il faut linéariser le système (). On se place au pas de temps \({t}_{i}\) (par simplification des écritures, on omet de préciser la dépendance des quantités au pas de temps). On approxime les quantités à l’itération de Newton \(n\) , par rapport à ces mêmes quantités à l’itération \(n-1\) :
On pose:
En développant tous les termes:
On a finalement le système suivant:
Avec:
\({K}^{int,n-1}={\partial}_{u}{L}^{int,n-1}\) la matrice tangente classique (voir 2.3.2 et 2.4.3 );
\(\nabla {\mathrm{C}}^{n-1}={\partial}_{u}{\mathrm{C}}^{n-1}\) la matrice jacobienne de la relation non-linéaire;
\(\nabla \nabla {\mathrm{C}}^{n-1}{\lambda}_{\mathit{nl}}^{n-1}={\partial}_{{u}^{2}}^{2}{\mathrm{C}}^{n-1}{\lambda}_{\mathit{nl}}^{n-1}\) la matrice Hessienne de la relation non linéaire
On retrouve le problème approché suivant, trouver les incréments \(\Delta u\in {ℝ}^{n}\) , \(\Delta {\lambda}_{l}\in ℝ\) et \(\Delta {\lambda}_{\mathit{nl}}\in ℝ\) du système en prédiction (linéarisation autour du pas de temps \(i\) ). Au § 2.3.5 , nous avons trouvé l’expression de l’incrément de chargement mécanique \(\Delta {L}_{i}^{\text{méca}}\) :
Pour les conditions limites linéaires, le paragraphe § 3.3.1 en a aussi donné l’expression.
Il convient de noter la différence entre \(\mathrm{B}\) , matrice des relations linéaires et \({\mathrm{C}}_{i}\) matrice des relations non-linéaires. En effet, la matrice \(\mathrm{B}\) est constante pendant tout le calcul, tandis que la matrice \({\mathrm{C}}_{i}\) change à chaque itération de Newton. C’est pour ça que cette dernière est indicée par \(i\) , contrairement à la matrice \(\mathrm{B}\) . De manière équivalente, pendant les itérations de Newton \(n\) , on doit résoudre le système suivant:
Remarques:
Comme dans le cas standard, il est possible de faire du quasi-Newton en ignorant la contribution de \(\nabla \nabla \mathrm{C}\) à la matrice tangente. Mais c’est un choix du développeur car l’utilisateur n’a pas accès au choix de réactualisation de ce terme contrairement au cas général;
Le problème devenant non-linéaire en termes de conditions limites, il est nécessaire de basculer avec une matrice réactualisée toutes les itérations de Newton (REAC_ITER=1). Un message d’alarme vous prévient si vous ne le faites pas;
Adaptation du critère de convergence#
Il est important de bien définir le critère d’arrêt à la fin de chaque itération de Newton aussi bien après la prédiction que pendant la correction. On rappelle (voir § 2.5 ) que dans Code_Aster, on peut le définir de plusieurs manières. On explicite ici l’adaptation de RESI_GLOB_MAXI et RESI_GLOB_RELA aux cas non-linéaires.
Le critère RESI_GLOB_MAXI consiste à vérifier que la norme infinie du résidu est inférieure à la valeur \(\gamma\) spécifiée par l’utilisateur. Il est nécessaire d’introduire la contribution des conditions limites non-linéaires:
Le critère RESI_GLOB_RELA (choisi par défaut) revient à vérifier que le résidu est suffisamment petit, comme précédemment, et ceci relativement à une quantité représentative du chargement.
Recherche linéaire#
La recherche linéaire exposé ici concerne la recherche linéaire en l’absence de pilotage.
Principe#
L’introduction de la recherche linéaire dans l’opérateur STAT_NON_LINE résulte d’un constat : la méthode de Newton avec matrice consistante ne converge pas dans tous les cas de figure, notamment lorsque l’on part trop loin de la solution. D’autre part, l’utilisation de matrices autres que la matrice tangente consistante peut, lorsqu’elles sont trop «souples», conduire à divergence. La recherche linéaire permet de se prémunir contre de telles divergences.
Elle consiste à considérer \((\delta {u}_{i}^{n},\delta {\lambda}_{i}^{n})\) , non plus comme l’incrément des déplacements et des multiplicateurs de Lagrange, mais comme une direction de recherche dans laquelle on va chercher à minimiser une fonctionnelle (l’énergie de la structure). On trouvera un pas d’avancement \(\rho\) dans cette direction, et l’actualisation des inconnues consistera à faire:
En l’absence de recherche linéaire (par défaut) le scalaire \(\rho\) est bien sûr égal à 1.
Minimisation d’une fonctionnelle#
Afin de mieux se convaincre du bien fondé de la recherche linéaire, on peut interpréter la méthode de Newton comme une méthode de minimisation d’une fonctionnelle (dans le cas où les matrices tangentes sont symétriques). Nous insistons sur le fait que les équations obtenues sont rigoureusement celles de la méthode de Newton exposée dans le [§ 2 ] et que seule la façon d’y parvenir est différente.
«Oublions» pour simplifier l’exposé la dualisation des conditions aux limites de Dirichlet et plaçons-nous dans l’hypothèse des petites déformations. On considère la fonctionnelle :
où la densité d’énergie libre \(\Phi\) permet de relier le tenseur des contraintes \(\sigma\) au tenseur des déformations linéarisées \(\varepsilon\) par la relation \(\varepsilon =\frac{\partial \Phi }{\partial \sigma }\) dans le cas de l’hyperélasticité (on généralise cette situation aux autres non linéarités dans la suite du document). La fonctionnelle \(J\) étant convexe, trouver le minimum de \(J\) est équivalent à annuler son gradient, soit:
Ce qui est exactement le Principe des Travaux Virtuels puisque :
Ainsi, résoudre les équations issues du Principe des Travaux Virtuels (base du problème formulé dans le [§ 1.3 ]) est équivalent à minimiser la fonctionnelle \(J\) qui représente l’énergie de la structure (énergie interne diminuée du travail des forces extérieures \(\mathrm{f}\) et \(t\) ).
Méthode de minimisation#
La minimisation se fait de façon itérative, classiquement en deux temps à chaque itération:
Calcul d’une direction de recherche \(\delta\) le long de laquelle on va chercher l’itéré suivant,
Calcul du meilleur pas d’avancement \(\rho\) dans cette direction : \({u}^{n+1}={u}^{n}+\rho .\delta\)
Dans un problème de minimisation, l’idée naturelle est d’avancer dans la direction opposée au gradient de la fonctionnelle, qui est localement la meilleure direction de descente puisque cette direction est normale aux lignes d’isovaleurs et dirigée dans le sens des valeurs décroissantes Figure 2-4.3-1
Figure 2 ‑ 4.3-1
Cependant, il est possible d’améliorer le choix de la direction de descente en utilisant cette méthode de gradient dans une nouvelle métrique. C’est ce qui va nous permettre de retrouver les équations classiques de la méthode de Newton. Prenons l’exemple simple d’un problème à deux variables \(x\) et \(y\) pour lequel la fonctionnelle a la forme d’une ellipse dont le minimum est en \((\frac{\alpha}{a},\frac{\beta}{b})\) :
En choisissant comme direction de descente l’inverse du gradient de \(J\) , on passe d’un itéré au suivant (raisonnons sur \(x\) seulement) par :
qui ne pointe pas vers la solution puisque la normale en un point d’une ellipse ne passe en général pas par le centre de l’ellipse
Figure 2‑ 4.3-2
Par contre, si l’on effectue un changement de variables pour que les isovaleurs de \(J\) deviennent des cercles :
Figure 2 ‑ 4.3-3
L’utilisation de la direction inverse du gradient de \(\stackrel{ˉ}{J}\) permet alors d’obtenir la solution en une itération:
Ainsi, l’utilisation de la méthode de gradient dans la nouvelle métrique permet une convergence immédiate. Dans un cas plus compliqué (fonctionnelle convexe mais différente d’une ellipse), la convergence n’est pas instantanée mais le changement de variables permet de réduire sensiblement le nombre d’itérations nécessaires.
Application à la minimisation de l’énergie#
Pour simplifier, on va se placer dans le cas linéaire discrétisé où la fonctionnelle \(J\) vaut :
On note \(K\) la matrice de rigidité de la structure, et \(L\) le vecteur des chargements imposés. Pour minimiser \(J\) , nous allons utiliser la même méthode de descente que précédemment en faisant au préalable un changement de variables tout à fait similaire. La matrice \(K\) étant symétrique définie positive, ses valeurs propres sont réelles positives : on peut donc définir la « racine carrée » de \(K\) que l’on notera \(\sqrt{K}\) (également symétrique). On pose \(\stackrel{ˉ}{u}=\sqrt{K}.u\) , la minimisation de \(J\) est alors équivalente à celle de :
En utilisant une décomposition par la diagonale:
Avec \(\mathrm{D}\) diagonale, donc:
Ce qui conduit à:
En prenant comme direction de descente la direction inverse du gradient de \(\stackrel{ˉ}{J}\) , on obtient:
Soit, en revenant aux variables initiales:
Ou encore :
On retrouve les équations de la méthode de Newton : la matrice \(K\) est la Hessienne de la fonctionnelle \(J\) (matrice de la dérivée seconde) et le second membre est la différence du chargement et des forces internes, appelé aussi résidu d’équilibre. Ainsi la méthode de Newton peut-être interprétée comme résultant de la minimisation de l’énergie de la structure via une méthode de gradient appliquée après un changement de métrique.
Détermination du pas d’avancement#
Revenons au problème réel, celui de la résolution de \({L}_{i}^{int}({u}_{i})={L}_{i}^{\text{ext}}\) . Ce problème peut être interprété comme la minimisation de la fonctionnelle suivante:
où \(W({u}_{i})\) correspond à la discrétisation, sur la base des fonctions de forme, de l’énergie interne de la structure:
On a calculé par la méthode de Newton un incrément de déplacement \(\delta {u}_{i}^{n}\) qui, dans le problème de minimisation, s’interprète comme une direction de recherche, d’après ce qui précède. On va calculer le pas d’avancement \(\rho\) dans cette direction permettant de minimiser la valeur de la fonctionnelle:
Pour trouver le minimum de cette fonction scalaire de \(\rho\) que l’on notera \(f(\rho )\) , on cherche le point où sa dérivée s’annule (cela revient à rendre orthogonaux le résidu final et la direction de recherche):
\(f'(\rho )\) est la projection du résidu sur la direction de recherche. Avec les notations du [§ 2.2 ] et en prenant en compte les réactions d’appui, l’équation scalaire à résoudre pour déterminer le pas d’avancement \(\rho\) , s’écrit :
À la fin de la recherche linéaire, on actualise les déplacements et paramètres de Lagrange par:
Le test de convergence porte:
Sur le nombre maximum d’itérations de recherche linéaire indiqué par l’utilisateur sous le mot‑clé ITER_LINE_MAXI du mot-clé facteur NEWTON (la valeur par défaut 0 inhibe la recherche linéaire, et \(\rho\) vaut alors 1),
Sur le critère RESI_LINE_RELA donné par \(f(\rho )\le \tau .f(0)\) , où \(\tau\) vaut par défaut 0.1.
La recherche linéaire est en quelque sorte une « assurance » permettant de se prémunir contre des divergences graves de la méthode de Newton. Lorsque la direction de recherche est « mauvaise » (si la matrice tangente est trop souple, par exemple), l’algorithme de recherche linéaire aboutit à une faible valeur de \(\rho\) , ce qui évite de trop s’éloigner de la solution. Il n’est pas nécessaire de faire beaucoup d’itérations dans la méthode de sécante (deux ou trois suffisent pour éviter les catastrophes) car chacune coûte assez cher (il faut recalculer les forces internes) et on n’a pas l’ambition de trouver à chaque itération de Newton la valeur de \(\rho\) vraiment optimale.
Calcul du coefficient de recherche linéaire#
Il existe deux méthodes pour calculer le \(\rho\) optimal dans la recherche linéaire .
Méthode sécante (METHODE=’CORDE’)#
Afin que la détermination de \(\rho\) ne soit pas trop coûteuse, on utilise une méthode de sécante dont le nombre maximum d’itérations est fixé par l’utilisateur. La méthode de sécante peut s’interpréter comme une méthode de Newton où la dérivée au point courant est approchée par la direction joignant le point courant et le point précédent:
Où l’on a noté \({g}^{p}=f'({\rho}^{p})\) . On part de \({\rho}^{0}=0\) et \({\rho}^{1}=1\) . La méthode de sécante a un ordre de convergence de l’ordre de 1.6. Elle se représente schématiquement sur la Figure 2-4.6.1-1 .
Figure 2 ‑ 4.6.1-1
Méthode mixte (METHODE=’MIXTE’)#
Cette méthode mélange plusieurs techniques de résolution pour être plus robuste. Elle consiste essentiellement en l’application d’une méthode de la sécante (voir paragraphe précédente) entre deux bornes pré-déterminées. Il s’agit en fait de l’application de la méthode de sécante avec des bornes variables. Voici l’algorithme utilisé:
On suppose que \(f'(0)>0\) . Si ce n’est pas le cas, on change le sens de la direction de descente (en examinant des \(\rho\) négatifs, ce qui revient à définir \(f'\) comme étant égale à \(-f'\)
On cherche un \({\rho}_{\max}\) positif tel que \(f'({\rho}_{\max})<0\) . La méthode est simplement itérative en faisant \({\rho}_{n+1}=3.{\rho}_{n}\) avec \({\rho}_{0}=1\) (étape de bracketage ou encadrement)
On a ainsi les deux nouvelles bornes entre lesquelles la fonction change de signe. Si on suppose que la fonction \(f'\) est continue, il existe donc une solution entre ces bornes.
On applique la méthode de la sécante sur cet intervalle : on part de \({\rho}^{0}=0\) et \({\rho}^{1}={\rho}_{\max}\) .
Cas particulier : la méthode de Newton-Krylov#
Il a été précisé plus haut que la recherche linéaire est réalisée simultanément sur les inconnues \(\mathrm{u}\) et \(\lambda\) comme le montre la formule () d’actualisation des variables. Or la fonctionnelle à minimiser ne présente pas de minimum en fonction des inconnues \((\mathrm{u},\lambda )\) , il s’agit en effet d’un Lagrangien qui présente un point selle, c’est-à-dire un minimum en \(\mathrm{u}\) et un maximum en \(\lambda\) (voir [R3.03.01]). Cette manière de faire n’est donc pas licite dans le cas général.
Cela dit, on peut montrer que dans le cas où le système en prédiction est résolu «exactement» (tout du moins à une précision numériquement satisfaisante), cette approche est licite. C’est généralement le cas dans l’usage habituel de Code_Aster .
Ce n’est par contre pas le cas dans le cadre de l’usage de la méthode de Newton-Krylov, où les systèmes linéaires sont justement résolus de manière volontairement inexacte. Dans cette situation, pour contourner le problème, seules les inconnues \(\mathrm{u}\) sont affectées par la recherche linéaire et la formule de mise à jour des variables devient donc:
Dans la mesure où la fonctionnelle à minimiser dispose bien d’un minimum en \(\mathrm{u}\) , la procédure de recherche linéaire est licite.
Dans la mesure où la fonctionnelle à minimiser dispose bien d’un minimum en \(\mathrm{u}\) , la procédure de recherche linéaire est licite.
Pilotage#
On se reportera à la documentation [R5.03.80].
Régularisation viscoélastique REGU_VISC#
En présence d’instabilités matérielles, résultant par exemple de l’endommagement, la solution ne dépend plus continûment du temps (ou des paramètres de chargement) mais exhibe des sauts. Le pilotage est une technique pour franchir ces sauts, en autorisant des évolutions non monotones du chargement pour retrouver la continuité. Mais ces techniques ne sont pas applicables lorsque le temps physique intervient, que ce soit à travers un historique de variables de commande (transitoire thermique, par exemple) ou dans le comportement du matériau (viscosité).
La régularisation visqueuse est une autre approche pour retrouver la continuité. Elle consiste à superposer aux contraintes issues de la loi de comportement \(\sigma\) une contrainte d’origine visqueuse \({\sigma}_{\mathit{vis}}\) de sorte à pénaliser suffisamment les évolutions très rapides (dont la limite est un saut de solution):
Il s’agit là d’une technique numérique, sans rapport a priori avec une réalité visqueuse du matériau. On décrit dans le fascicule [R5.03.34] la loi de comportement viscoélastique qui gouverne la contrainte visqueuse \({\sigma}_{\mathit{vis}}\) .
Bibliographie#
Ibrahimbegovic – «Nonlinear Solid Mechanics: Theoretical Formulations and Finite Element Solution Methods» – Springer – 2009
M.A. Crisfield – «Non-Linear Finite Element Analysis of Solids and Structures : Essentials» – Wiley Professional Software – 1996
M.A. Crisfield – «Non-Linear Finite Element Analysis of Solids and Structures : Advanced Topics» – Wiley Professional Software – 1997
J.C. Simo, T.J.R. Hughes – «Computational inelasticity» – Springer – 2000
«Éléments d’analyse et de résolution numérique des relations de l’élasto-plasticité» – EDF – Bulletin de la Direction des Études et Recherches – Série C – N° 3 1986 – p. 57 - 89.
J.F. Maître– «Analyse numérique» – cours polycopié de l’ENTPE.
Shi, M. A. Crisfield – “Combining arc-length control and line searches in path following” – Communications in Numerical Methods in Engineering, Vol. 11, 793-803 (1995).
C.T. Kelley - «Iterative Methods for Solving Linear and Nonlinear Equations», Vol. 16 in Frontiers in Applied Mathematics,Siam, Philadelphia, 1995.
Historique des versions du document#
Version Aster |
Auteur(s) ou contributeur(s), organisme |
Description des modifications |
5 |
|
|
7.4 |
P.Badel, J.Laverne,N. Tardieu EDF R&D AMA |
|
8.5 |
M.Abbas EDF R&D AMA |
Mise à jour des notations, ajout des variables de commande. |
9.2 |
M.Abbas EDF R&D AMA |
Mise à jour du § 2.2. |
9.4 |
M.Abbas, D.Markovic EDF R&D AMA |
Ajout de la matrice tangente-sécante et de la recherche linéaire mixte. |
10.1 |
M.Abbas EDF R&D AMA |
Mise à jour du paragraphe sur les grandes déformations. |
11.1 |
N. Tardieu EDF R&D AMA |
Ajout de la méthode de Newton-Krylov. |
11.3 |
M. David EDF R&D AMA |
Précisions concernant le critère RESI_REFE_RELA |