u2.08.04 Notice de calcul de charges critiques et de modes de flambement d’Euler et de calcul de flambage#

Résumé :

L’objectif de cette documentation est de présenter un guide méthodologique pour le calcul de charges critiques et de modes de flambement d’Euler d’une part et l’analyse de flambage non linéaire (instabilités) d’une structure d’autre part. On y aborde principalement les deux fonctionnalités de Code_Aster :

  • l’analyse de flambement linéaire, dite d’Euler, au travers de CALC_MODES (avec le mot-clé TYPE_RESU=‘MODE_FLAMB’),

  • le calcul de l’évolution quasistatique (opérateur STAT_NON_LINE) de la structure qui présente des non-linéarités géométriques et comportementales, dont on cherche les instabilités (option CRIT_STAB), un point limite, voire la réponse post-critique.

La première étape est, généralement, un calcul de flambement d’Euler, qui permettra de connaître les modes de flambement et les charges critiques correspondantes. Du point de vue du concepteur, la connaissance du premier mode et de sa charge critique est souvent suffisante, afin de se définir une marge de fonctionnement par rapport au chargement imposé : le coefficient multiplicateur entre le chargement imposé et la charge critique la plus faible donne la marge de sécurité.

Remarques

  • La connaissance du premier mode de flambement d’Euler peut aussi servir d’indication pour optimiser la gestion du calcul incrémental non linéaire mené par la suite. En effet, à l’approche de la charge critique, on peut alors décider de modifier le pilotage ou de réduire le pas de temps, voire d’augmenter le nombre d’itérations de vérification de l’équilibre dans la méthode de résidu, à chaque pas de charge.

  • L’allure du mode de flambement d’Euler peut aussi servir pour imposer un défaut géométrique initial sur la structure, afin de s’assurer, entre autre, que le calcul non linéaire incrémental de flambage bifurquera bien sur ce mode.

L’analyse d’Euler étant par définition linéaire, elle ne permet pas de prendre en compte des relations de comportement inélastiques, du contact unilatéral ou l’aspect suiveur des forces exercées. Il est alors nécessaire de faire un calcul incrémental non linéaire, qui en quasi-statique s’appuiera sur la commande STAT_NON_LINE de Code_Aster . C’est la méthode classique incrémentale par résidu en équilibre.

Les points particuliers de son utilisation seront abordés par la suite, en particulier l’utilisation de l’analyse de stabilité non linéaire avec le mot-clé CRIT_STAB (qui est aussi disponible dans DYNA_NON_LINE pour les études dynamiques).

On trouvera un tableau des cas-tests disponibles sur ce type de simulations au § u2.08.04-_tude_non_lin_aire_quasi-statique_de_la_structure de ce document.

Étude non linéaire quasistatique de la structure#

Cette étape se justifie si la structure présente de fortes non-linéarités, dont l’analyse de flambement d’Euler ne peut tenir compte. L’opérateur de résolution des problèmes non linéaires en quasistatique se nomme STAT_NON_LINE [bib7].

Ces non-linéarités peuvent être liées au matériau qui peut avoir un comportement élastoplastique [bib8], comme dans l’exemple qui va suivre. La prise en compte du contact, voire du frottement, est une autre source de non-linéarités. On peut aussi citer le cas des chargements suiveurs, comme la pression ([bib1] et [bib2] pour les éléments de type coques volumiques), qui nécessitent une approche non linéaire dans Code_Aster.

On peut distinguer deux types d’analyses de stabilité non linéaires, pouvant se combiner.

D’une part la généralisation de l’analyse de flambement d’Euler présentée précédemment : on va faire du suivi du trajet de chargement pour atteindre des points limites ou des charges critiques avec les modes associés au cours du calcul non linéaire incrémental. Cela se traduit par une analyse de stabilité sur les matrices de raideur réactualisées. Ce type d’analyse se fait souvent sur une structure sans défaut initial.

D’autre part, on peut tenir compte de défauts introduits sur le modèle parfait (en insérant des imperfections géométriques initiales), afin de « forcer » la bifurcation de la branche de solutions et de faire du suivi de la branche stable pour analyser la réponse post-critique. Le pilotage du trajet de chargement quasistatique doit être traité de manière adaptée (opérateur STAT_NON_LINE, cf. [U4.51.03]).

