r5.03.80 Méthodes de pilotage du chargement

Résumé:

Ce document décrit les méthodes de pilotage du chargement disponibles dans Code_Aster (par un degré de liberté, par longueur d’arc, par incrément de déformation et par prédiction élastique). Elles introduisent une inconnue supplémentaire, l’intensité de la part pilotable du chargement, et une équation supplémentaire, la contrainte de pilotage. Ces méthodes permettent en particulier de calculer la réponse d’une structure qui présenterait des instabilités, aussi bien d’origines géométrique (flambement) que matériau (adoucissement).

Principe des méthodes de pilotage du chargement

De manière générale, les fonctionnalités de pilotage disponibles dans Code_Aster permettent de déterminer l’intensité d’une partie du chargement pour satisfaire une contrainte portant sur les déplacements. Leur emploi est limité à des simulations pour lesquelles le temps ne joue pas de rôle physique, ce qui exclut a priori les problèmes dynamiques ou visqueux. Il est également incompatible avec les problèmes exprimant une condition d’unilatéralité (contact et frottement). On peut distinguer plusieurs gammes d’utilisation auxquelles répondent autant de méthodes de pilotage (mot clé facteur PILOTAGE):

• Contrôle physique des forces par le déplacement d’un point de la structure: pilotage par degré de liberté imposé (TYPE=”DDL_IMPO”);

• Suivi d’instabilités géométriques (flambement), la réponse de la structure pouvant exhiber des snap-back «doux»: pilotage par longueur d’arc (TYPE=”LONG_ARC”);

• Suivi d’instabilités matériaux (en présence de lois de comportement adoucissantes), la réponse de la structure pouvant exhiber des snap-back «brutaux»: pilotage par la prédiction élastique (TYPE= “PRED_ELAS”) ou plus généralement par l’incrément de déformation, (TYPE=”DEFORMATION”);

• Calcul des charges limites des structures (TYPE=”ANA_LIM”), cf. [ref:R7.07.01 <R7.07.01>].

Plus précisément, les méthodes de pilotage reposent sur les deux idées suivantes.

D’une part, on considère que le chargement (forces extérieures et déplacements donnés) se décompose additivement en deux termes, l’un connu et imposé par l’utilisateur (\({L}_{\text{impo}}^{\text{méca}}\) et \({u}_{\text{impo}}^{d}\) ) et l’autre (\({L}_{\text{pilo}}^{\text{méca}}\) et \({u}_{\text{pilo}}^{d}\) ) dont seule la direction est connue, son intensité \(\eta\) devenant une nouvelle inconnue du problème:

(2847)\[\begin{split}\lbrace \begin{array}{c}{L}^{\text{méca}}={L}_{\text{impo}}^{\text{méca}}+\eta .{L}_{\text{pilo}}^{\text{méca}}\\ {u}^{d}={u}_{\text{impo}}^{d}+\eta .{u}_{\text{pilo}}^{d}\end{array}\end{split}\]

D’autre part, afin de pouvoir résoudre le problème, on lui associe une nouvelle équation qui porte sur les déplacements et qui dépend de l’incrément de temps: c’est la contrainte de pilotage, qui s’exprime par:

(2848)\[P(\Delta U)=\Delta \tau \text{avec}P(0)=0\]

où \(\Delta \tau\) est indirectement une donnée utilisateur qui s’exprime via le pas de temps courant \(\Delta t\) et un coefficient de pilotage (COEF_MULT) tel que:

(2849)\[\Delta \tau =\frac{\Delta t}{\text{COEF\_MULT}}\]

La condition \(P(0)=0\) est requise afin d’obtenir un incrément de déplacement d’autant plus petit que le pas de temps est petit. Finalement, les inconnues du problème deviennent les déplacements \(u\) , les multiplicateurs de Lagrange \(\lambda\) associés aux conditions aux limites et l’intensité du chargement piloté \(\eta\) , baptisée ETA_PILOTAGE. Le système non linéaire à résoudre s’écrit dorénavant:

(2850)\[\begin{split}\lbrace \begin{array}{cc}{L}^{int}(u)+{B}^{T}.\lambda & ={L}_{\text{impo}}^{\text{méca}}+\eta .{L}_{\text{pilo}}^{\text{méca}}\\ B.u& ={u}_{\text{impo}}^{d}+\eta .{u}_{\text{pilo}}^{d}\\ P(\Delta u)& =\Delta \tau \end{array}\end{split}\]

Remarques:

• A l’heure actuelle, les chargements suiveurs (i.e. qui dépendent des déplacements) et les conditions de Dirichlet de type «* DIDI » ne sont pas pilotables.

• Le chargement ne dépend plus directement du temps mais résulte de la résolution de tout le système non linéaire. Cela implique que la part pilotée du chargement ne doit pas dépendre du temps physique, mais correspond à un effort que l’on ajuste pour satisfaire une contrainte cinématique supplémentaire.*

Résolution du système global

L’introduction d’une nouvelle équation ne perturbe pas outre mesure la méthode de résolution du système non linéaire. En effet, on procède comme en [R5.03.01], c’est-à-dire que la résolution est incrémentale. On note le pas de temps \(i\) en indice et l’itération de Newton \(n\) en exposant. Le problème non-linéaire est alors résolu en deux temps:

• Une phase de prédiction qui donne une première estimation des déplacements et des multiplicateurs de Lagrange notée \((\Delta {u}_{i}^{0},\Delta {\lambda}_{i}^{0})\) ;

