r7.05.01 Critères de stabilité mécanique

Résumé:

Ce document présente les différents critères de stabilité mécanique disponibles dans Code_Aster . On peut les classer suivant deux catégories:

• Le critère de stabilité associé aux systèmes conservatifs, qui se présente comme la généralisation du critère d’Euler basé sur l’analyse de la matrice de raideur globale réactualisée.

• Le critère de stabilité associé aux systèmes dissipatifs, qui doit tenir compte des contraintes d’irréversibilité liées à la dissipation d’énergie.

Ces critères sont utilisés pour distinguer dans les problèmes quasistatiques, les solutions numériques instables issues du calcul d’équilibre effectué dans la méthode des éléments finis (dérivée première de l’énergie nulle mais dérivée seconde négative) des solutions physiques, stables, pour lesquelles la dérivée seconde de l’énergie est positive.

Les critères présentés dans ce document sont directement transposables au cadre de la dynamique, mais comme ils ne tiennent pas compte ni de la matrice de masse ni de celle d’amortissement, on ne peut parler de critère de stabilité dynamique au sens classique (par exemple, d’amortissement devenant négatif ou nul).

Ces critères sont appelés au sein des opérateurs STAT_NON_LINE et DYNA_NON_LINE, pour pouvoir être évalués à chaque pas de la résolution incrémentale quasi-statique ou transitoire dynamique non linéaire.

Stabilité d’un système dissipatif

Définition de la stabilité d’un système dissipatif

Lorsque l’on s’intéresse à des phénomènes dissipatifs (cas des matériaux plastiques ou fragiles…) on ajoute dans l’expression de l’énergie totale Φ un terme représentant la partie dissipée. Les dégradations irréversibles sont alors associées à des quantités, le plus souvent scalaires, comme l’endommagement ou la plasticité. Le critère de flambement présenté précédemment ne tient pas compte de l’irréversibilité. Il est suffisant mais non nécessaire pour justifier de la stabilité [bib19]. Lorsque le critère de flambement est atteint, l’unicité de la solution n’est plus garantie, mais on ne peut pas conclure directement sur la stabilité. On note u la variable d’état dérivant de la partie réversible de l’énergie et α la variable d’état dérivant de la partie irréversible du phénomène mécanique étudié. Le critère de stabilité consiste alors à vérifier la positivité de la dérivée seconde de l’énergie dans la direction de l’accroissement de la variable α [bib20]. Considérons une perturbation admissible ( ν,b≥0 ). Le critère de stabilité revient alors à vérifier que l’on retrouve toujours l’inégalité:

\(\phi (u,\alpha )\le \phi (u+v,\alpha +b)\) éq 2.1‑1

Du point de vue mathématique, cela revient à vérifier que la fonction Φ réalise un minimum local en ( u , α ).

On s’intéresse plus particulièrement au cas des modèles d’endommagement. On considère ( u , α ), un étatde la structure Ωvérifiant l’équilibre :

\(\forall v\in {C}^{0},{\int}_{\Omega}(\frac{\partial \phi }{\partial u}(u,\alpha ).v)\mathit{dx}=0\) éq 2.1‑2

et le critère d’endommagement, tenant compte de l’irréversibilité de la variable d’état α:

\(\frac{\partial \phi }{\partial \alpha }(u,\alpha )\ge 0\mathit{et}\forall \beta \in {C}^{0}\ge 0,{\int}_{\Omega}(\frac{\partial \phi }{\partial \alpha }(u,\alpha ).\beta )\mathit{dx}=0\) éq 2.1‑3

On obtient rapidement, avec un développement de Taylor à l’ordre 2 [bib20], l’équivalence entre le critère 2.1-1 et le critère suivant, écrit sur la dérivée seconde de Φ:

\({D}^{2}\phi (u,\alpha )\ge 0\) éq 2.1‑4

où D2 est l’opérateur de dérivée seconde.

Écriture dans le cadre de la méthode des éléments finis

Du point de vue des éléments finis et en conservant les notations précédentes, ce critère s’écrit comme la positivité du quotient de Rayleigh sous contraintes d’inégalités suivant :