Bien évidemment, ce suivi de solution post-critique peut être initié par l’analyse aux valeurs propres des matrices de raideur réactualisées, en particulier pour bien détecter la bifurcation et définir le défaut alors introduit en se basant sur le mode de flambement observé.

Analyse de stabilité sur matrices de raideur réactualisées#

Que ce soit en quasistatique (opérateur STAT_NON_LINE) ou en dynamique (opérateur DYNA_NON_LINE), Code_Aster permet de mener des analyses incrémentales de stabilité au sens du flambage sur les matrices de raideur courantes. Ces étapes de calcul sont gérées par le mot-clé facteur CRIT_STAB avec l’option TYPE = “FLAMBEMENT” ( cf. [U4.51.03] et [U4.53.01]).

Remarques pour l’analyse de stabilité Dans Code_Aster, on se restreint aujourd” hui au cas où la matrice de raideur réactualisée est symétrique, car on appelle le solveur modal de type SORENSEN, cf. [U4.52.02]. On traite donc des bifurcations des branches d’équilibre, c’est-à-dire une perte d’unicité. Cela restreint donc les familles de modèles de comportement que l’on peut traiter (règle d’écoulement dite « associée »), ainsi qu’aux cas de chargements suiveurs tels qu’il existe un potentiel : par exemple pour une enceinte fermée sous pression. Pour un chargement suiveur enlevant le caractère conservatif du système donc conduisant à une matrice de raideur non symétrique, l’instabilité devient dynamique.

Remarques pour l’analyse dynamique transitoire non linéaire

L’analyse de stabilité des trajectoires repose sur l’étude du système linéarisé des équations du mouvement et du problème aux valeurs propres associé. L’instabilité dynamique peut se produire si on exhibe une valeur propre complexe : * avec une partie réelle strictement positive et sans partie imaginaire : le système est dit instable par divergence, la trajectoire possédant alors une évolution exponentiellement croissante non-oscillante ; * complexe avec une partie réelle strictement positive : le système est dit instable par flottement, la trajectoire possédant alors une évolution exponentiellement croissante et oscillante. C’est par exemple le cas du crissement produit par un solide glissant et frottant sur un autre ; * pour des lois de comportement irréversible à règle d’écoulement non associée la perte d’hyperbolicité peut advenir à l’apparition d’une paire de valeurs propres complexes conjuguées à partie imaginaire négative liées à la non-symétrie du tenseur tangent : il s’agit d’une instabilité de flottement. Si les valeurs propres deviennent négatives, il s’agit d’une instabilité de divergence. En interaction fluide-structure, on identifie aussi ces deux types d’instabilité : instabilité de flottement lorsque l’amortissement du système s’annule puis devient négatif – il est nécessaire pour cela qu’il y ait un écoulement fluide – instabilité de flambement lorsqu’une pulsation propre du système s’annule.

L’option CRIT_STAB agit en dynamique transitoire non linéaire dans Code_Aster comme pour le cas quasistatique : on mène toujours une analyse de stabilité (donc ne se basant que sur l’étude des matrices de raideur symétriques), pas une analyse d’instabilité dynamique (amortissement devenant négatif).

Avec la modélisation fluide-structure couplée \((u,p,\phi )\) [R4.02.02], il faut modifier la matrice de raideur assemblée (ainsi que la raideur géométrique quand elle est utilisée). Pour cela, il faut renseigner les mot-clés suivants, sous CRIT_STAB:

  • MODI_RIGI = “OUI”,

  • DDL_EXCLUS=(“PHI”,”PRES”,”DH”,).

La liste des degrés de liberté exclus doit comporter tous les types de degrés de liberté liés au modèle fluide: dans l’exemple du cas-test FDNV100 on a donc le potentiel PHI , la pression PRES et le déplacement vertical au niveau de la surface libre DH *. Si l’on ne fait pas ce traitement, alors l’appel à* CRIT_STAB va planter pour cause de matrice singulière et aucune stratégie de décalage ne saurait surmonter cela.