• Une phase de correction de Newton \((\delta {u}_{i}^{n},\delta {\lambda}_{i}^{n})\) qui vient corriger cette première estimation;

Avec un nombre \({n}_{\mathrm{CV}}\) suffisant d’itérations de Newton, on obtient un résultat convergé [1]_ :

(2851)\[\begin{split}\lbrace \begin{array}{}{u}_{i}^{\text{convergé}}={u}_{i-1}+\Delta {u}_{i}^{0}+\sum_{j=1}^{n={n}_{\mathrm{CV}}}\delta {u}_{i}^{j}\\ {\lambda}_{i}^{\text{convergé}}={\lambda}_{i-1}+\Delta {\lambda}_{i}^{0}+\sum_{j=1}^{n={n}_{\mathrm{CV}}}\delta {\lambda}_{i}^{j}\end{array}\end{split}\]

Le principe est donc de linéariser le système () que l’on écrit au pas de temps \(i\) :

(2852)\[\begin{split}\lbrace \begin{array}{cc}{L}_{i}^{int}+{B}^{T}.{\lambda}_{i}& ={L}_{\text{impo},i}^{\text{méca}}+{\eta}_{i}.{L}_{\text{pilo},i}^{\text{méca}}\\ B.{u}_{i}& ={u}_{\text{impo},i}^{d}+{\eta}_{i}.{u}_{\text{pilo},i}^{d}\\ P(\Delta {u}_{i})& =\Delta {\tau}_{i}\end{array}\end{split}\]

Il n’y aura pas de linéarisation par rapport à la variable de pilotage \({\eta}_{i}\) . De la sorte, on préserve toute la méthodologie de réactualisation de l’opérateur tangent déjà mise en œuvre pour les calculs sans pilotage. De plus, la structure «bande» de la matrice tangente est conservée. Le schéma de résolution n’est donc plus strictement une méthode de Newton (ce qui n’est pas gênant car nous verrons que l’équation de pilotage est, au plus, de degré deux). En prédiction, si on linéarise le système () sans l’équation de pilotage, par rapport au temps autour de \(({u}_{i-1},{\lambda}_{i-1})\) , on obtient (voir tout le développement dans [R5.03.01]):

(2853)\[\begin{split}\lbrace \begin{array}{cc}{\mathrm{K}}_{i-1}.\Delta {\mathrm{u}}_{i}^{0}+{\mathrm{B}}^{T}.\Delta {\lambda}_{i}^{0}& ={\mathrm{L}}_{\text{impo},i}^{\text{méca}}+{\eta}_{i}.{\mathrm{L}}_{\text{pilo},i}^{\text{méca}}-{\mathrm{Q}}_{i-1}^{T}.{\sigma}_{i-1}+\Delta {\mathrm{L}}_{i}^{\text{varc}}\\ \mathrm{B}.\Delta {\mathrm{u}}_{i}^{0}& ={\mathrm{u}}_{\text{impo},i}^{d}+{\eta}_{i}.{\mathrm{u}}_{\text{pilo},i}^{d}-\mathrm{B}.{\mathrm{u}}_{i-1}\end{array}\end{split}\]

En correction, on linéarise le système () toujours sans l’équation de pilotage, par rapport au temps, mais autour de \(({u}_{i}^{n},{\lambda}_{i}^{n})\) , on a :

(2854)\[\begin{split}\lbrace \begin{array}{cc}{\mathrm{K}}_{i}^{n-1}.\delta {\mathrm{u}}_{i}^{n}+{\mathrm{B}}^{T}.\delta {\lambda}_{i}^{n}& ={\mathrm{L}}_{\text{impo},i}^{\text{méca}}+{\eta}_{i}.{\mathrm{L}}_{\text{pilo},i}^{\text{méca}}-{\mathrm{L}}_{i}^{int,n-1}-{\mathrm{B}}^{T}.{\lambda}_{i}^{n-1}\\ \mathrm{B}.\delta {\mathrm{u}}_{i}^{n}& ={\eta}_{i}.{\mathrm{u}}_{\text{pilo},i}^{d}\end{array}\end{split}\]

On va réunir les deux systèmes en une écriture commune, afin de simplifier l’exposé. Le système à résoudre s’écrit finalement:

(2855)\[\begin{split}\left[\begin{array}{cc}{K}_{i}^{n-1}& {B}^{T}\\ B& 0\end{array}\right].\left\lbrace \begin{array}{}{\delta u}_{i}^{n}\\ {\delta \lambda }_{i}^{n}\end{array}\right\rbrace =\left\lbrace \begin{array}{}{L}_{\text{impo},i}^{n-1}\\ {u}_{\text{impo},i}\end{array}\right\rbrace +{\eta}_{i}.\left\lbrace \begin{array}{}{L}_{\text{pilo},i}\\ {u}_{\text{pilo},i}\end{array}\right\rbrace\end{split}\]

On peut faire plusieurs remarques:

• La matrice \({K}_{i}^{n-1}\) dépend à la fois du pas de temps courant et éventuellement de l’itération de Newton précédente. Les différentes manières de la construire (quasi-Newton, élastique, sécante, cohérente, tangente en vitesse, etc.) sont décrites dans [R5.03.01];