\(\forall (v,b\ge 0)\mathrm{\ne }0,{Q}_{\mathit{rc}}=\frac{{(v,b)}^{t}.{K}^{T}(v,b)}{{(v,b)}^{T}.(v,b)}\ge 0\) éq 2.2‑1

La façon la plus classique de vérifier qu’une fonction est positive consiste à calculer son minimum et à s’assurer de sa positivité. On présente dans la partie qui suit, l’algorithme de minimisation sous contraintes d’inégalités programmé dans Code_Aster.

Algorithme d’optimisation sous contraintes d’inégalités

Parmi les algorithmes disponibles dans la littérature, celui qui s’est avéré le plus robuste et également le plus simple à programmer est la méthode des puissances, à laquelle on ajoute, à chaque itération, la projection des degrés de liberté \(b\) sur l’ensemble des valeurs positives, de sorte à vérifier la contrainte unilatérale imposée [bib21]. L’algorithme est construit en deux étapes, de la façon suivante:

1. La méthode des puissances est très utilisée en mathématique appliquée pour la recherche des valeurs propres maximales d’une matrice M . Un décalage sur la plus grande valeur propre λm :N=* **λmId-M permet ensuite d’effectuer la recherche des plus petits modes propres. On présente ci-dessous l’algorithme sous sa forme initiale. Soit, sous la forme utilisée pour la recherche de la valeur propre maximale d’une matrice symétrique A:

   \(\begin{array}{c}\text{Soit}P\text{la projection d'un vecteur sur sa partie positive.}\\ \text{Initalisation :}{(v,b\ge 0)}_{0}/\mathrm{\parallel }{(v,b\ge 0)}_{0}\mathrm{\parallel }=1.\\ \text{Pour}k\ge 1\text{}\text{Faire :}\\ \text{}(w,c)=A{(v,b)}_{k-1},\\ \text{}{(v,b)}_{k}=(w,P(c)),\\ \text{}{(v,b)}_{k}={(v,b)}_{k}/\mathrm{\parallel }{(v,b)}_{k}\mathrm{\parallel }\\ \text{}\text{Si :}\mathrm{\parallel }{(v,b)}_{k}-{(v,b)}_{k-1}\mathrm{\parallel }<\epsilon \text{alors Fin}\\ \text{}\text{Sinon :}k=k+1\\ \text{}\text{Fin Si}\\ \text{Fin Pour}\end{array}\)

Figure 2.3-a : Schéma de l’algorithme de recherche du maximum sous contrainte pour une matrice A

1. Sa limite de fiabilité se situe aux alentours du millier de degrés de liberté. Pour s’affranchir de cette limite et être en mesure de traiter des problèmes industriels, on utilise la méthode de réduction disponible dans Sorensen [bib22], qui consiste à projeter le problème sur une base constituée des n plus petits modes propres de la structure. En tenant compte de cette étape on récrit l’algorithme de la façon suivante, où K est toujours l’opérateur tangent, Q est l’opérateur de projection et B n la projection dans l’espace réduit :

   \(\begin{array}{c}\text{Réduction :}{K}^{T}={Q}^{T}{B}_{n}Q\\ \text{Initialisation :}{(v,b\ge 0)}_{0};\text{}{z}_{0}=Q{(v,b)}_{0}/(Q\mathrm{\parallel }{(v,b)}_{0}\mathrm{\parallel })=1,\\ \text{Pour}k\ge 1\text{}\text{Faire :}\\ \text{}(w,c)={Q}^{T}{B}_{n}{z}_{k-1},\\ \text{}{(v,b)}_{k}=(w,P(c)),\\ \text{}{z}_{k}=Q{(v,b)}_{k},\\ \text{}{z}_{k}={z}_{k}/\mathrm{\parallel }{z}_{k}\mathrm{\parallel }\\ \text{}\text{Si :}\mathrm{\parallel }{z}_{k}-{z}_{k-1}\mathrm{\parallel }<\epsilon \text{alors Fin}\\ \text{}\text{Sinon déflation}\\ \text{}\text{Fin si}\\ \text{Fin pour}\end{array}\)

Figure 2.3-b : Schéma de l’algorithme de recherche de maximum sous contraintes d’inégalité avec méthode de projection

Implémentation dans le code

L’étude de stabilité est déclenchée à l’appel de la commande CRIT_STAB de l’opérateur STAT_NON_LINE, sous la condition TYPE = “STABILITE”. Les grandeurs, ou degrés de liberté, vérifiant la condition unilatérale d’irréversibilité sont déclarés dans une liste DDL_STAB=(“”,””,…).

Exemple d’application : Cas de la barre en traction uniforme

Le principale cas de référence retrouvé dans la bibliographie est l’étude de la stabilité de la solution homogène d’une barre endommagée sous l’effet d’un chargement de traction uniforme [bib23] :

Figure 2.5-a: Représentation de la barre en traction uniforme

Résultats analytiques de stabilité

On démontre [bib23] qu’il existe deux types de solutions au problème étudié :

• La solution homogène, uniformément endommagée.

• Les solutions localisées qui concentrent l’endommagement de la barre sur une zone précise.

L’étude de stabilité de la solution homogène revient donc à vérifier que les niveaux d’énergie atteints, en considérant de petites perturbations de la solution ayant la forme d’une localisation, sont toujours supérieurs à celui obtenu à partir de la solution homogène.

Partant de la formulation en gradient d’endommagement suivante de l’énergie [bib24] :

\(\varphi ={\int}_{\Omega}\left(\frac{1}{2}{\left(1-\alpha \right)}^{2}{E}_{0}\epsilon {(u)}^{2}+\frac{{{\sigma}_{M}}^{2}}{{E}_{0}}\alpha +\frac{{E}_{0}{l}^{2}}{2}\nabla \alpha .\nabla \alpha \right)\mathit{dx}\) éq 2.5.1‑1

la comparaison des niveaux d’énergie permet alors de tracer le diagramme de stabilité de la solution homogène en fonction du rapport entre la longueur de la barre L , la longueur interne du modèle l et du chargement appliqué Ut [bib23] :

Figure 2.5.1-a: Diagramme de stabilité analytique de la barre en traction

Résultats de stabilité obtenus avec Code_Aster

En utilisant l’algorithme d’optimisation développé dans le code (2.3), on retrouve un diagramme de stabilité similaire à 5% près sur le chargement à partir duquel l’instabilité est détectée[bib25]:

Figure 2.5.2-a: Diagramme de stabilité de la barre en traction obtenu avec Code_Aster

Des études paramétriques sur l’influence du nombre de fréquences propres calculées et sur lesquelles on s’appuie pour l’étude de stabilité montre qu’il est nécessaire d’en prendre une trentaine pour des problèmes discrétisés à plus ou moins 100 000 ddls mais qu’il devient important d’en considérer une bonne centaine pour les problèmes de plus grande tailles. Ce qui entraîne des coûts de calcul CPU plus importants, l’algorithme étant gourmand en temps. C’est pourquoi il n’est réellement déclenché que lorsque l’étude de flambement montre que l’unicité de la solution n’est plus assurée.

Lorsque l’algorithme montre que le critère de stabilité n’est plus vérifié (CHAR_STAB négatif), le vecteur minimisant le quotient de Rayleigh (éq 2.2-1) est appelé direction d’instabilité. On le retrouve comme résultat sous la nomination MODE_STAB(par analogie avec MODE_FLAMB pour le critère de flambement). La perturbation de la solution courante par cette direction d’instabilité la déstabilise et permet de bifurquer vers une solution stable. Dans l’exemple présenté ici, la direction d’instabilité est la localisation del’endommagement sur l’une des deux extrémités de la barre.

Conclusion

Code_Aster permet d’effectuer des études de stabilité sur des problèmes dissipatifs tels que les problèmes de plasticité ou d’endommagement. L’algorithme utilisé se base sur la méthode des puissances à laquelle s’ajoute la projection des degrés de liberté, dans l’espace respectant les contraintes unilatérales d’irréversibilité. Son application n’est réellement mise en oeuvre que lorsque le critère d’unicité est violé. Dans le cas contraire, l’unicité est suffisante pour garantir la stabilité. Si l’on détecte à un pas de temps du calcul éléments finis la perte de stabilité de la solution, l’algorithme fournit la direction à ajouter comme perturbation pour retrouver une solution bifurquée stable.

Bibliographie

• [1] C. CHAVANT, A. COMBESCURE, J. DEVOS, A. HOFFMANN, Y. MEZIERE: « Flambage élastique et plastique des coques minces », Cours IPSI, 1982.

• [2] S. DURANT, A. COMBESCURE: « Analyse de bifurcation en grandes déformations élasto‑plastiques : formulation et validation élémentaires », Rapport interne n°190, LMT‑Cachan, 1997.

• [3] N. GREFFET: « Simulation couplée fluide-structure appliquée aux problèmes d’instabilité non linéaire sous écoulement », Thèse de doctorat, LMT, ENS-Cachan, 2001.

• [4] R. HILL : « A general theory for uniqueness and stability in elastic-plastic solids », J. Mech. Phys. Solids, Vol. 6, 236-249, 1958.

• [5] T.J.R. HUGHES, W.K. LIU, I. LEVIT : « Nonlinear dynamic finite element analysis of shells, Nonlinear Finite Element Analysis in Structural Mechanics », W. Wunderlich, E. Stein & K.J.Bathe editors, Berlin, Springer, 151-168, 1981.

• [6] G. LOOSS : « Elementary stability and bifurcation theory », Springer-Verlag, 1990.

• [7] A. LEGER, A. COMBESCURE, M. POTIER-FERRY: « Bifurcation, flambage, stabilité en Mécanique des structures », Cours IPSI, 1998.

• [8] A. LEGER: « Bifurcation, flambage, stabilité en mécanique des structures », Note EDF-DER HI‑74/98/024/0.

• [9] J. SHI : « Computing critical points and secondary paths in nonlinear structural stability analysis by finite element method », Computer & Structures, Vol. 58, n°1, 203-220, 1996.

• [10] J.C. WOHLEVER, T.J. HEALEY : « A group theoretic approach to the global bifurcation analysis of an axially compressed cylindrical shell », Comp. Meth. In Applied Mech. And Engrg., Vol. 122, 315-349, 1995.

• [11] «Efforts extérieurs de pression en grands déplacements» [R3.03.04].

• [12] «Éléments vibro-acoustiques» [R4.02.02].

• [13] «Algorithme de résolution pour le problème généralisé» [R5.01.01].

• [14] «Algorithme de calcul du problème quadratiques de valeurs propres» [R5.01.02].

• [15] «Algorithme non linéaire quasi-statique» [R5.03.01].

• [16] «Intégration des relations élasto-plastiques» [R5.03.02].

• [17] «Modèle de Rousselier en grandes déformations» [R5.03.06].

• [18] «Notice de calcul au flambage» [U2.08.04].

• [19] Q-S. NGUYEN: « Stabilité et mécanique non linéaire », HERMES Science Publications, 2000.

• [20] A. BENALLAL et J-J. MARIGO: « Bifurcation and stability issues in gradient theories with softening », Modeling Simul. Mater. Sci. Eng., 15: 283-195, 2007.

• [21] A. PINTO DA COSTA et A. SEEGER: »Numerical resolution of cone-constrained eigenvalue problems », Computational and Applied Mathematics, 28(1): 37-61, 2009.

• [22] «Solveurs modaux et résolution du problème généralisé (GEP)» [R5.01.01].

• [23] K. PHAM, H. AMOR, J-J. MARIGO et C. MAURINI: « Gradient damage models and their use to approximate brittle fracture », International Journal of Damage Mechanics, vol. 20, no. 4: 618-652, 2011.

• [24] «Loi d’endommagemet régularisée quadratique ENDO_CARRE » [R5.03.26]

• [25] «Validation de l’algorithme d’optimisation sous contrainte d’inégalités de l’option DDL_STAB » [V6.02.138].

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
6	N. GREFFET, J.M. PROIX, L. SALMONA EDF- R&D/AMA	Texte initial
04/08/12	N. GREFFET EDF-R&D/AMA,	Généralisation à la dynamique