En quasi-statique, ce problème ne se pose pas car la modélisation fluide-structure couplée n’a alors pas de sens.

Ce mot-clé permet de déclencher le calcul, à la fin de chaque incrément de temps, d’un critère de stabilité. Ce critère est utile pour déceler, au cours du chargement, le point à partir duquel on perd la stabilité (par flambage par exemple).

Ce critère est calculé de la façon suivante : à la fin d’un pas de temps, en petites perturbations, on résout \(det\left({K}^{T}+\mu .{K}^{g}\right)=0\) . \({K}^{T}\) est la matrice tangente cohérente à cet instant. \({K}^{g}\) est la matrice de rigidité géométrique, calculée à partir du champ de contraintes à cet instant et :math:`mu ` le coefficient de charge critique.

En pratique, le chargement est instable si \(|\mu |<1\) (en fait \(0<\mu <1\) ). On calcule les valeurs propres par la méthode de Sorensen ( cf. CALC_MODES [U4.52.02]). Ceci peut être assez coûteux pour les problèmes de grande taille.

Le mot-clé CHAR_CRIT permet de gagner du temps en ne faisant qu’un test de Sturm dans la bande de fréquence fournie. Si on trouve au moins une fréquence, alors on calcule réellement les valeurs des charges critiques dans cet intervalle.

Pour les grands déplacements et les grandes déformations, on résout \(det\left({K}^{T}+\mu .I\right)=0\) car \({K}^{T}\) contient alors \({K}^{g}\) .

Le critère est alors un critère d’instabilité : quand :math:`mu ` change de signe (donc passe par 0) le chargement est instable.

Remarques concernant la matrice de raideur géométrique

L’analyse de stabilité sera plus précise si on dispose de la matrice de raideur géométrique : en effet, en son absence, certaines instabilités ne peuvent être détectées comme, par exemple, dans les cas purement élastiques. Entre des modèles suffisamment réalistes, il faut donc privilégier l’utilisation, pour des analyses de stabilité précises, de modèles permettant le calcul de cette raideur géométrique.

Le mot-clé NB_MODE désigne le nombre de charges critiques à calculer. Souvent la première suffit mais il peut y avoir des modes multiples.

On stocke le mode propre correspondant à la plus petite charge critique (en valeur absolue) dans l’objet résultat, sous le nom MODE_FLAMB. Ce mode propre peut être extrait et visualisé (comme un champ de déplacements ou un mode propre classique). Il est normalisé à 1 sur la plus grande composante de déplacement. Toutes les charges critiques calculées sont affichées dans le fichier .mess.

En pratique, afin de limiter le surcoût de calcul, on conseille d’optimiser les appels à CRIT_STAB. On peut pour cela utiliser les mots-clefs LIST_INST/INST/PAS_CALCUL dans le mot-clef facteur CRIT_STAB. On peut ainsi préciser à quels pas de temps on calculera les modes de flambement.

En complément, il est judicieux de n’utiliser CRIT_STAB que sur les intervalles de temps où l’on soupçonne la possibilité d’instabilités.

Enfin, si l’on veut une très bonne évaluation des charges critiques, il convient de bien raffiner le pas de temps à l’approche de cette zone. Ce conseil est d’autant plus pertinent en quasistatique car l’utilisateur a alors souvent recours à des pas de temps plus grands qu’en dynamique.

Il est possible d’arrêter proprement (base réutilisable en poursuite) un calcul avec STAT_NON_LINE ou DYNA_NON_LINE lors de la détection d’une instabilité. Ce n’est pas le fonctionnement par défaut où Code_Aster va tenter de continuer à résoudre le problème, et si on a la convergence, alors cela signifie que l’on a réussi à suivre une des branches de solution.

Pour gérer cet arrêt ( cf. cas-test ERREU10 [V1.01.308]), il faut avoir au préalable utilisé DEFI_LIST_INST avec les arguments suivants :

ECHEC=_F(EVENEMENT=”INSTABILITE”,ACTION=”ARRET”,),

Ce qui signifie qu’en cas d’événement de type instabilité, l’action déclenchée sera l’arrêt.