• On a supposé que les chargements extérieurs étaient linéaires (ils ne dépendent que du pas de temps). On ne considère donc pas les chargements suiveurs tels que la pression ou la force centrifuge bien que du point de vue théorique, ça ne pose pas de difficultés. Par contre, le matériau peut être décrit avec un comportement non-linéaire, ce qui implique que \({\mathrm{L}}_{i}^{int,n-1}\) dépend de l’itération de Newton (résultat de la linéarisation des efforts intérieurs).

• Les conditions limites de Dirichlet sont toujours linéaires, ce qui permet à la matrice \(B\) d’être constante sur tout le transitoire.

• Avec le bon choix de la matrice tangente et du second membre, formellement, il y a équivalence entre \((\delta {u}_{i}^{n=0},\delta {\lambda}_{i}^{n=0})\) et l’incrément en prédiction \((\Delta {u}_{i}^{0},\Delta {\lambda}_{i}^{0})\) .

On peut maintenant exprimer les corrections de déplacements \(\delta {u}_{i}^{n}\) et de multiplicateurs de Lagrange \(\delta {\lambda}_{i}^{n}\) en fonction de \({\eta}_{i}\) moyennant la résolution du système linéaire () par rapport à chacun des deux seconds membres. C’est-à-dire que l’on sépare les deux solutions:

(2856)\[\begin{split}\left\lbrace \begin{array}{}{\delta u}_{i}^{n}\\ {\delta \lambda }_{i}^{n}\end{array}\right\rbrace =\left\lbrace \begin{array}{}{\delta u}_{\text{impo},i}^{n}\\ {\delta \lambda }_{\text{impo},i}^{n}\end{array}\right\rbrace +{\eta}_{i}.\left\lbrace \begin{array}{}{\delta u}_{\text{pilo},i}^{n}\\ {\delta \lambda }_{\text{pilo},i}^{n}\end{array}\right\rbrace\end{split}\]

Ces deux solutions correspondent au découplage des deux chargements:

(2857)\[\begin{split}\left\lbrace \begin{array}{}{L}_{i}^{n-1}\\ {u}_{i}\end{array}\right\rbrace =\left\lbrace \begin{array}{}{L}_{\text{impo},i}^{n-1}\\ {u}_{\text{impo},i}\end{array}\right\rbrace +{\eta}_{i}.\left\lbrace \begin{array}{}{L}_{\text{pilo},i}\\ {u}_{\text{pilo},i}\end{array}\right\rbrace\end{split}\]

On résout indépendamment [2] :

(2858)\[\begin{split}\left[\begin{array}{cc}{K}_{i}^{n-1}& {B}^{T}\\ B& 0\end{array}\right].\left\lbrace \begin{array}{}{\delta u}_{\text{impo},i}^{n}\\ {\delta \lambda }_{\text{impo},i}^{n}\end{array}\right\rbrace =\left\lbrace \begin{array}{}{L}_{\text{cst},i}^{n-1}\\ {u}_{\text{cst},i}\end{array}\right\rbrace\end{split}\]

Et:

(2859)\[\begin{split}\left[\begin{array}{cc}{K}_{i}^{n-1}& {B}^{T}\\ B& 0\end{array}\right].\left\lbrace \begin{array}{}{\delta u}_{\text{pilo},i}^{n}\\ {\delta \lambda }_{\text{pilo},i}^{n}\end{array}\right\rbrace =\left\lbrace \begin{array}{}{L}_{\text{pilo},i}^{n-1}\\ {u}_{\text{pilo},i}\end{array}\right\rbrace\end{split}\]

On peut maintenant substituer la correction de déplacement \(\Delta {u}_{i}^{n}\) dans l’équation de contrôle du pilotage du système; il en résulte une équation scalaire en \({\eta}_{i}\) :

(2860)\[\tilde{P}({\eta}_{i})\underset{\text{déf.}}{=}P(\Delta {\mathrm{u}}_{i}^{n-1}+\delta {\mathrm{u}}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {\mathrm{u}}_{\text{pilo},i}^{n})=\Delta {\tau}_{i}\]

La méthode de résolution de cette équation dépend de la nature du contrôle de pilotage adopté cf. [§ 7 ]. Finalement, il ne reste plus qu’à réactualiser les inconnues déplacements et multiplicateurs de Lagrange:

(2861)\[\begin{split}\lbrace \begin{array}{}\Delta {u}_{i}^{n}=\Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {u}_{\text{pilo},i}^{n}\\ {\lambda}_{i}^{n}={\lambda}_{i}^{n-1}+\delta {\lambda}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {\lambda}_{\text{pilo},i}^{n}\end{array}\end{split}\]

[1]

La notion de résultat « convergé » est plus amplement détaillée dans [R5.03.01].

[2]

En pratique, la factorisation de la matrice n’est faite qu’une fois et l’on résout simultanément les deux systèmes.

Équation de contrôle du pilotage

Pilotage par contrôle d’un degré de liberté des déplacements : DDL_IMPO

Pour ce premier type de pilotage, la fonction \(P\) se borne à extraire un degré de liberté de l’incrément de déplacement. En particulier, il s’agit donc d’une fonction linéaire:

(2862)\[P(\Delta {u}_{i}^{n})=\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n}\right\rbrace =\Delta {\tau}_{i}\]

Où le vecteur nodal \(\langle S\rangle\) est un vecteur de sélection qui est nul partout excepté pour le degré de liberté à extraire où il vaut un:

(2863)\[\langle S\rangle =\langle 0\dots 0\dots \underset{\text{noeud n, ddl i}}{1}\dots 0\dots 0\rangle\]

L’équation () se réduit alors également à une équation linéaire:

(2864)\[{\eta}_{i}=\frac{\Delta {\tau}_{i}-\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n-1}\right\rbrace -\langle S\rangle .\left\lbrace \delta {u}_{\text{impo},i}^{n}\right\rbrace }{\langle S\rangle .\left\lbrace \delta {u}_{\text{pilo},i}^{n}\right\rbrace }\]

On notera qu’il n’y a pas de solution lorsque la correction de déplacement piloté \(\delta {u}_{\text{pilo},i}^{n}\) ne permet pas d’ajuster le degré de liberté requis, ce qui peut arriver si, par erreur, on bloque le degré de liberté en question. Dans ce cas, le code s’arrêtera en ECHEC DU PILOTAGE.

Pilotage par longueur d’arc : LONG_ARC

Une autre forme de pilotage très largement utilisée consiste à contrôler la norme de l’incrément de déplacement (par rapport à certains nœuds et certaines composantes): on parle alors de pilotage par longueur d’arc cylindrique, voir Bonnet and Wood [bib1]. Plus précisément, la fonction \(P\) s’exprime par:

(2865)\[P(\Delta {u}_{i}^{n})={\parallel \Delta {u}_{i}^{n}\parallel }_{\text{<S>}}=\sqrt{(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n}\right\rbrace ).(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n}\right\rbrace )}=\Delta {\tau}_{i}\]

Où, à nouveau, le vecteur nodal \(\langle S\rangle\) permet de sélectionner les degrés de liberté employés pour le calcul de la norme (il vaut un pour les degrés de liberté sélectionnés et zéro ailleurs). Dans ce cas là on remarquera que \(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n}\right\rbrace\) n’est pas un scalaire, mais un vecteur correspondant au produit des composantes de \(\langle S\rangle\) par celles de \(\left\lbrace \Delta {u}_{i}^{n}\right\rbrace\) . Dans ce cas, l’équation de pilotage se réduit à une équation du second degré:

(2866)\[\begin{split}\begin{array}{}{\eta}_{i}^{2}.\left[(\langle S\rangle .\left\lbrace \delta {u}_{\text{pilo},i}^{n}\right\rbrace ).(\langle S\rangle .\left\lbrace \delta {u}_{\text{pilo},i}^{n}\right\rbrace )\right]+\\ 2.{\eta}_{i}.\left[(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}\right\rbrace ).(\langle S\rangle .\left\lbrace \delta {u}_{\text{pilo},i}^{n}\right\rbrace )\right]+\\ \left[(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}\right\rbrace ).(\langle S\rangle .\left\lbrace \Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}\right\rbrace )-\Delta {\tau}_{i}^{2}\right]=0\end{array}\end{split}\]

Cette équation peut ne pas admettre de solution. Dans ce cas, on choisit la valeur \({\eta}_{i}\) qui minimise le polynôme (). On vérifie alors bien \(\tilde{P}({\eta}_{i})>\Delta {\tau}_{i}\) . Dans le cas contraire, elle admet deux racines (ou une racine double). On choisit celle des deux qui minimise:

• soit la norme de l’incrément de déplacement \(∥\Delta {\mathrm{u}}_{i}^{n}∥\) ;

• soit le résidu d’équilibre \(R({u}_{i},{t}_{i})={L}_{i}^{int}-{L}_{i}^{\text{méca}}\) ;

• soit l’angle formé par \(\Delta {u}_{i-1}\) et \(\Delta {u}_{i}^{n}\) (où \(\Delta {u}_{i-1}\) est l’incrément de déplacement solution du pas de temps précédent), c’est-à-dire celle qui maximise le cosinus de cet angle dont l’expression est:

(2867)\[\cos(\Delta {u}_{i-1},\Delta {u}_{i}^{n})=\frac{\langle \Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {u}_{\text{pilo},i}^{n}\rangle .\left\lbrace \Delta {u}_{i-1}\right\rbrace }{∥\Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n}+{\eta}_{i}.\delta {u}_{\text{pilo},i}^{n}∥.∥\Delta {u}_{i-1}∥}\]

Pour calculer l’angle, on n’utilise que les inconnues correspondant aux déplacements, à l’exclusion de toutes les autres (en particulier les rotations dans les éléments de structures).

On remarque que ce dernier critère angulaire sélectionne la mauvaise solution lorsqu’on passe d’une phase de charge à une phase de décharge ou vice-versa, c’est à dire lorsque le coefficient de pilotage \({C}_{i}\) change de signe. La valeur \({C}_{i-1}\) est donc mise en mémoire et actualisée à chaque pas de temps convergé, son initialisation pour chaque appel à STAT_NON_LINE étant transparente pour l’utilisateur : \({C}_{i-1}\) est ainsi initialisé lors de la lecture de l’état initial, \({C}_{i}\) ayant été inclus dans les paramètres de calcul à archiver dans le concept résultat. Si \({C}_{i}{C}_{i-1}>0\) , on maximise le critère () . Dans le cas contraire, on le minimise afin de sélectionner l’incrément de déplacement en sens opposé de l’incrément précédent.

Pilotage par l’incrément de déformation : DEFORMATION

Le pilotage par incrément de déformation consiste à exiger que l’incrément de déformation du pas courant reste proche en direction de la déformation au début du pas de temps, et ce pour au moins un point de Gauss de la structure. On requiert ainsi qualitativement qu’ a minima un point de la structure conserve le mode de déformation qu’il avait au préalable (par exemple, traction dans une direction donnée).

Cas des petites déformations