Sous CRIT_STAB, les options facultatives suivantes permettent de piloter ce critère d’arrêt:

  • PREC_INSTAB pour définir la précision (adimensionnelle) du critère d’arrêt,

  • SIGNE pour spécifier les valeurs critiques à considérer.

Le deuxième mot-clé ne sert que lorsque la matrice de raideur géométrique est utilisée. Par défaut on considère la solution comme étant instable si la charge critique devient comprise entre 1 et -1, mais on peut, au besoin, ne prendre en compte que la partie positive ou négative de cet intervalle.

Sans matrice de raideur géométrique, l’instabilité sera détectée lorsqu’une valeur propre de la matrice de raideur globale assemblée, soit :

  • tendra vers 0 (avec une précision relative donnée par PREC_INSTAB),

  • changera de signe.

Par exemple, dans le cas d’un réservoir rempli d’eau sous séisme, on peut commencer le calcul incrémental ou transitoire avec une charge critique valant 0,8 (analyse avec raideur géométrique) : ce qui signifie que le réservoir flamberait si on imposait une dépression valant 0,8 fois la pression hydrostatique imposée (le valeur positive de la charge critique correspond à une inversion de sens du chargement considéré). Donc si on ne spécifie rien, le calcul serait considéré instable et s’arrêterait. Comme, dans ce cas, on fait l’hypothèse qu’il n’y aura pas dépressurisation (par exemple par vidange), alors on dédouane l’intervalle \([0,1]\) dans l’analyse de stabilité. Donc le problème deviendra instable si la charge critique atteint l’intervalle \([-1,0]\) .

Pour les analyses en évolution monotone ce type de raisonnement se conçoit aisément, ce qui est bien le cas pour la partie statique du chargement, mais pour la partie dynamique, le chargement étant cyclique, et sauf à disposer d’informations spécifiques, il est plus sûr et plus conservatif de garder l’option par défaut et donc de considérer la structure instable si la charge critique devient inférieure à 1 en valeur absolue.

Lors d’un arrêt sur instabilité, le calcul va s’arrêter en fermant proprement la base: l’utilisateur pourra l’exploiter en poursuite.

Suivi de solution instable#

Pour l’étude d’une structure potentiellement instable ou susceptible de connaître un point limite, qui risque donc de rencontrer une bifurcation en solution au cours de l’évolution du chargement, il est souvent utile de pouvoir choisir une branche de solution particulière (souvent la solution physique quand elle est définie a priori sans ambiguïtés). Pour cela, l’utilisateur peut avoir à introduire un défaut initial qui va « forcer » la structure à bifurquer sur la branche de solution particulière.

Plusieurs méthodes existent pour définir ce défaut.

  • L’une des plus adaptée est de pré-déformer légèrement la structure suivant l’allure du mode de flambement correspondant à la branche que l’on veut suivre. L’amplitude de cette pré-déformation doit être faible, par exemple moins de 1/10ème de l’épaisseur pour une structure mince. L’idéal étant de trouver le défaut minimal qui est compatible avec une performance satisfaisante de l’algorithme de résidu en équilibre. En effet, un défaut trop faible peut entraîner une difficulté de convergence du résidu, principalement dans le cas d’un pilotage en effort. Ce défaut initial peut fort judicieusement être construit à partir du mode d’instabilité calculé avec l’option CRIT_STAB de STAT_NON_LINE / DYNA_NON_LINE ( cf. paragraphe précédent). Ce mode tient alors compte de toutes les non-linéarités introduites dans le modèle complet. L’alternative plus économique est d’utiliser le mode de flambement d’Euler, mais qui correspond au cas linéaire.

  • Le défaut géométrique peut aussi être défini par mesures expérimentales de la pièce réelle dont la géométrie ne saurait être parfaite.

  • Le défaut peut aussi prendre la forme d’une perturbation du chargement (défaut d’alignement, rajout d’un chargement localisé, …) ou des caractéristiques mécaniques du matériau (affaiblissement local du module de Young, par exemple). Il peut néanmoins être alors plus difficile d’adapter le défaut au mode de flambage désiré, surtout si la structure présente des modes relativement voisins.