Pour les petites déformations, on peut rendre compte de cette exigence moyennant le choix de la fonction de pilotage suivante:

(2868)\[P(\Delta {u}_{i}^{n})=\underset{g}{\text{Max}}\frac{{\varepsilon}_{g,i-1}}{\parallel {\varepsilon}_{g,i-1}\parallel }.\Delta {\varepsilon}_{g,i}=\Delta {\tau}_{i}\]

Où l’indice \(g\) balaie les points de Gauss de la structure (ou seulement les mailles précisées par le mot-clé GROUP_MA dans PILOTAGE) et où la déformation en un point de Gauss se déduit du vecteur nodal des déplacements via l’emploi des matrices «partie symétrique du gradient des fonctions de forme» \({B}_{g}\) (à ne pas confondre avec la matrice des conditions de Dirichlet):

(2869)\[{\varepsilon}_{g,i-1}={B}_{g}\cdot {u}_{i-1}\]

Le contrôle du pilotage en fonction de \({\eta}_{i}\) s’écrit alors:

(2870)\[P(\Delta {\mathrm{u}}_{i}^{n})=\underset{g}{\text{Max}}({A}_{g}^{(0)}+{\eta}_{i}\cdot {A}_{g}^{(1)})=\underset{g}{\text{Max}}({L}_{g}({\eta}_{i}))=\Delta {\tau}_{i}\]

Avec les deux termes:

(2871)\[{A}_{g}^{(0)}=\frac{{\varepsilon}_{g,i-1}}{\parallel {\varepsilon}_{g,i-1}\parallel }\cdot {B}_{g}(\Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n})\]

Et:

(2872)\[{A}_{g}^{(1)}=\frac{{\varepsilon}_{g,i-1}}{\parallel {\varepsilon}_{g,i-1}\parallel }\cdot {B}_{g}(\delta {u}_{\text{pilo},i}^{n})\]

Ce pilotage ne dépend pas de la loi de comportement, pourvu qu’elle utilise le tenseur des petites déformations (option DEFORMATION=”PETIT” ou DEFORMATION=”PETIT_REAC”).

Cas des grandes déformations

En présence de grandes déformations, on peut généraliser la fonction de pilotage () en employant des déformations de Green-Lagrange (mesure lagrangienne des déformations dans la configuration initiale):

(2873)\[P(\Delta {\mathrm{u}}_{i}^{n})=\underset{g}{\text{Max}}\frac{{\mathrm{E}}_{g,i-1}}{\mathrm{\parallel }{\mathrm{E}}_{g,i-1}\mathrm{\parallel }}.\Delta {\mathrm{E}}_{g,i}\]

La mesure des déformations \(E\) est écrite en fonction du tenseur gradient de transformation \(F\) :

(2874)\[\mathrm{E}=\frac{1}{2}\cdot ({\mathrm{F}}^{T}.\mathrm{F}-\mathrm{I})\]

Toutefois, on n’aboutirait plus comme précédemment à une fonction affine par morceaux. Pour y remédier, on décide de procéder à une linéarisation de \(\Delta {E}_{g,i}\) par rapport à \(\Delta {u}_{i}^{n-1}\) . \(P\) a alors une expression similaire à () avec:

(2875)\[{A}_{g}^{(0)}=\frac{{E}_{g,i-1}}{\parallel {E}_{g,i-1}\parallel }\cdot \text{sym}\left[{F}_{g,i}^{T}.{\nabla}_{g}(\Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n})\right]\]

Et:

(2876)\[{A}_{g}^{(1)}=\frac{{E}_{g,i-1}}{\parallel {E}_{g,i-1}\parallel }\cdot \text{sym}\left[{F}_{g,i}^{T}.\delta {u}_{\text{pilo},i}^{n}\right]\]

Où \({\nabla}_{g}(u)\) désigne le gradient (non symétrisé) des déplacements \(u\) évalué au point de Gauss d’indice \(g\) . Tout comme dans le cas précédent, ce pilotage ne dépend pas de la loi de comportement, pourvu qu’elle utilise le tenseur des grandes déformations.

Résolution de l’équation non-linéaire de pilotage

La fonction \(P(\Delta {u}_{i}^{n})\) est convexe et linéaire par morceaux. Elle admet généralement aucune, une ou deux solutions, cf. figure (). Lorsqu’elle n’admet pas de solutions, on arrête l’algorithme de STAT_NON_LINE (échec du pilotage): si l’utilisateur utilise la subdivision automatique du pas de temps, ce dernier sera alors activé.

Dans le cas où elle admet deux solutions, trois possibilités sont offertes: on choisit celle des deux qui minimise:

• soit la norme de l’incrément de déplacement \(∥\Delta {\mathrm{u}}_{i}^{n}∥\) ;

• soit le résidu d’équilibre \(R({u}_{i},{t}_{i})={L}_{i}^{int}-{L}_{i}^{\text{méca}}\) ;

• soit l’angle formé par \(\Delta {u}_{i-1}\) et \(\Delta {u}_{i}^{n}\) ;

Pour résoudre l’équation (), on propose l’algorithme présenté ci-dessous. Il est basé sur la construction d’intervalles emboîtés: les bornes du dernier d’entre eux sont les solutions de l’équation et, comme annoncé précédemment, on choisit celle qui conduit au \(\tilde{u}(\eta )\) le plus proche de \({u}_{i-1}\) . Cet algorithme, rapide, s’appuie sur la résolution de \(G\) équations scalaires linéaires, où \(G\) désigne le nombre total de points de Gauss. L’algorithme peut prendre fin prématurément lorsque l’un des intervalles est vide, ce qui signifie que l’équation () n’admet pas de solutions.

	Initialisation de l’intervalle	\({I}_{0}=]-\infty ,+\infty [\)
	Boucle sur les points de Gauss \(g\)