Remarque

Dans certains cas, même sur le problème non perturbé, le chargement est tel qu’il provoque la bifurcation désirée.

Un des autres points particuliers, liés à l’instabilité, est le choix de la technique de pilotage de l’algorithme STAT_NON_LINE, cf. [U4.51.03]. En effet, le pilotage classique en effort n’est plus adapté car il ne peut capter une branche instable de solution. De même, à l’approche d’un point limite, la convergence avec le pilotage en effort deviendra de plus en plus difficile, la matrice de rigidité tangente devenant singulière. Il est alors nécessaire de réduire l’incrément de charge et d’augmenter le nombre maximal d’itération pour continuer le calcul.

On peut aussi servir de la possibilité de s’arrêter proprement en cas d’instabilité ( cf. paragraphe précédent) pour gérer en poursuite la bifurcation sur la branche de solution choisie en initiant la suite du calcul par une perturbation suivant ce mode d’instabilité.

Il existe des techniques de pilotage [bib9] permettant de contourner ces difficultés numériques. Parmi les méthodes proposées par Code_Aster , celle dite par longueur d’arc [bib12] (option TYPE=’LONG_ARC’ du mot clé PILOTAGE dans STAT_NON_LINE), qui est la plus adaptée pour les instabilités de type flambage, dans le cas de snap-backs éventuels « doux » [bib13]. Pour les cas de snap-backs plus brutaux (claquage), M.A.Crisfield propose une variante [bib13], non disponible dans Code_Aster .

D’autres méthodes existent, comme celle de E.Riks [bib14] (non disponible non plus dans Code_Aster), qui traite aussi le cas dynamique.

Si l’on ne veut qu’obtenir le point limite, y compris avec une bonne précision, un pilotage en chargement peut suffire, à condition de bien gérer les paramètres de pas d’incrément de charge (penser à utiliser la commande DEFI_LIST_INST) et de nombre d’itérations maximal autorisé (ITER_GLOB_MAXI de CONVERGENCE). Il peut aussi être utile, à l’approche du point limite, de ne plus employer la matrice tangente réactualisée pour le solveur, puisqu’elle est quasi-singulière. On peut alors se contenter de ne pas réactualiser cette matrice à chaque calcul (paramètres REAC_INCR et REAC_ITER) ou, dans le pire des cas, adopter la matrice élastique de base (PREDICTION=’ELASTIQUE’ et MATRICE=’ELASTIQUE’ du mot clé NEWTON).

Voici un exemple d’utilisation de STAT_NON_LINE pour un calcul élastoplastique en grands déplacements ([bib4] pour les éléments employés, qui sont de type coques volumiques), avec pilotage en efforts :

RESU = STAT_NON_LINE ( …

EXCIT = ( _F( CHARGE = CONDLIM ,

TYPE_CHARGE = “FIXE_CSTE” , ) ,

_F( CHARGE = PESA ,

TYPE_CHARGE = “FIXE_CSTE” , ) ,

_F( CHARGE = PRESPH ,

FONC_MULT = FONCMUL2 ,

TYPE_CHARGE = “SUIV” , ) ,

_F( CHARGE = PRESPS1 ,

FONC_MULT = FONCMUL ,

TYPE_CHARGE = “SUIV” , ) , ) ,

… )

Remarques

  • On utilise la matrice tangente réactualisée à chaque calcul, en autorisant le sous-découpage du pas de charge.

  • Les pressions imposées sont des efforts suiveurs (TYPE_CHARGE=”SUIV”).

  • Dans le cas d’une modélisation en éléments massifs, le tenseur de déformation recommandé en grands déplacements est ‘SIMO_MIEHE’.

Si on veut remplacer le pilotage en effort par une méthode par longueur d’arc, il suffit de rajouter:

RESU = STAT_NON_LINE ( …

PILOTAGE = _F( GROUP_NO = “G” ,

TYPE = “LONG_ARC” ,

NOM_CMP = ( “DY” , ) ,

COEF_MULT = 7. ) ,

…)