(2.0)	La pente est nulle	\(\text{Si}{A}_{g}^{(1)}=0\)
(2.0.1)	Si \({A}_{g}^{(0)}>\Delta \tau\) échec
	Sinon on continue
(2.1)	Racine de la fonction linéaire active	\({\eta}_{g}\text{tel que}{L}_{g}({\eta}_{g})=\Delta \tau\)
(2.2)	Construction de l’intervalle suivant
(2.2.1)	Si la fonction linéaire active est croissante	\(\text{Si}{A}_{g}^{(1)}>0\mathrm{\Rightarrow }{I}_{g}={I}_{\text{g-1}}\cap ]-\infty ,{\eta}_{g}]\)
(2.2.2)	Si la fonction linéaire active est décroissante	\(\text{Si}{A}_{g}^{(1)}<0\Rightarrow {I}_{g}={I}_{\text{g-1}}\cap [{\eta}_{g},+\infty [\)
(2.3)	Arrêter si l’intervalle est vide
	Les solutions sont les bornes de l’intervalle[3]_	\(\eta \in \text{Fr}({I}_{G})\mathrm{\Rightarrow }\underset{g}{\max}{L}_{g}(\eta )=\Delta \tau\)

Algorithme de résolution de l’équation affine par morceaux

Pilotage par la prédiction élastique : PRED_ELAS

Algorithme

Si le pilotage par l’incrément de déformation s’avère suffisant pour suivre des solutions dissipatives dans la plupart des instabilités matériaux, l’existence de solutions n’est néanmoins pas prouvée. On lui préfère alors une méthode de pilotage fondée sur la prédiction élastique pour laquelle l’existence de solutions est démontrée mais qui, en revanche, est spécifique à chaque loi de comportement (implantée uniquement pour les lois ENDO_SCALAIRE, ENDO_FRAGILE, ENDO_ISOT_BETON, ENDO_ORTH_BETON, BETON_DOUBLE_DP, VMIS_ISOT_* et celles relatives aux éléments de joints [R5.06.09]). Plus précisément, lorsque la loi de comportement est gouvernée par un seuil, on définit \(P\) comme le maximum sur tous les points de Gauss de la valeur de la fonction seuil dans le cas d’un essai élastique (réponse incrémentale élastique du matériau).

Ainsi, considérons que l’état du matériau est décrit par la déformation \(\varepsilon\) et un ensemble de variables internes \(a\) . Appelons respectivement \(\sigma (\varepsilon ,a)\) et \(A(\varepsilon ,a)\) les contraintes et les forces thermodynamiques associées à \(a\) . Supposons en outre que les lois d’évolution de \(a\) soient gouvernées par un seuil \(f(A,\varepsilon ,a)\) :

(2877)\[f(A,\varepsilon ,a)\le 0\text{et}\lambda \cdot f(A,\varepsilon ,a)=0\]

Et une fonction d’écoulement \(G(A,\varepsilon ,a)\) :

(2878)\[\dot{a}=\lambda G\cdot (A,\varepsilon ,a)\]

Une telle formulation englobe la plupart des modèles de comportement dissipatifs et indépendants de la vitesse de chargement. La fonction seuil vaut alors pour un essai élastique:

(2879)\[{f}^{\text{el}}(\varepsilon )=f(A(\varepsilon ,{a}_{i-1}),\varepsilon ,{a}_{i-1})\]

On simplifie le problème en linéarisant \({f}^{\text{el}}\) par rapport à \(\varepsilon\) au voisinage du point \({\varepsilon}^{\text{*}}\) tel que \({f}^{\text{el}}=\Delta \tau\) :

(2880)\[{f}_{L}^{\text{el}}(\varepsilon )\underset{\text{def}}{=}{f}^{\text{el}}({\varepsilon}^{\text{*}})+{(\frac{\partial f}{\partial A}\cdot \frac{\partial A}{\partial \varepsilon }+\frac{\partial f}{\partial \varepsilon })\mid }_{{\varepsilon}^{\text{*}}}\cdot (\varepsilon -{\varepsilon}^{\text{*}})\]

Le choix de \({\varepsilon}^{\text{*}}\) permet de résoudre le problème exact, comme montré dans [bib3]. Il nécessite donc de résoudre un problème non-linéaire local pour chaque point de Gauss, puis l’algorithme global se fonde sur des droites, ce qui lui assure une grande efficacité. Il existe zéro, une ou deux solutions au problème local \({f}^{\text{el}}=\Delta \tau\) . Dans le cas d’absence de solution au problème local, le problème global ne peut avoir de solution: on arrête. Sinon (cas à une ou deux solutions), on linéarise autour de chacune des solutions.

Finalement, la fonction de contrôle du pilotage est définie comme le maximum de \({f}_{L}^{\text{el}}\) par rapport à tous les points de Gauss \(g\) , fonction qui ne dépend que de \(\Deltau\) :

(2881)\[P(\Delta {u}_{i}^{n})=\underset{g}{\text{Max}}{f}_{L}^{\text{el}}({\varepsilon}_{g,i-1}+{B}_{g}\cdot \Delta {u}_{i}^{n})=\Delta {\tau}_{i}\]

Finalement, l’équation de contrôle du pilotage s’écrit:

(2882)\[P(\Delta {u}_{i}^{n})=\underset{g}{\text{Max}}({A}_{g}^{(0)}+{\eta}_{i}\cdot {A}_{g}^{(1)})=\underset{g}{\text{Max}}({L}_{g}({\eta}_{i}))=\Delta {\tau}_{i}\]

Avec:

(2883)\[{A}^{(0)}={f}^{\text{el}}({\varepsilon}^{\text{*}})+{(\frac{\partial f}{\partial A}\cdot \frac{\partial A}{\partial \varepsilon }+\frac{\partial f}{\partial \varepsilon })\mid }_{{\varepsilon}^{\text{*}}}\cdot ({B}_{g}({u}_{i-1}+\Delta {u}_{i}^{n-1}+\delta {u}_{\text{impo},i}^{n})-{\varepsilon}^{\text{*}})\]

Et:

(2884)\[{A}^{(1)}={(\frac{\partial f}{\partial A}\cdot \frac{\partial A}{\partial \varepsilon }+\frac{\partial f}{\partial \varepsilon })\mid }_{{\varepsilon}^{\text{*}}}\cdot ({B}_{g}\cdot \delta {u}_{\text{pilo},i}^{n})\]

Expressions en fonction des lois de comportement

On est ainsi ramené à un problème identique à celui du pilotage par l’incrément de la déformation. On emploie bien entendu les mêmes algorithmes de résolution que ceux présentés dans le paragraphe précédent (voir § 11 ).

Les caractéristiques de l’équation de pilotage dépendent de la loi de comportement car l’expression de la fonction seuil élastique \({f}^{\text{el}}\) est elle-même une donnée du problème non-linéaire de comportement:

• Pour les éléments à discontinuité interne PLAN_ELDI et AXIS_ELDI, la loi de comportement est de type CZM_EXP et l’expression de la loi de pilotage se trouvera dans la documentation [R7.02.14]. La fonction seuil élastique s’écrit \({f}_{\mathit{el}}({\sigma}_{n},{\sigma}_{t})=\sqrt{\text{<}{\sigma}_{n}{\text{>}}_{+}^{2}+{\sigma}_{t}^{2}}-R(\kappa )\) ;

• Pour les éléments de joint *_JOINT, les lois de comportement sont des lois CZM régularisées CZM_EXP_REG et CZM_LIN_REG (voir [R7.02.11]), et pour les éléments de type interface *_INTERFACE, les lois de comportement sont des lois CZM non-régularisées comme CZM_OUV_MIX et CZM_TAC_MIX. L’expression de la loi de pilotage se trouvera dans la documentation [R7.02.11]. La fonction seuil élastique s’écrit \({f}_{\mathit{el}}\mathrm{\simeq }\underset{\mathit{pt}\mathit{gauss}}{\max}\left\lbrace \frac{\mathrm{\parallel }{\delta}^{0}+\eta {\delta}^{1}\mathrm{\parallel }-{\kappa}^{-}}{{G}_{c}/{\sigma}_{c}+{\kappa}^{-}}\right\rbrace =\frac{\Delta \tau }{C}\) ;

• Pour la loi ENDO_ISOT_BETON ([R7.01.04]), il y a une particularité. Au lieu de chercher un paramètre de pilotage \(\eta\) qui fasse sortir le critère d’une valeur \(\mathrm{\Delta \tau }\) avec l’endommagement issu du pas de temps précédent, on cherche un paramètre \(\eta\) qui nous ramène sur le critère avec un endommagement augmenté de \(\mathrm{\Delta \tau }\) . La fonction seuil élastique s’écrit donc \({\tilde{f}}^{\text{el}}(\eta ,{d}^{-})=\mathrm{\Delta \tau }\Rightarrow {\tilde{f}}^{\text{el}}(\eta ,{d}^{-}+\mathrm{\Delta \tau })=0\) ;

• Pour la loi ENDO_ORTH_BETON, voir [R7.01.09];

• Pour les autres lois d’endommagement de type élastique comme ENDO_SCALAIRE ou ENDO_FRAGILE(avec formulation locale ou à gradients), tout comme pour la loi ENDO_ISOT_BETON, la fonction seuil élastique est choisie de manière à chercher un paramètre de pilotage contrôlé par la valeur d’endommagement et non pas par la sortie du critère. La fonction seuil élastique s’écrit:

(2885)\[\begin{split}\begin{array}{c}{\tilde{f}}^{\text{el}}(\eta )=\frac{1}{2}({e}_{0}+\eta {e}_{1})\cdot E\cdot ({e}_{0}+\eta {e}_{1})-s\\ \lbrace \begin{array}{cccc}\text{loi locale}& {e}_{0}={\epsilon}_{0}& {e}_{1}={\epsilon}_{1}& s=k({d}^{-})\\ \text{gradient d'endommagement}& {e}_{0}={\epsilon}_{0}& {e}_{1}={\epsilon}_{1}& s=k({d}^{-})-c\Delta {d}^{-}\\ \text{déformation régularisée}& {e}_{0}={\stackrel{ˉ}{\epsilon}}_{0}& {e}_{1}={\stackrel{ˉ}{\epsilon}}_{1}& s=k({d}^{-})\end{array}\end{array}\end{split}\]

Illustration de l’algorithme de maximisation

Afin d’illustrer notre propos et notamment de comprendre pourquoi un tel choix de \({\varepsilon}^{\text{*}}\) permet de résoudre le problème exact, on prend l’exemple d’une loi de comportement cohésive CZM_EXP_REG ou CZM_LIN_REG (cf [R7.02.11]) implémentée sur un élément de joint (cf [R3.06.09]). Le saut de déplacement \(\delta\) est la variable primale (notée génériquement \(\epsilon\) dans ce qui précède), et \(\kappa\) la variable interne (notée génériquement \(A\) ). La fonction seuil s’écrit alors :

(2886)\[{f}^{\mathrm{èl}}(\delta )=\frac{{\parallel \delta \parallel }_{\mathrm{èq}}-{\kappa}_{i-1}}{{\kappa}_{\mathrm{REF}}}\text{avec}{\kappa}_{i-1}=\underset{v\in [0,{t}_{i-1}]}{\max}{\parallel \delta (v)\parallel }_{\mathrm{èq}};{\parallel \delta \parallel }_{\mathrm{èq}}=\sqrt{{\langle {\delta}_{n}\rangle }_{+}^{2}+{\delta}_{\tau}^{2}}\]

Sans détailler la démonstration, on peut montrer que \({\parallel \delta \parallel }_{\mathrm{èq}}\) vérifie une inégalité triangulaire. Ceci nous permet de prouver la convexité de la fonction :

(2887)\[{f}^{\mathrm{èl}}(\eta )={f}^{\mathrm{èl}}(\delta (\eta ))=\frac{{\parallel {\delta}_{i-1}+\Delta {\delta}_{i}^{n}+d{\delta}^{\mathrm{cte}}+\eta d{\delta}^{\mathrm{pilo}}\parallel }_{\mathrm{èq}}-{\kappa}_{i-1}}{{\kappa}_{\mathrm{REF}}}\]

Nous avons tracé en figure une représentation graphique qualitative des fonctions \({f}^{\mathrm{èl}}\) pour les différents points de Gauss des éléments de joint, ainsi que les solutions du problème exact de maximisation global sur l’ensemble des points de Gauss :

\(\underset{g}{\max}{f}_{g}^{\mathrm{èl}}(\eta )\stackrel{\mathrm{dèf}}{=}\underset{g}{\max}{f}^{\mathrm{èl}}({\delta}_{g}(\eta ))=\Delta {\tau}_{i}\)

Il apparaît alors que les solutions \(\eta\) de ce problème exact sont identiques aux solutions obtenues avec l’algorithme de résolution de l’équation affine par morceaux, c’est-à-dire les solutions de l’équation () \(\underset{g}{\max}{A}_{1}^{(g)}\eta +{A}_{0}^{(g)}=\Delta {\tau}_{i}\) , avec dans notre exemple :

\({A}_{1}^{(g)}=\frac{\partial {f}_{g}^{\mathrm{èl}}}{\partial \eta }({\eta}_{g}^{\text{*}});{A}_{0}^{(g)}=\Delta \tau -{A}_{1}^{(g)}{\eta}_{g}^{\text{*}}\) où \({\eta}_{g}^{\text{*}}\) est une solution de l’équation locale \({f}_{g}^{\mathrm{èl}}(\eta )=\Delta {\tau}_{i}\)

Nous renvoyons une fois de plus le lecteur à [3] pour une démonstration rigoureuse.

Figure 4.6.3-a : Résolution de l’équation de pilotage, représentation graphique qualitative

Choix de solution

Dans le cas où le problème global possède deux solutions, les algorithmes de choix décrits en sont réutilisés. En général, la sélection par résidu est retenue en raison de sa robustesse. Cependant, le fait de choisir la solution qui minimise le résidu ne nous garantit pas la convergence de la méthode de Newton. Certains cas ont pu être exhibés pour lesquels la solution de résidu initial maximal converge tandis que ce n’est pas le cas de la solution de résidu initial minimal.

En règle générale pour les calculs pilotés le pas de temps ne remplit plus le rôle de chargement, mais juste un paramètre de suivi de la solution. Bien que les solutions mécaniques soient identiques, on peut ainsi obtenir une différence significative au même pas de pilotage pour les problèmes avec les instabilités importantes, en changeant par exemple la machine d’exécution. Ce qui est illustré sur la figure , ou lors du snap-back dans la réponse globale force-déplacement on observe le décalage d’instant de calcul entre les deux courbes.

Figure 4.6.4-a : Illustration de décalage de l’instant de pilotage sur une courbe force-déplacement pour un problème fortement instable en fonction de la précision machine

[1]

On note \(\eta \in \text{Fr}({I}_{G})`le fait que :math:\)eta`appartient aux bornes de l’intervalle \({I}_{G}\)

Bibliographie

1. Bonnet J. and Wood R. D. [1997]. Nonlinear continuum mechanics for finite element analysis. Cambridge university press.

2. Bonnans J. F., Gilbert J.-C., Lemaréchal C. et Sagastizabal C. [1997]. Optimisation numérique: aspects théoriques et pratiques. Mathématiques et Applications, 27, ed. Springer.

3. Lorentz E. and Badel P. [2001]. A new path-following constraint for strain-softening finite element simulations. International journal for numerical methods in engineering. 2004, vol. 60, no 2, pp. 499 – 526.

4. Shi J. and Crisfield M. A. [1995]. Combining arc-length control and line searches in path following. Comm. Num. Meth. Eng. 11, pp. 793-803.