Remarques

  • Dans Code_Aster, on ne peut pas piloter de forces suiveuses.

  • Pour le pilotage par longueur d’arc, il est, en général, recommandé que GROUP_NOcontienne toute la structure.

Pour finir, citons deux articles de Crisfield qui donnent une bonne vision générale des problèmes et méthodes liés aux calculs non linéaires pouvant présenter divers types d’instabilités ([bib15] et [bib11]).

La documentation [U2.06.11] montre un exemple d’utilisation de CRIT_STAB pour l’étude du comportement d’un réservoir métallique.

Cas-tests de flambement et de flambage#

Voici quelques cas-tests de Code_Aster traitant du flambement et du flambage :

Charges critiques et modes de flambement d’Euler en statique#

SDLS504

Flambement latéral d’une poutre (déversement)

référence analytique

[V2.03.504]

SDLS505

Flambement d’une enveloppe cylindrique sous pression externe

référence analytique

[V2.03.505]

SSLL103

Flambement élastique d’une cornière

référence analytique

[V3.01.103]

SSLL105

Flambement élastique d’une structure en L

référence avec différents codes

[V3.01.105]

SSLL403

Flambement d’une poutre sous l’effet de son poids propre

référence analytique

[V3.01.403]

SSLL404

Flambement d’une arche

référence analytique

[V3.01.404]

SSLS110

Stabilité d’une plaque carrée comprimée

référence analytique

[V3.03.110]

SSLS125

Flambement d’un cylindre libre sous pression externe

référence analytique, voir aussi SDLS505

[V3.03.125]

SSLV139

Flambement d’une plaque circulaire soumise à une force de compression uniformément répartie sur son contour

référence analytique

[V3.04.139]

Charges critiques et modes de flambement d’Euler en statique non linéaire#

SSNL123

Flambement d’une poutre Multi-Fibres en compression

référence analytique

[V6.02.123]

SSNL126

Flambement élastoplastique d’une poutre droite en compression

référence analytique

[V6.02.126]

Stabilité, flambage et post-flambage en statique non linéaire#

Stabilité, flambage en dynamique non linéaire#

SDNV106

Analyse aux valeurs propres dans DYNA_NON_LINE

stabilité et modes vibratoires

[V5.03.106]

SDND105

Choc d’un point matériel contre une paroi avec flambage plastique

DYNA_VIBRA + mot clé FLAMBAGE et DYNA_NON_LINE + comportement CHOC_ENDO_PENA sur un élément discret ; référence analytique

[V5.01.105]

Bibliographie#

  1. Efforts extérieurs de pression en grands déplacements [R3.03.04].

  2. Pression suiveuse pour les éléments de coques volumiques [R3.03.07].

  3. Éléments finis de coques volumiques [R3.07.04].

  4. Éléments de coques volumiques en non linéaire géométrique [R3.07.05].

  5. Algorithme de résolution pour le problème généralisé [R5.01.01].

  6. Paramètres modaux et norme des vecteurs propres [R5.01.03].

  7. Algorithme non linéaire quasi statique [R5.03.01].

  8. Intégration des relations de comportement élasto-plastique de Von Mises [R5.03.02].

  9. Méthodes de pilotage du chargement [R5.03.80].

  10. Méthode multifrontale [R6.02.02].

  11. M.A. Crisfield, G. Jelenic, Y. Mi, H.-G. ZhongG & Z. Fan : S ome aspects of the non‑linear finite element method, Finite Elements in Analysis and Design, Vol. 27, 19-40, 1997.

  12. M.A. Crisfield : A fast incremental iterative solution procedure that handles snap through , Computers & Structures, Vol. 13, 55-62, 1981.

  13. H.-B. Hellweg & M.A. Crisfield : A new arc-length method for handling sharp snap-backs , Computers & Structures, Vol. 66, 705-709, 1998.

    1. Riks, C.C. Rankin & F.A. Brogan : On the solution of mode jumping phenomena in thin-walled shell structures , Comp. Meth. In Applied Mech. And Engrg., Vol. 1367, 59-92, 1996.

    1. Shi & M.A. Crisfield : Combining arc-length and line searches in path-following , Comm. Numer. Meth. Engrg, Vol. 11, 793-803, 1995.