# Bifurcations primaires et secondaires d'un cylindre sous compression axiale

Philippe Le Grognec, Anh Le Van

► To cite this version:

Philippe Le Grognec, Anh Le Van. Bifurcations primaires et secondaires d'un cylindre sous compression axiale. Revue Européenne des Éléments Finis, 2003, 12 (1), pp.7-41. 10.3166/reef.12.7-41. hal-01005292

HAL Id: hal-01005292

https://hal.science/hal-01005292v1

Submitted on 6 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Bifurcations primaires et secondaires d'un cylindre sous compression axiale

Philippe Le Grognec — Anh Le van

Université de Nantes
Laboratoire de Génie Civil de Nantes-Saint Nazaire
2, rue de la Houssinière BP 92208
F-44322 Nantes Cedex 3

{Philippe.Le-Grognec, Anh.Le-van}@physique.univ-nantes.fr

RÉSUMÉ. Cet article présente une modélisation par éléments finis du comportement en post-flambement des structures de type coques minces élastoplastiques en grandes rotations. Une méthode de longueur d'arc généralisée et des méthodes de continuation appropriées ont été mises en œuvre, permettant la gestion automatique des points limites et points de bifurcation. L'application de la modélisation proposée au problème du cylindre sous compression axiale met en évidence, au-delà des modes axisymétriques et des modes « diamant », des modes transitoires correspondant à des branches bifurquées secondaires.

ABSTRACT. An elastoplastic thin shell finite element model with finite rotations is presented in this paper in order to compute the post-buckling behaviour of shell-type structures. A generalized arc-length method and appropriate continuation methods have been implemented to deal with both limit and bifurcation points. A particular attention is devoted to cylindrical shells under axial compression producing axisymmetric and diamond modes, but also transition modes corresponding to secondary bifurcating branches.

MOTS-CLÉS : Bifurcation, flambement, élastoplasticité, cylindres comprimés axialement.

KEYWORDS: Bifurcation, buckling, elastoplasticity, axially compressed cylinders.
# 1. Introduction

L'utilisation courante dans l'industrie des structures de type coque justifie l'attention particulière qu'on leur accorde dans la littérature. On s'intéresse dans ce travail au flambement et au comportement en post-flambement d'une coque cylindrique en régime élastoplastique. Plus précisément, notre étude portera sur la compression axiale de la coque, qui représente, avec la compression radiale [GUS 00], l'un des plus critiques parmi les différents types de sollicitation conduisant à l'instabilité.

Nous commencerons par retracer brièvement au paragraphe 2 les différents travaux dans la littérature concernant le flambement du cylindre sous compression axiale, tant sur le plan théorique que numérique ou expérimental.

Au paragraphe 3, nous explicitons l'équation de bifurcation du problème de cylindre sous compression axiale et nous en déduisons la charge critique et le mode de bifurcation analytiques dans le cas axisymétrique et pour des conditions aux limites simples. La plasticité est décrite par la théorie des matériaux standard généralisés [HAL 75] en évolution quasi-statique, avec la décomposition additive de la déformation totale [GRE 65],[GRE 71] et avec le critère de von Mises et un écrouissage isotrope linéaire.

La mise en œuvre numérique par éléments finis est proposée au paragraphe 4. Elle s'appuie sur une formulation lagrangienne totale en plasticité, avec un écrouissage mixte (isotrope non linéaire et cinématique linéaire). Afin de tenir compte des non-linéarités matérielles et géométriques, nous adoptons des hypothèses cinématiques appropriées permettant de déduire de manière cohérente, à partir du cadre tridimensionnel, un modèle de coque mince capable de décrire à la fois les grands déplacements et les grandes rotations. Après validation du modèle grâce aux résultats théoriques obtenus au paragraphe 3, la résolution numérique permet de connaître le comportement post-critique lointain d'un cylindre sous compression axiale. Des techniques de longueur d'arc et de continuation adaptées aux très fortes non-linéarités rencontrées sont développées, permettant le passage des points limites en force ou en déplacement et le franchissement des points de bifurcation avec la possibilité de bifurquer sur une branche choisie. L'objectif est d'explorer une large gamme de modes de bifurcation du cylindre sous compression axiale en faisant varier la géométrie, les propriétés matérielles et les conditions aux limites : les modes axisymétriques, les modes de type « diamant », et des solutions transitoires appelées modes « patte d'éléphant », correspondant à des branches bifurquées secondaires décrites dans Goto et Zhang [GOT 99].

# 2. Etat de l'art

Un cylindre sous compression axiale a ceci de critique, en commun notamment avec une sphère sous pression externe, qu'il fait l'objet d'un flambement dit instable, c'est-à-dire d'une diminution de la charge juste après la bifurcation. La valeur critique correspondante constitue alors la charge maximale applicable à la structure. Par
ailleurs, la réponse du cylindre est bien plus sensible aux imperfections que dans les cas de flambement stable de poutres ou de plaques comprimées.

Dans la suite, les principaux résultats théoriques, numériques et expérimentaux disponibles dans la littérature seront brièvement présentés, permettant de discuter de l'influence des différents paramètres géométriques et matériels, puis de juger de l'effet des imperfections.

## 2.1. Connaissances théoriques

### 2.1.1. Influence de la loi de comportement

#### 2.1.1.1. Elasticité

Parmi l'ensemble des modèles de coque existants, la théorie de Donnell [DON 76] est la plus fréquemment utilisée pour l'analyse de problèmes de flambement de cylindres sous compression axiale. Elle permet aisément, le plus souvent par le biais de méthodes énergétiques, de calculer la valeur critique du chargement en élasticité pour un mode axisymétrique ou non axisymétrique faisant apparaître simultanément des ondes axiales et circonférencielles. La charge critique élastique, identique pour ces deux types de mode, est fréquemment donnée sous la dimension d'une contrainte :

$$ \sigma_c^e = \frac{Et}{R} \frac{1}{\sqrt{3(1 - \nu^2)}} \quad [1] $$

On distingue en effet, excepté pour les cylindres extrêmement longs qui flambent suivant un mode de poutre d'Euler, deux types de comportement post-critique théorique. L'équation de Donnell régissant de manière découplée la partie radiale $w$ du mode de flambement s'écrit :

$$ \frac{Et^3}{12(1 - \nu^2)} \nabla^8 w + N_{xx} \nabla^4 w_{,xx} + \frac{Et}{R^2} w_{,xxxx} = 0 \quad [2] $$

où $t$ représente l'épaisseur de la coque, $R$ le rayon du cylindre, $N_{xx}$ l'effort axial par unité de longueur de bord, $\nabla$ est l'opérateur biharmonique et $w$ dépend des positions longitudinale $x$ et angulaire $y$ du point considéré.

Elle fournit, dans le cas le plus général représentatif des coques minces, une solution modale non axisymétrique. Pour les coques relativement plus épaisses, un comportement axisymétrique est obtenu, solution alors de l'équation simplifiée (après quatre intégrations en $x$) :

$$ \frac{Et^3}{12(1 - \nu^2)} w_{,xxxx} + N_{xx} w_{,xx} + \frac{Et}{R^2} w = 0 \quad [3] $$

Déjà en élasticité, on met en évidence de larges écarts entre la théorie et les résultats expérimentaux. Le mode de bifurcation observé expérimentalement est souvent
moins régulier que les expressions exhibées en théorie, mais surtout la valeur critique expérimentale chute parfois à moins de 50 % de la valeur analytique. Ces différences, dues aux imperfections inévitable dans le processus expérimental, et particulièrement aux défauts géométriques, peuvent s'expliquer et se quantifier par la théorie de Koiter, tout au moins pour ce qui est de l'élasticité. De manière équivalente à Koiter, l'utilisation d'une imperfection modale (axisymétrique sinusoïdale) certainement la plus défavorable [BUD 72] permet d'établir des coefficients de sécurité pour la charge critique dépendant de la taille de l'imperfection et de juger de la stabilité sur la branche bifurquée. L'étude réalisée par Amazigo et Budiansky [AMA 72] complète ces analyses en discutant de l'influence du caractère local ou global des imperfections.

### 2.1.1.2. Plasticité

Pour une coque assez épaisse, le flambement a toutes les chances de se produire dans le domaine plastique. On va caractériser ce flambement plastique par la nouvelle contrainte critique principalement dans le cas d'un mode axisymétrique.

La diversité des valeurs critiques obtenues jusqu'à ce jour est inhérente à l'adoption de multiples théories de plasticité et en particulier à la confrontation entre la théorie de déformation et la théorie incrémentale. Sans ignorer les nombreux résultats issus de la théorie de déformation, on s'intéresse davantage aux principaux résultats de la théorie incrémentale qui est celle retenue dans notre formulation, et particulièrement à ceux de Batterman développés de manière détaillée dans [BAT 65]. A partir des équations d'équilibre en vitesses, Batterman obtient la contrainte critique :

$$ \sigma_c^p = \frac{Et}{R} \sqrt{\frac{4h}{3[(5 - 4\nu)E + 4h(1 - \nu^2)]}} \quad [4] $$

avec $h$ le module d'écrouissage.

Les résultats obtenus par la théorie de déformation s'écartent assez nettement de cette dernière valeur. Généralement, la théorie de déformation apporte des solutions en meilleur accord avec les résultats expérimentaux, malgré l'absence de prise en compte de la possibilité de décharge contraire à la réalité expérimentale. Ces écarts entre les deux théories et cette inadéquation apparente de la théorie incrémentale à rendre compte du comportement réel de ce type de structure vont de pair avec la sensibilité des résultats aux imperfections. Après avoir remarqué que bon nombre de problèmes exhibaient, déjà en élasticité, une telle sensibilité aux défauts de toutes sortes, on constate que la plasticité ne fait qu'accroître cette sensibilité, voire la fait apparaître pour certains problèmes globalement insensibles en élasticité.

Le cylindre sous compression axiale en plasticité connaît donc une grande sensibilité aux imperfections. On montre que les défauts de forme sont les principaux responsables des erreurs observées, notamment sur les valeurs critiques. Hutchinson [HUT 74], dans son analyse asymptotique du post-flambement élastoplastique, étend les résultats de Koiter en élasticité concernant l'influence de divers types d'imperfections géométriques. Il est naturel que pour les cylindres sous compression axiale,
l'effet des défauts de forme quels qu'ils soient se fasse davantage sentir pour les cylindres les plus minces.

D'autres sources d'erreurs sont reliées aux conditions aux limites expérimentales ou aux chargements imparfaits. Enfin, le dernier type d'erreurs commises peut provenir de l'hétérogénéité du matériau utilisé et des contraintes résiduelles souvent mal connues.

Signalons qu'en présence d'un tel affrontement entre les théories de déformation et incrémentale, Ariaratnam et Dubey [ARI 69] cherchent une explication dans la forme de la surface d'écoulement plastique et choisisissent de conserver la théorie incrémentale en remplaçant le critère de von Mises qui avait pratiquement l'exclusivité jusqu'alors, pour obtenir de meilleures prédicitions de la charge critique. Christoffersen et Hutchinson, en 1979, proposent, quant à eux, une nouvelle théorie de surface de charge à coins (« corner theory ») constituant un compromis entre les deux théories précédentes et qui sera largement utilisée par la suite, lors de développements numériques.

Mentionnons enfin qu'à l'inverse du comportement élastique, la contrainte critique en plasticité correspondant à un mode non axisymétrique est différente des valeurs obtenues pour les modes axisymétriques. Cette nouvelle valeur, obtenue en particulier dans [BAT 66] pour la théorie incrémentale, est toujours supérieure aux précédentes et tend rapidement, pour un module tangent suffisamment petit, vers la valeur suivante :

$$ \sigma_c^{\text{lim}} = \frac{2Et}{3R} \sqrt{\frac{3}{(1+\nu)(5-4\nu)}} \quad [5] $$

### 2.1.2. Influence des conditions aux limites

Les contraintes critiques [1] et [4] ont été données dans le cas d'un appui simple aux extrémités du cylindre. On peut aisément montrer qu'un encastrement des bases doit fournir la même valeur critique et le même type de mode. De façon à compléter la description et la compréhension de l'influence des conditions aux limites, Batterman [BAT 68] étudie analytiquement le cas de bords libres, se plaçant ainsi dans la situation la plus défavorable parmi toutes les conditions aux limites envisageables.

Il obtient alors, pour une coque cylindrique semi-infinie avec un module tangent constant, une contrainte critique diminuée de moitié par rapport au cas classique des appuis simples :

$$ \sigma_c^p = \frac{Et}{R} \sqrt{\frac{h}{3[(5 - 4\nu)E + 4h(1 - \nu^2)]}} \quad [6] $$

Pour un module dépendant de l'état de contraintes, la relation précédente est implicite et fournit une contrainte supérieure à la moitié de la valeur calculée pour un cylindre simplement appuyé. Hoff montre en outre que, pour une coque de longueur finie, le résultat est légèrement inférieur au résultat de la coque semi-infinie, mais ne dépend pratiquement pas de la longueur. Enfin, le mode de flambement axisymétrique
prend la forme d'une onde dont l'amplitude est maximale aux extrémités et de plus en plus amortie au fur et à mesure que l'on s'en éloigne.

### 2.1.3. Bilan

Une analyse plus fine de l'effet des paramètres géométriques et matériels, des conditions aux limites et du mode de chargement sur la réponse non linéaire complète du cylindre sous compression axiale nécessite le recours à un outil numérique de type éléments finis. Les résultats analytiques précédents permettront de valider le programme numérique utilisé dans ce travail. Ensuite, la prise en compte des imperfections et d'autres facteurs dans la modélisation numérique peut nous renseigner davantage sur le comportement réel de ces structures, en particulier jusqu'à un stade post-critique avancé.

Avant de passer aux travaux numériques et expérimentaux disponibles dans la littérature, signalons l'existence d'autres approches théoriques, plus géométriques, traitant certains aspects du comportement d'un cylindre sous compression axiale.

D'un point de vue mathématique, la déformée post-critique peut être vue comme une transformation isométrique de la géométrie initiale du cylindre. On peut de cette façon, sans connaître l'évolution précise de la géométrie déformée, accéder à la forme flambée avancée et à la charge limite applicable, dans le cas d'un flambement global ou éventuellement localisé.

D'une manière un peu différente, certains travaux, dont [SIN 99] et [GUP 00], entreprennent de représenter le comportement d'un cylindre sous compression axiale par la formation cyclique de plis (engendrant des rotules plastiques) et examinent, dans les deux cas axisymétrique et non axisymétrique, la position et la forme de ces plis successifs. La comparaison des résultats obtenus avec les observations expérimentales est très satisfaisante si on se contente d'une analyse limite du comportement post-critique lointain.

## 2.2. Principaux résultats numériques et expérimentaux antérieurs

Aflak [AFL 88] a mené une étude expérimentale poussée doublée d'une analyse numérique du problème, couvrant divers types d'imperfections et de conditions aux limites, et il en déduit des règles générales de dimensionnement pour les cylindres sous compression axiale. Plus tard, Limam [LIM 91] s'intéresse plus particulièrement à l'étude de la combinaison de chargements (compression/pression interne). Il montre expérimentalement que la pression interne dans le cylindre a un effet stabilisateur et conduit à une meilleure adéquation des résultats avec la théorie. De plus, il met l'accent sur le caractère plus néfaste d'une imperfection modale non axisymétrique par rapport à un défaut axisymétrique. Enfin, il traite le cas plus complexe de coques cylindriques raidies.
Dans le même esprit, Mandal et Calladine [MAN 00], Lancaster, Calladine et Palmer [LAN 00] ont montré, grâce à un procédé expérimental adapté, que la présence d'imperfections sous la forme de contraintes localisées est responsable en grande partie de l'abaissement de la valeur critique. En effet, la relative liberté en déplacement longitudinal des extrémités lors du flambement de cylindres sous leur poids propre diminue considérablement la sensibilité à ce type d'imperfections.

Les études précédentes montrent bien le rôle primordial des imperfections dans le comportement d'un cylindre sous compression axiale et surtout dans la valeur critique de flambement associée. Cependant, nous écartons dans notre étude la prise en compte de ces imperfections pour ne porter notre attention que sur l'analyse du comportement post-critique lointain du cylindre parfait.

En théorie, on suppose très souvent une loi de comportement linéarisée afin de simplifier la résolution des équations, alors qu'une telle loi peut s'avérer impropre même à la description du comportement pré-critique. Ce n'est que récemment que Peek [PEE 00] détermine analytiquement le comportement post-critique immédiat du cylindre sous compression axiale en plasticité finie (et en axisymétrie). Il remarque que la plasticité finie, si elle est prise en compte dans la solution pré-critique, a un rôle stabilisateur. Le premier point de bifurcation obtenu est suivi (de très près) d'un point limite sur la branche bifurquée, correspondant au maximum de la charge applicable, si bien que cette bifurcation perd le caractère instable qu'on lui affectait jusqu'alors.

En fait, de nombreux résultats numériques ont déjà anticipé cette solution et offrent également une panoplie d'informations qui s'accordent souvent avec les résultats d'es sais expérimentaux. La concordance est le plus souvent qualitative dans la mesure où les résultats numériques mettent en valeur l'influence d'un (ou de plusieurs) para mètre(s) étudié(s) conformément aux observations expérimentales.

Dans cet esprit, Tvergaard [TVE 83] se base sur des résultats expérimentaux an térieurs (issus des travaux de Batterman, Lee, Horton et Johnson dans les années soixante) pour valider son approche numérique basée sur la théorie phénoménologique de Christoffersen et Hutchinson. Il s'appuie sur une géométrie initiale comportant un défaut axisymétrique sinusoïdal et cherche l'éventuelle occurence d'une bifurcation vers un mode non axisymétrique. Il retrouve, pour les coques assez minces, la pré sence d'un point de bifurcation juste avant le maximum de la charge. En revanche, pour les coques plus épaisses, la charge maximale supportée par la structure coïncide avec le phénomène de localisation de la déformée axisymétrique aux deux extrémités du cylindre en deux demi-ondes axiales extérieures. La bifurcation est donc retardée ou peut ne pas avoir lieu (les coques très épaisses tombent en effet en ruine axisymé trique). Quand la bifurcation précède la localisation, on parle de mode « diamant » en raison de la présence simultanée d'ondes axiales et circonférencielles sur l'ensemble de la structure. Dans le cas contraire, on parle de mode « patte d'éléphant », la perte d'axisymétrie suivant la localisation donnant naissance à des « doigts » régulièrement répartis sur les deux extrémités du cylindre. Tvergaard explique le retard de la bifur cation par le raidissement de la structure dû à la localisation. Il montre en outre l'effet du rapport géométrique du rayon sur l'épaisseur $R/t$ qui, lorsque celui-ci diminue, re-
tarde la bifurcation et abaisse l'ordre circonférenciel de la solution. Enfin, il constate qu'un faible écrouissage favorise la ruine axisymétrique, ou tout au moins l'obtention d'un mode « patte d'éléphant ».

Récemment, l'analyse de cette transition entre mode axisymétrique et mode « diamant » est menée sur un cylindre parfait et avec une loi élastoplastique incrémentale par Goto et Zhang [GOT 99]. Les conditions aux limites s'avèrent déterminantes pour obtenir les différents types de modes précédents. Goto et Zhang retrouvent, dans le cas du cylindre encastré aux extrémités, et pour une épaisseur moyenne, la succession de deux points de bifurcation, le premier (sur la branche fondamentale) au voisinage du point limite et le second (sur la branche primaire à charge décroissante) débouchant sur un mode localisé non axisymétrique. Plus le rapport $R/t$ est élevé, plus ce mode secondaire « patte d'éléphant » survient tôt. A la limite, pour un cylindre suffisamment mince, le second point de bifurcation remplace le premier et empêche toute localisation, donnant lieu ainsi à un mode « diamant ». Inversement, pour une épaisseur trop importante, le second point de bifurcation n'apparaît pas, tout au moins dans le domaine de prospection, et la ruine a toutes les chances d'être axisymétrique. L'augmentation du rapport $R/t$ entraîne l'accroissement du nombre d'ondes circonférencielles. De façon moins évidente, le rapport $l/R$ de la longueur sur le rayon, lorsqu'il augmente, avance également la seconde bifurcation sans pour autant modifier l'ordre circonférenciel de la déformée. Dans le cas d'extrémités simplement appuyées, la localisation, si elle a lieu, ne conduit jamais à une seconde bifurcation, compte tenu de la position extrême de l'onde localisée. Enfin, la condition de bords libres produit une seule bifurcation non axisymétrique directement localisée. L'ordre circonférenciel de la solution augmente encore naturellement avec le rapport $R/t$.

## 3. Solution analytique du flambement élastoplastique d'un cylindre sous compression axiale

Dans ce paragraphe, on applique au cylindre sous compression axiale l'analyse tridimensionnelle de la bifurcation et de la post-bifurcation élastoplastiques, basée sur l'utilisation de développements asymptotiques.

### 3.1. Théorie de la bifurcation plastique

Rappelons d'abord les principaux résultats de la théorie de bifurcation plastique dans le cadre tridimensionnel. Beaucoup de formulations ont été proposées en matière de plasticité en grandes transformations et il n'est pas facile d'en choisir une entre toutes dans la mesure où chacune a ses propres avantages et inconvénients. En fait, le sujet suscite aujourd'hui encore beaucoup de débats aussi bien sur la cinématique, sur les lois de comportement que sur les algorithmes numériques. Dans ce travail, nous choisisons une formulation lagrangienne totale basée sur la décomposition additive de Green-Naghdi [GRE 65],[GRE 71] et sur la théorie des matériaux standard généralisés [HAL 75].
Considérons un corps tridimensionnel occupant le domaine de référence Ω à l'instant initial et soumis aux grandes transformations quasi-statiques dans le domaine élastoplastique. Nous allons écrire les équations du problème en termes de variables lagrangiennes relatives à cette configuration. Le matériau est caractérisé en élasticité par (E,ν) ou (λ,μ). La plasticité est détectée par le critère de von Mises avec un écrouissage mixte (isotrope non linéaire et cinématique linéaire) :

$$f(\bar{\bar{\Sigma}}, \bar{\bar{A}}, A) = \sqrt{\frac{3}{2}(\bar{\bar{\Sigma}}^d - \bar{\bar{A}}) : (\bar{\bar{\Sigma}}^d - \bar{\bar{A}}) - \sigma_0 - A} \quad [7]$$

où $\bar{\bar{\Sigma}}^d$ désigne le déviateur du second tenseur des contraintes de Piola-Kirchhoff $\bar{\bar{\Sigma}}$ (symétrique) et $\sigma_0$ est la contrainte limite élastique. Les forces thermodynamiques sont représentées ici, à titre d'illustration, par un scalaire noté $A$ et un tenseur d'ordre deux noté $\bar{\bar{A}}$ supposé être un déviateur pur.

En régime plastique, la déformation de Green-Lagrange est décomposée de manière additive selon Green et Naghdi [GRE 65],[GRE 71] :

$$\bar{\bar{E}} = \bar{\bar{E}}^e + \bar{\bar{E}}^p \quad [8]$$

$\bar{\bar{E}}^e$ (resp. $\bar{\bar{E}}^p$) est appelée de façon formelle la déformation élastique (resp. plastique). En fait, ces tenseurs ne sont pas définis cinématiquement en fonction du déplacement, puisque $\bar{\bar{E}}^p$ sera définie par la loi d'évolution [11] alors que $\bar{\bar{E}}^e$ l'est simplement par [8], réalisant par construction la décomposition additive de la déformation. Comme l'ont souligné Green et Naghdi [GRE 65], le tenseur « élastique » $\bar{\bar{E}}^e$ n'est pas un tenseur de déformation élastique classique. S'il est ainsi appelé, c'est uniquement par analogie avec le cas des petites transformations. Dans cette théorie, seules les déformations totale et plastique $\bar{\bar{E}}$ et $\bar{\bar{E}}^p$ sont considérées comme des variables indépendantes, la troisième déformation $\bar{\bar{E}}^e$ étant liée aux deux premières par la relation de décomposition additive [NAG 75],[CAS 80],[NAG 90].

Comme dans [GRE 65] et [CAS 80], on suppose d'emblée que la déformation plastique $\bar{\bar{E}}^p$ est symétrique et objective. Par conséquent, $\bar{\bar{E}}^e$ l'est aussi.

### REMARQUE.
— La décomposition cinématique constitue la première étape de toute étude en plasticité finie, elle ne peut se faire qu'au prix de certaines hypothèses physiques faites a priori sur le comportement de l'élément de matière étudié au cours de la déformation et reste à l'heure actuelle un problème ouvert. De nombreuses décompositions cinématiques ont été proposées dans la littérature, mais notre objectif n'est pas d'en élire la meilleure si tant est que la tâche soit possible. Dans ce travail, nous avons retenu la décomposition additive [8] de Green-Naghdi, dans la mesure où elle permet d'aboutir à une formulation analogue à celle en petites transformations et d'étendre à moindre coût les algorithmes connus dans ce cadre. La décomposition la plus couramment utilisée dans la littérature est celle multiplicative du gradient de la transformation de Lee [LEE 67],[LEE 69] utilisant la notion de configuration intermédiaire relâchée, mais il en existe d'autres telles que la théorie de Besseling et
van der Giessen [GIE 89a],[GIE 89b],[BES 94] utilisant la notion de configuration de référence naturelle, la décomposition additive de la vitesse eulérienne de déformation de Nemat-Nasser [NEM 79],[NEM 82]. Green et Naghdi ont essayé de concilier le point de vue de Lee et le leur dans [GRE 71]. Voir aussi une synthèse des différentes formulations sur ce sujet dans [NAG 90].

Avec la décomposition additive [8], le développement de la thermodynamique est analogue à celui de la plasticité en petites transformations. Nous supposons que l'énergie libre volumique W s'écrit comme la somme de l'énergie élastique $W^e$, fonction quadratique de la déformation élastique $\bar{\bar{E}}^e$, et de l'énergie d'écrouissage $W^\alpha$, fonction des variables d'écrouissage $\bar{\bar{\alpha}}$ et $\alpha$, conjuguées respectives des forces thermodynamiques $\bar{\bar{A}}$ et $A$ :

$$W(\bar{\bar{E}}^e, \bar{\bar{\alpha}}, \alpha) = W^e(\bar{\bar{E}}^e) + W^\alpha(\bar{\bar{\alpha}}, \alpha) = \frac{1}{2}\bar{\bar{E}}^e : \bar{\bar{D}}^e : \bar{\bar{E}}^e + W^\alpha(\bar{\bar{\alpha}}, \alpha) \quad [9]$$

où $\bar{\bar{D}}^e$ désigne le tenseur d'ordre quatre d'élasticité supposé constant. Les lois d'état expriment le second tenseur des contraintes de Piola-Kirchhoff $\bar{\bar{\Sigma}}$ et les forces thermodynamiques $\bar{\bar{A}}$ et $A$ comme les dérivées de l'énergie libre par rapport à la déformation élastique et aux variables d'écrouissage, respectivement :

$$\bar{\bar{\Sigma}} = \frac{\partial W}{\partial \bar{\bar{E}}^e} = \bar{\bar{D}}^e : \bar{\bar{E}}^e \qquad \bar{\bar{A}} = \frac{\partial W}{\partial \bar{\bar{\alpha}}} \qquad A = \frac{\partial W}{\partial \alpha} \quad [10]$$

Si les déformations élastiques restent petites durant tout le processus, il est licite de décrire la loi de comportement élastique par le modèle de Saint-Venant-Kirchhoff : $D^e_{ijkl} = \lambda\delta_{ij}\delta_{kl} + \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj})$ ($\delta_{ij}$ : symbole de Kronecker).

En théorie des matériaux standard généralisés, la loi de normalité exprime que la déformation plastique $\bar{\bar{E}}^p$ et les variables d'écrouissage $\bar{\bar{\alpha}}$ et $\alpha$ dérivent d'un potentiel de dissipation qui n'est autre que la fonction indicatrice du domaine d'élasticité dans l'espace des contraintes généralisées $(\bar{\bar{\Sigma}}, \bar{\bar{A}}, A)$. On en déduit alors les lois d'évolution :

$$\dot{\bar{\bar{E}}}^p = \dot{\Lambda} \frac{\partial f}{\partial \bar{\bar{\Sigma}}} \qquad \dot{\bar{\bar{\alpha}}} = -\dot{\Lambda} \frac{\partial f}{\partial \bar{\bar{A}}} \qquad \dot{\alpha} = -\dot{\Lambda} \frac{\partial f}{\partial A} \quad [11]$$

où $\dot{\Lambda}$ désigne le multiplicateur plastique. Les variables $\bar{\bar{E}}^p$, $\bar{\bar{\alpha}}$ et $\alpha$ étant lagrangiennes, leurs vitesses sont des dérivées ordinaires par rapport au temps cinématique.

La condition de consistance fournit l'expression du multiplicateur plastique :

$$\dot{\Lambda} = \frac{\frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e : \dot{\bar{\bar{E}}}^e}{\bar{\bar{h}} + \frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}}} \quad [12]$$
où $\bar{h} = \frac{\partial f}{\partial \bar{A}} : \frac{\partial^2 W}{\partial \bar{\alpha}^2} : \frac{\partial f}{\partial \bar{A}} + \left(\frac{\partial f}{\partial A}\right)^2 \frac{\partial^2 W}{\partial \alpha^2}$ est le module d'écrouissage global. On en déduit l'expression du tenseur tangent élastoplastique matériel :

$$ \bar{\bar{D}}^p = \frac{\partial \bar{\bar{\Sigma}}}{\partial \bar{\bar{E}}} = \bar{\bar{D}}^e - \frac{\bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}} \otimes \frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e}{\bar{h} + \frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}}} \quad [13] $$

que l'on écrit encore sous la forme :

$$ \bar{\bar{D}}^p = \bar{\bar{D}}^e - \bar{N} \otimes \bar{N} \quad [14] $$

avec :

$$ \bar{N} = \frac{\bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}} \quad \bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e}{\sqrt{\bar{h} + \frac{\partial f}{\partial \bar{\bar{\Sigma}}} : \bar{\bar{D}}^e : \frac{\partial f}{\partial \bar{\bar{\Sigma}}}}} \quad [15] $$

A partir du tenseur tangent élastoplastique matériel [13], on peut construire le tenseur nominal :

$$ \bar{\bar{K}}^p = \frac{\partial \bar{\bar{\Pi}}}{\partial \bar{\bar{F}}} = \bar{\bar{F}} \cdot \frac{\partial \bar{\bar{\Sigma}}}{\partial \bar{\bar{E}}} \cdot \bar{\bar{F}}^T + \left( \bar{\bar{I}} \cdot \bar{\bar{\Sigma}} \right)^T \quad [16] $$

Dans la relation précédente, $\bar{\bar{F}}$ est le gradient de la transformation, $\bar{\bar{\Pi}}$ le premier tenseur des contraintes de Piola-Kirchhoff (non symétrique), $\bar{\bar{I}}$ désigne le tenseur unitaire d'ordre quatre et l'exposant $T$ la transposition majeure d'un tenseur d'ordre quatre. De manière analogue à [14], on peut mettre la dernière expression sous la forme :

$$ \bar{\bar{K}}^p = \bar{\bar{K}}^e - \bar{\bar{M}}^T \otimes \bar{\bar{M}} \quad [17] $$

avec :

$$ \bar{\bar{K}}^e = \bar{\bar{F}} \cdot \bar{\bar{D}}^e \cdot \bar{\bar{F}}^T + \left( \bar{\bar{I}} \cdot \bar{\bar{\Sigma}} \right)^T \qquad \bar{\bar{M}} = \bar{\bar{N}} \cdot \bar{\bar{F}}^T \quad [18] $$

On suppose alors qu'on connaît une courbe d'équilibre fondamentale $\lambda \mapsto \vec{u}_f(\lambda)$, solution en déplacement du problème élastoplastique sous chargement $\lambda$ croissant ($\dot{\lambda} > 0$), et qu'à un instant critique $t_c$, auquel correspond $\lambda_c = \lambda(t_c)$, il existe une solution bifurquée notée $\vec{u}$ (il n'y a pas de confusion possible entre le paramètre $\lambda$, facteur de charge représentant l'amplitude de la sollicitation appliquée à la structure, et la constante $\lambda$ de Lamé en élasticité !).

On étudie cette branche bifurquée par son développement asymptotique :

$$ \begin{gathered} \lambda = \lambda_c + \lambda_1 \xi + o(\xi) \\ \vec{u} = \vec{u}_f(\lambda) + \xi \vec{X} + o(\xi) \end{gathered} \quad [19] $$
Sous les hypothèses du module tangent et de valeur propre simple, le chargement critique $\lambda_c$ et le mode de bifurcation $\vec{X}$ sont donnés par l'équation de bifurcation [HUT 72],[NGU 00] :

$$ \forall \delta \vec{u}, \quad \int_{\Omega} \bar{\nabla}^T \delta \vec{u} : \bar{\bar{K}}^i (\vec{u}_f(\lambda_c)) : \bar{\nabla} \vec{X} d\Omega = 0 \quad [20] $$

où $\bar{\bar{K}}^i = \bar{\bar{K}}^e$ ou $\bar{\bar{K}}^p$ suivant qu'on se trouve respectivement dans le domaine élastique ou dans la zone plastique.

Dans le cas d'une bifurcation avec décharge naissante (situation correspondante à l'hypothèse dite de charge totale), on peut déterminer la valeur de $\lambda_1$ notée alors $\lambda_m$ :

$$ \begin{aligned} \lambda_m &= \min\{\lambda \mid \forall \vec{x} \in \Omega_f^p(t_c), \lambda \bar{\bar{M}}_c : \bar{\nabla} \vec{u}_f, \lambda(\lambda_c) + \bar{\bar{M}}_c : \bar{\nabla} \vec{X} \ge 0 \} \\ &= \max_{\vec{x} \in \Omega_f^p(t_c)} \left( - \frac{\bar{\bar{M}}_c : \bar{\nabla} \vec{X}}{\bar{\bar{M}}_c : \bar{\nabla} \vec{u}_f, \lambda(\lambda_c)} \right) \end{aligned} \quad [21] $$

où $\Omega_f^p(t_c)$ désigne la zone plastique à l'instant critique $t_c$ de la solution fondamentale.

## 3.2. Application au cas du cylindre sous compression simple

On se place maintenant dans le cas particulier d'un état de compression simple, de façon à pouvoir expliciter les composantes du tenseur tangent matériel. On considère un cylindre circulaire de longueur $l$, de rayon $R$ et d'épaisseur $t$. L'une des extrémités du cylindre reste en contact avec le plan $x = 0$ (les déplacements radiaux et les rotations sont libres) tandis qu'en $x = l$ s'applique une contrainte de compression axiale $\Sigma$ répartie uniformément. On prendra $\lambda = -\Sigma > 0$. Enfin, l'écrouissage est pris purement isotrope linéaire, de module $h$ constant (faire, dans [7], $A = h\alpha$).

### 3.2.1. Solution fondamentale

En supposant les déformations pré-critiques suffisamment petites, hypothèse vérifiée en pratique, on peut simplifier l'expression [16] en y remplaçant le tenseur $\bar{\bar{F}}$ par le tenseur identité :

$$ \frac{\partial \bar{\bar{\Pi}}}{\partial \bar{\bar{F}}} = \frac{\partial \bar{\bar{\Sigma}}}{\partial \bar{\bar{E}}} + \left( \bar{\bar{I}} \cdot \bar{\bar{\Sigma}} \right)^T \quad [22] $$

ce qui revient à écrire le tenseur tangent élastoplastique nominal de la manière suivante :

$$ \bar{\bar{K}}^p = \bar{\bar{D}}^p + \left( \bar{\bar{I}} \cdot \bar{\bar{\Sigma}} \right)^T \quad [23] $$

De même, on peut confondre les tenseurs $\bar{\bar{M}}$ et $\bar{\bar{N}}$ dans [18] :

$$ \bar{\bar{M}}(= \bar{\bar{N}} \cdot \bar{\bar{F}}^T) = \bar{\bar{N}} \quad [24] $$
Le tenseur des contraintes $\bar{\bar{\Sigma}}$ s'écrit dans une base orthonormée appropriée :

$$ \bar{\bar{\Sigma}} = \begin{bmatrix} \Sigma & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad [25] $$

Son déviateur prend la forme :

$$ \bar{\bar{\Sigma}}^d = \begin{bmatrix} \frac{2}{3}\Sigma & 0 & 0 \\ 0 & -\frac{1}{3}\Sigma & 0 \\ 0 & 0 & -\frac{1}{3}\Sigma \end{bmatrix} \Rightarrow \left\| \bar{\bar{\Sigma}}^d \right\|^2 = \frac{2}{3}\Sigma^2 \qquad [26] $$

On déduit alors de [13] le tenseur tangent élastoplastique matériel :

$$ \bar{\bar{D}}^p = \bar{\bar{D}}^e - \frac{6\mu^2 \bar{\bar{\Sigma}}^d \otimes \bar{\bar{\Sigma}}^d}{h + 3\mu \left\| \bar{\bar{\Sigma}}^d \right\|^2} \qquad [27] $$

soit :

$$ \begin{aligned} D_{1111}^p &= \lambda + 2\mu - \frac{4\mu^2}{h+3\mu} \\ D_{2222}^p &= D_{3333}^p = \lambda + 2\mu - \frac{\mu^2}{h+3\mu} \\ D_{1122}^p &= D_{1133}^p = \lambda + \frac{2\mu^2}{h+3\mu} \\ D_{2233}^p &= \lambda - \frac{\mu^2}{h+3\mu} \\ D_{1212}^p &= D_{1313}^p = D_{2323}^p = \mu \end{aligned} \qquad [28] $$

Les autres termes sont nuls ou se déduisent de [28] par les symétries majeure et mineures du tenseur $\bar{\bar{D}}^p$.

Dans le cas particulier de la compression simple ($\Sigma < 0$), l'écriture du tenseur $\bar{\bar{N}}$ se simplifie comme suit :

$$ \bar{\bar{N}} = \frac{\mu}{\sqrt{h + 3\mu}} \begin{bmatrix} -2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \qquad [29] $$

D'une part, l'état de compression simple [25] produit respectivement les états de déformation totale et de déformation plastique suivants :

$$ \bar{\bar{E}} = \begin{bmatrix} E_{xx} & 0 & 0 \\ 0 & E_{\theta\theta} & 0 \\ 0 & 0 & E_{rr}(= E_{\theta\theta}) \end{bmatrix} \qquad \bar{\bar{E}}^p = \begin{bmatrix} E_{xx}^p & 0 & 0 \\ 0 & -\frac{1}{2}E_{xx}^p & 0 \\ 0 & 0 & -\frac{1}{2}E_{xx}^p \end{bmatrix} \qquad [30] $$

dans la base cylindrique $(x, \theta, r)$, où :

$$ \begin{aligned} E_{xx} &= -\frac{\lambda}{E} - \frac{\lambda - \sigma_0}{h} \\ E_{\theta\theta} &= \frac{\nu\lambda}{E} + \frac{\lambda - \sigma_0}{2h} \end{aligned} \qquad [31] $$
D'autre part, le déplacement d'un point courant du cylindre s'exprime dans la base cylindrique en fonction des déplacements longitudinal $u$, orthoradial $v$ et radial $w$ d'un point de la surface moyenne et d'un paramètre $\rho$ décrivant la position dans l'épaisseur $(-\frac{t}{2} \le \rho \le \frac{t}{2})$ :

$$ \vec{u} = \begin{cases} u(x, \theta) - \rho w, & x \ne 0 \\ v(x, \theta) - \frac{\rho}{R}(w, \theta - v) & x \ne 0 \\ w(x, \theta) & x = 0 \end{cases} \qquad [32] $$

Les conditions aux limites de bords libres conduisent à chercher une solution fondamentale comportant un champ de déplacement longitudinal qui ne dépend que de la variable $x$, un déplacement radial constant sur l'ensemble du cylindre et un champ orthoradial identiquement nul :

$$ \vec{u}_f = \begin{cases} u_f(x) \\ 0 \\ w_f \end{cases} \qquad [33] $$

La relation entre le déplacement et la déformation donne :

$$ \begin{aligned} E_{xx} &= u_{f,x} \\ E_{\theta\theta} &= \frac{w_f}{R} \end{aligned} \qquad [34] $$

En égalant [31] et [34], on obtient la solution fondamentale :

$$ \vec{u}_f = \begin{cases} -\frac{\lambda}{E}x - \frac{\lambda - \sigma_0}{h}x & \\ 0 & \\ R\left(\frac{\nu\lambda}{E} + \frac{\lambda - \sigma_0}{2h}\right) & \end{cases} \qquad [35] $$

### 3.2.2. Equation de bifurcation

Le calcul de la charge critique et du mode de bifurcation passe par la résolution de l'équation de bifurcation [20]. On construit le mode de bifurcation à partir de trois champs $(\mathcal{U}, \mathcal{V}$ et $\mathcal{W})$, conformément à la cinématique [32] :

$$ \vec{X} = \begin{cases} \mathcal{U}(x, \theta) - \rho \mathcal{W}, & x \ne 0 \\ \mathcal{V}(x, \theta) - \frac{\rho}{R}(\mathcal{W}, \theta - \mathcal{V}) & x \ne 0 \\ \mathcal{W}(x, \theta) & x = 0 \end{cases} \qquad [36] $$

Le champ de déplacement virtuel adopté est de la même forme :

$$ \delta \vec{u} = \begin{cases} \delta u(x, \theta) - \rho \delta w, & x \ne 0 \\ \delta v(x, \theta) - \frac{\rho}{R}(\delta w, \theta - \delta v) & x \ne 0 \\ \delta w(x, \theta) & x = 0 \end{cases} \qquad [37] $$
Le gradient du mode de bifurcation, tout comme celui du déplacement virtuel, s'écrit :

$$ \bar{\nabla} \vec{X} = \begin{bmatrix} U,_{x} - \rho W,_{xx} & \frac{1}{R}(U,_{\theta} - \rho W,_{x\theta}) & -W,_{x} \\ V,_{x} - \frac{\rho}{R}(W,_{\theta x} - V,_{x}) & \frac{1}{R}(V,_{\theta} - \frac{\rho}{R}(W,_{\theta\theta} - V,_{\theta}) + W) & \frac{1}{R}(V - W,_{\theta}) \\ W,_{x} & \frac{1}{R}(W,_{\theta} - V + \frac{\rho}{R}(W,_{\theta} - V)) & 0 \end{bmatrix} \quad [38] $$

Dans le cas d'une structure mince de type coque, on considère que la contrainte normale $\Sigma_{rr}$ est nulle partout (alors que la cinématique choisie induit, de manière contradictoire, une déformation $E_{rr}$ nulle) et on remanie par conséquent les relations de comportement tridimensionnelles afin de prendre en compte cette condition de contraintes planes $\Sigma_{rr} = 0$. La loi de comportement bidimensionnelle qui en résulte utilise les tenseurs réduits suivants :

$$ \begin{aligned} C_{ijkl}^{e} &= D_{ijkl}^{e} - \frac{D_{33kl}^{e}}{D_{333}^{e}} D_{ij33}^{e}, \\ C_{ijkl}^{p} &= D_{ijkl}^{p} - \frac{D_{33kl}^{p}}{D_{333}^{p}} D_{ij33}^{p}, \quad (i, j) \neq (3, 3), (k, l) \neq (3, 3) \end{aligned} \quad [39] $$

Nous ne développerons pas l'écriture fastidieuse de l'équation de bifurcation [20]. Notons simplement qu'après intégration dans l'épaisseur, les intégrations par parties appropriées en $x$ et $\theta$ nous conduisent aux équations locales suivantes :

$$ (\alpha - \lambda_c)tU,_{xx} + \frac{\mu}{R^2}tU,_{\theta\theta} + \frac{\beta+\mu}{R}tV,_{x\theta} + \frac{\beta}{R}tW,_{x} = 0 $$

$$ \frac{\beta+\mu}{R}tU,_{x\theta} - \frac{\mu}{R^2}\frac{t^3}{12}V + (\mu - \lambda_c)tV,_{xx} + \frac{\mu-\lambda_c}{R^2}\frac{t^3}{12}V,_{xx} + \frac{\gamma}{R^2}tV,_{\theta\theta} + \frac{\gamma}{R^4}\frac{t^3}{12}V,_{\theta\theta} \\ + \frac{\mu}{R^2}\frac{t^3}{12}W,_{\theta} - \frac{\beta+2\mu-\lambda_c}{R^2}\frac{t^3}{12}W,_{xx\theta} - \frac{\gamma}{R^4}\frac{t^3}{12}W,_{\theta\theta\theta} = 0 $$

$$ \frac{\beta}{R}tU,_{x} + \frac{\gamma}{R^2}tV,_{\theta} + \frac{\mu}{R^4}\frac{t^3}{12}V,_{\theta} - \frac{\beta+2\mu-\lambda_c}{R^2}\frac{t^3}{12}V,_{xx\theta} - \frac{\gamma}{R^4}\frac{t^3}{12}V,_{\theta\theta\theta} + \frac{\gamma}{R^2}tW + \lambda_c tW,_{xx} \\ - \frac{\mu}{R^4}\frac{t^3}{12}W,_{\theta\theta} + (\alpha - \lambda_c)\frac{t^3}{12}W,_{xxxx} + \frac{2\beta+4\mu-\lambda_c}{R^2}\frac{t^3}{12}W,_{xx\theta\theta} + \frac{\gamma}{R^4}\frac{t^3}{12}W,_{\theta\theta\theta\theta} = 0 \quad [40] $$

où on a noté en abrégé les modules de [39] :

$$ \begin{aligned} \alpha &= C_{1111}^{p} = \frac{2\mu}{D}(E + 4h) \\ \beta &= C_{1122}^{p} = C_{2211}^{p} = \frac{4\mu}{D}(E + 2\nu h) \\ \gamma &= C_{2222}^{p} = \frac{8\mu}{D}(E + h) \end{aligned} \quad [41] $$

avec $D = 2\mu(5 - 4\nu) + 4h(1 - \nu)$.

Dans le cas axisymétrique où les déplacements orthoradiaux sont nuls et où toute grandeur est indépendante de la variable $\theta$, il ne reste plus que deux parmi les trois équations [40], qui s'écrivent après omission des termes négligeables :

$$ \begin{aligned} & \alpha t U,_{x} + \frac{\beta}{R} t W = 0 \\ & \alpha \frac{t^3}{12} W,_{xxxx} + \lambda_c t W,_{xx} + \frac{\gamma}{R^2} t W + \frac{\beta}{R} t U,_{x} = 0 \end{aligned} \quad [42] $$
La combinaison de ces deux équations mène à une équation portant sur la seule composante radiale du mode :

$$ \alpha \frac{t^3}{12} \mathcal{W}_{xxxx} + \lambda_c t \mathcal{W}_{xx} + \left( \gamma - \frac{\beta^2}{\alpha} \right) \frac{t}{R^2} \mathcal{W} = 0 \quad [43] $$

On retrouve une équation semblable à celle bien connue en élasticité.

### 3.2.3. Résolution dans le cas des bords libres

L'équation de bifurcation [43] s'applique directement au cas du cylindre à bords libres. La solution fondamentale [35] est très simple, par contre les calculs de bifurcation sont compliqués et on donne sans détailler la valeur de la contrainte critique :

$$ \lambda_c^l = \frac{Et}{R} \sqrt{\frac{h}{3[(5 - 4\nu)E + 4h(1 - \nu^2)]}} \quad [44] $$

où l'exposant $l$ rappelle qu'il s'agit du cas de bords libres. Le mode de bifurcation quant à lui n'est pas simple à obtenir, il faut passer au cas suivant des bords appuyés pour davantage de résultats analytiques.

### 3.2.4. Résolution dans le cas des bords appuyés

On considère ici un cylindre à bords appuyés : aux deux extrémités du cylindre les déplacements radiaux sont bloqués, seules les rotations y sont libres. On peut de nouveau utiliser l'équation de bifurcation [43] pour calculer la charge critique et le mode de bifurcation, moyennant l'approximation qui consiste à prendre comme solution fondamentale celle du cas des bords libres. Cette approximation, motivée par les faibles déformations pré-critiques, sera validée par les résultats numériques obtenus à la section 4.3.1.

#### 3.2.4.1. Charge critique et mode de bifurcation

De même qu'en élasticité, on est amené à supposer un mode de bifurcation sinusoïdal :

$$ \mathcal{W}^a = \sin \frac{n\pi x}{l} \quad (n \in \mathbb{N}) \quad [45] $$

de façon à respecter naturellement les conditions aux limites (l'exposant $a$ se rapporte aux bords appuyés). En remplaçant $\mathcal{W}^a$ dans [43], on obtient la valeur critique sous la forme :

$$ \lambda_c^{a(n)} = \alpha \frac{t^2 n^2 \pi^2}{12 l^2} + \left( \gamma - \frac{\beta^2}{\alpha} \right) \frac{1}{R^2} \frac{l^2}{n^2 \pi^2} \quad [46] $$

Le mode recherché s'obtient en minimisant $\lambda_c^{a(n)}$ avec $n \in \mathbb{N}$. Pour un cylindre suffisamment long, la minimisation sur $\mathbb{R}$ donne une bonne approximation de la valeur critique :

$$ \lambda_c^a = 2\lambda_c^l = \frac{Et}{R} \sqrt{\frac{4h}{3[(5 - 4\nu)E + 4h(1 - \nu^2)]}} \quad [47] $$
valeur qui correspond à :

$$n^a = \frac{l}{\pi} \sqrt[4]{\frac{48h[(5 - 4\nu)E + 4h(1 - \nu^2)]}{R^2 t^2 (E + 4h)^2}} \quad [48]$$

Le mode complet s'obtient en reprenant le système [42] :

$$\left\{
\begin{array}{l}
\mathcal{U}^a = \frac{l}{n^a \pi R} \frac{2(E+2\nu h)}{E+4h} \left(\cos \frac{n^a \pi x}{l} - 1\right) \\
\mathcal{V}^a = 0 \\
\mathcal{W}^a = \sin \frac{n^a \pi x}{l}
\end{array}
\right.
\quad [49]$$

### 3.2.4.2. Calcul de $\lambda_1^a$

Le paramètre $\lambda_1^a$ noté $\lambda_m^a$ dans le cas d'une bifurcation avec décharge naissante est donné par [21] :

$$\lambda_m^a = \max_{\Omega} \left( - \frac{\bar{M}_c : \bar{\nabla} \vec{X}^a}{\bar{M}_c : \bar{\nabla} \vec{u}_f,_{\lambda}(\lambda_c)} \right) \quad [50]$$

Ici $\Omega_f^p(t_c) = \Omega$, volume de référence du cylindre, et :

$$\bar{M}_c = \frac{\mu}{\sqrt{h + 3\mu}} \begin{bmatrix} -2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad [51]$$

$$\bar{\nabla} \vec{X}^a = \begin{bmatrix} \left(-\frac{2(E+2\nu h)}{R(E+4h)} + \rho \frac{n^2 \pi^2}{l^2}\right) \sin \frac{n\pi x}{l} & 0 & -\frac{n\pi}{l} \cos \frac{n\pi x}{l} \\ 0 & \frac{1}{R} \sin \frac{n\pi x}{l} & 0 \\ \frac{n\pi}{l} \cos \frac{n\pi x}{l} & 0 & 0 \end{bmatrix} \quad [52]$$

$$\bar{\nabla} \vec{u}_f,_{\lambda}(\lambda_c) = \begin{bmatrix} -\left(\frac{1}{E} + \frac{1}{h}\right) & 0 & 0 \\ 0 & \frac{\nu}{E} + \frac{1}{2h} & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad [53]$$

On obtient alors :

$$\lambda_1^a = \frac{2Eh[5E + 4h + 8\nu h + \sqrt{48h[(5 - 4\nu)E + 4h(1 - \nu^2)]}]}{R(E + 4h)(5E + 4h + 2\nu h)} \quad [54]$$

On remarque, au passage, que le cas des bords encastrés se traite de manière similaire à celui des bords appuyés.
# 4. Modélisation numérique du flambement élastoplastique d'un cylindre sous compression axiale

## 4.1. Cinématique, déformations et contraintes

Le modèle de coque élastoplastique en grandes rotations développé dans ce travail est issu de celui décrit dans [BRA 97] auquel nous avons apporté quelques modifications.

La position de la coque est définie par le champ de vecteurs position des points de la surface moyenne et un champ de vecteurs directeurs unitaires indiquant en chaque point la direction de la fibre suivant l'épaisseur. Le vecteur directeur $\vec{T}$ dans la configuration de référence est choisi normal à la surface moyenne, contrairement à celui dans la configuration déformée $\vec{t}$ qui ne le sera pas nécessairement, conformément à la formulation de Mindlin-Reissner.

Le déplacement d'un point $(\xi, \eta, \zeta)$, où $\xi$ et $\eta$ représentent les coordonnées sur la surface moyenne et $\zeta$ la position dans l'épaisseur, est noté $\vec{D}(\xi, \eta, \zeta) = \vec{U}(\xi, \eta) + \zeta(\vec{t}(\xi, \eta) - \vec{T}(\xi, \eta))$ avec $\vec{U}(\xi, \eta)$ le déplacement du point correspondant appartenant à la surface moyenne.

Le mouvement du voisinage d'un point est caractérisé par le tenseur gradient de la transformation $\vec{F}$ défini en tout point par $\vec{F} = \vec{a}_{\alpha} \otimes \vec{A}^{\alpha} + \vec{t} \otimes \vec{T} + \zeta \vec{t}_{,\alpha} \otimes \vec{A}^{\alpha}$ où $\{\vec{A}_i\}_{i=1,2,3}$ et $\{\vec{a}_i\}_{i=1,2,3}$ définissent les bases locales covariantes dans les configurations de référence et actuelle et $\{\vec{A}^i\}_{i=1,2,3}$ est la base contravariante duale de $\{\vec{A}_i\}_{i=1,2,3}$, relativement à la surface moyenne, négligeant ainsi la variation de métrique dans l'épaisseur.

L'évolution du vecteur directeur est donnée par $\vec{t} = \vec{R} \cdot \vec{T}$ où le tenseur $\vec{R} = \exp \vec{\Theta} = \sum_{k=0}^{+\infty} \frac{\vec{\Theta}^k}{k!}$ est orthogonal ($\vec{R}^{-1} = \vec{R}^T$) et définit la rotation finie autour de $\vec{\theta}$, vecteur adjoint du tenseur antisymétrique $\vec{\Theta}$. La formule de Rodrigues fournit une expression simple de $\vec{R}$ :

$$ \vec{R} = \cos \theta \vec{I} + \frac{\sin \theta \vec{\Theta}}{\theta} + \frac{1 - \cos \theta}{\theta^2} \vec{\theta} \otimes \vec{\theta} \quad [55] $$

où $\theta = \|\vec{\theta}\|$ désigne l'amplitude de la rotation.

Le tenseur des déformations de Green-Lagrange $\vec{\bar{E}}$ s'écrit :

$$ \vec{\bar{E}} = \frac{1}{2}(\vec{F}^T \cdot \vec{F} - \vec{I}) = \frac{1}{2}(\vec{a}_i \cdot \vec{a}_j - \vec{A}_i \cdot \vec{A}_j) \vec{A}^i \otimes \vec{A}^j \quad [56] $$
Si on décompose $E_{ij} = E_{ij}^{(0)} + \zeta E_{ij}^{(1)} + \zeta^2 E_{ij}^{(2)}$ et que l'on néglige le terme du second ordre en $\zeta$, il reste, compte tenu de $\vec{a}_{\alpha} = \vec{A}_{\alpha} + \vec{U}_{,\alpha}$ :

$$
\begin{align*}
E_{\alpha\beta}^{(0)} &= \frac{1}{2}(\vec{A}_{\alpha} \cdot \vec{U}_{,\beta} + \vec{A}_{\beta} \cdot \vec{U}_{,\alpha} + \vec{U}_{,\alpha} \cdot \vec{U}_{,\beta}) \\
E_{\alpha 3}^{(0)} &= E_{3\alpha}^{(0)} = \frac{1}{2}(\vec{A}_{\alpha} \cdot (\vec{t} - \vec{T}) + \vec{U}_{,\alpha} \cdot \vec{t}) \tag{57} \\
E_{\alpha\beta}^{(1)} &= \frac{1}{2}(\vec{A}_{\alpha} \cdot (\vec{t}_{,\beta} - \vec{T}_{,\beta}) + \vec{A}_{\beta} \cdot (\vec{t}_{,\alpha} - \vec{T}_{,\alpha}) + \vec{U}_{,\alpha} \cdot \vec{t}_{,\beta} + \vec{U}_{,\beta} \cdot \vec{t}_{,\alpha})
\end{align*}
$$

$E_{33}^{(0)}$ et $E_{33}^{(1)}$ sont nuls du fait de l'inextensibilité du vecteur directeur.

L'utilisation du principe des travaux virtuels nécessite de connaître les variations de ces composantes :

$$
\begin{align*}
\delta E_{\alpha\beta}^{(0)} &= \frac{1}{2}((\vec{A}_{\alpha} + \vec{U}_{,\alpha}) \cdot \delta \vec{U}_{,\beta} + (\vec{A}_{\beta} + \vec{U}_{,\beta}) \cdot \delta \vec{U}_{,\alpha}) \\
\delta E_{\alpha 3}^{(0)} &= \delta E_{3\alpha}^{(0)} = \frac{1}{2}((\vec{A}_{\alpha} + \vec{U}_{,\alpha}) \cdot \delta \vec{t} + \delta \vec{U}_{,\alpha} \cdot \vec{t}) \tag{58} \\
\delta E_{\alpha\beta}^{(1)} &= \frac{1}{2}((\vec{A}_{\alpha} + \vec{U}_{,\alpha}) \cdot \delta \vec{t}_{,\beta} + (\vec{A}_{\beta} + \vec{U}_{,\beta}) \cdot \delta \vec{t}_{,\alpha} + \delta \vec{U}_{,\alpha} \cdot \vec{t}_{,\beta} + \delta \vec{U}_{,\beta} \cdot \vec{t}_{,\alpha})
\end{align*}
$$

Le second tenseur des contraintes de Piola-Kirchhoff de [10] est décomposé dans
la base covariante :

$$
\bar{\Sigma} = \Sigma^{\alpha\beta} \vec{A}_{\alpha} \otimes \vec{A}_{\beta} + \Sigma^{\alpha 3} \vec{A}_{\alpha} \otimes \vec{T} + \Sigma^{3\alpha} \vec{T} \otimes \vec{A}_{\alpha} + \Sigma^{33} \vec{T} \otimes \vec{T} \quad [59]
$$

4.2. Résolution par éléments finis

La discrétisation du principe des travaux virtuels conduit à un système matriciel
non linéaire que l’on résout de manière classique par le schéma itératif de Newton-
Raphson. A chaque itération, l’intégration locale fournit notamment la matrice des
contraintes Σ et l’opérateur tangent consistant ∂Σ/∂E en fonction de l’incrément de
déformation, en contraintes planes. La discrétisation en temps des relations en vi-
tesses se fait avec le schéma d’Euler implicite. L’évaluation du vecteur force interne
ainsi que de la matrice raideur nécessite des intégrations dans le volume, qui sont ef-
fectuées par un procédé d’intégration numérique dans l’épaisseur. La programmation
a été effectuée avec le critère de von Mises pour un écrouissage mixte (isotrope non
linéaire et cinématique linéaire). La correction plastique décrite dans [SIM 86] permet
de prendre en compte correctement la condition de contraintes planes.

On évalue la fonction seuil $f$ connaissant les prédictions élastiques $\eta^E$ et $p^E$ ; si $f$
est négative, $C_{n+1}^p = \frac{\partial \Sigma_{n+1}}{\partial E_{n+1}} = C^e$, sinon, on évalue le multiplicateur plastique $\lambda$ en
résolvant :

$$
f = \frac{1}{2} \eta_{n+1}^{T} P \eta_{n+1} - \frac{1}{3} (\sigma_{0} + K(p_{n+1}))^{2} = 0 \quad [60]
$$

avec :

$$
\begin{align*}
\eta_{n+1} &= \frac{1}{1+\frac{2}{3}\lambda H} \Gamma(\lambda) C^{-1} \eta^E , & \Gamma(\lambda) &= (C^{-1} + \frac{\lambda}{1+\frac{2}{3}\lambda H} P)^{-1} \tag{61} \\
p_{n+1} &= p_n + \lambda \sqrt{\frac{2}{3} \eta_{n+1}^T P \eta_{n+1}}
\end{align*}
$$
et on met à jour les autres quantités, dont la matrice tangente consistante :

$$C_{n+1}^{p} = \frac{\partial \Sigma_{n+1}}{\partial E_{n+1}} = \Gamma(\lambda) - \frac{(\Gamma(\lambda)P\eta_{n+1})(\Gamma(\lambda)P\eta_{n+1})^T}{\eta_{n+1}^T P^T \Gamma(\lambda) P\eta_{n+1} + \mu} \quad [62]$$

avec :

$$\begin{aligned} \mu &= \frac{2}{3} \frac{\gamma_1}{\gamma_2} (K'(p_{n+1})\gamma_1 + H\gamma_2)\eta_{n+1}^T P\eta_{n+1} \\ \gamma_1 &= 1 + \frac{2}{3}\lambda H, \quad \gamma_2 = 1 - \frac{2}{3}\lambda K'(p_{n+1}) \end{aligned} \quad [63]$$

Afin de pouvoir franchir les points limites (extrema de force ou de déplacement) qui sont relativement fréquents dans les problèmes d'instabilité, nous avons mis en place une méthode de longueur d'arc quadratique généralisant l'approche de [LAM 92] en permettant de gérer à la fois les forces imposées et les déplacements imposés.

Notre méthode s'appuie sur la méthode quadratique classique, développée par Crisfield [CRI 83],[CRI 84], consistant à introduire, à l'itération $i + 1$, la condition :

$$\Delta U_{i+1}^T \Delta U_{i+1} + \Delta \lambda_{i+1}^2 \psi^2 F_r^T F_r = \Delta l^2 \quad [64]$$

où $U$ est le vecteur des déplacements nodaux, $F_r$ le vecteur de référence des forces nodales équivalentes, $\lambda$ le facteur de charge, $\psi$ un facteur d'échelle, $l$ la longueur d'arc et $\Delta$ se rapporte aux grandeurs incrémentales.

Etant donnée la correction en charge $\delta\lambda_i$ effectuée à l'itération $i + 1$, la correction similaire à apporter en déplacement se décompose en deux parties :

$$\delta U_i = -K_T^{-1} R(U_i, \lambda_{i+1}) = -K_T^{-1} (R(U_i, \lambda_i) - \delta\lambda_i F_r) = \delta U_{iR} + \delta\lambda_i \delta U_{iF} \quad [65]$$

avec $K_T$ la matrice raideur tangente et $R$ le résidu à l'itération courante.

La relation [64] conduit à une équation du second degré en $\delta\lambda_i$ qui a parfois des racines complexes. Lam et Morley [LAM 92] proposent alors de projeter le résidu sur le vecteur force externe en décomposant, à l'itération $i$, le résidu $R_i$ en sa projection $G_i = \frac{R_i^T F_r}{F_r^T F_r}$ et une partie orthogonale à $F_r$ qu'on notera $H_i = R_i - G_i F_r$. Le problème des racines complexes est ensuite résolu en réduisant au mieux la contribution du terme orthogonal $H_i$ responsable de leur présence.

La correction itérative [65] est tout d'abord remplacée par :

$$\delta U_i = \delta U_{iH} + x_i \delta U_{iF} \quad [66]$$

avec $\delta U_{iH} = -K_T^{-1} H_i$, ce qui nous permet d'obtenir une nouvelle équation du second degré en $x_i$. En cas de racines complexes, il est alors possible de pondérer, durant la correction, le terme responsable relatif au résidu orthogonal.

Nous envisageons ensuite de généraliser ce pilotage par longueur d'arc au cas des déplacements imposés. Dans cette configuration, on définit différemment la notion de longueur d'arc en faisant intervenir le vecteur déplacement imposé de référence $U_r$ à la place du vecteur $F_r$ dans [64]. Ce passage des forces imposées aux déplacements
imposés modifie également la forme générale des équations d'équilibre et l'écriture de la correction itérative en déplacement [65]. La racine à chaque itération, le paramètre d'échelle et la prédiction à tous les incréments sont choisis de façon appropriée, ce qui permet de traiter de manière unifiée les deux cas de forces et déplacements imposés.

Par ailleurs, des techniques spécifiques ont été développées à l'instar de Riks [RIK 79],[RIK 91] et Seydel [SEY 94], de manière à pouvoir choisir une branche arbitraire de bifurcation sur une structure sans imperfection.

Nous retenons globalement une méthode dite indirecte (peu coûteuse) pour le traitement des points de bifurcation rencontrés durant la démarche de calcul incrémentale. On distingue alors quatre étapes fondamentales à notre approche :

- on commence par détecter, à chaque nouvel incrément, un éventuel changement de nature de la matrice raideur tangente, caractéristique du passage d'un point critique, à l'aide d'une « fonction test » telle que le déterminant de la matrice ou ses pivots issus d'une factorisation de Crout ;

- on utilise ensuite un paramètre de « raideur scalaire équivalente », traduisant en quelque sorte la pente de la courbe force-déplacement associée, afin de préciser la nature du point critique détecté (point limite ou point de bifurcation) ;

- s'il s'agit d'un point de bifurcation, on le situe au mieux par interpolation sur les incréments l'encadrant ;

- après la détermination d'un point de bifurcation, on peut choisir de continuer d'évoluer sur la branche actuelle ou bien de bifurquer. Le cas échéant, la description de la branche bifurquée n'est possible qu'en guidant volontairement la solution vers un nouveau point de cette branche, moyennant le choix d'une prédiction et d'une correction adéquates. A l'instar de Seydel [SEY 94], on procède alors, au point de bifurcation, à une prédiction basée sur une approximation modale en déplacement et à une correction classique.

## 4.3. Résultats numériques

Les résultats numériques obtenus s'orientent vers deux objectifs. Les premiers calculs, relatifs à des solutions axisymétriques, serviront de comparaison avec les résultats analytiques développés au paragraphe 3. Ils permettront ainsi une validation de notre programme quant à sa possibilité de modélisation d'un tel comportement critique de cylindre sous compression axiale. Dans un deuxième temps, l'approche numérique va permettre d'aller plus loin, en obtenant des solutions non axisymétriques, faisant apparaître, entre autres, des modes « diamant » et des modes secondaires « patte d'éléphant ». Ces derniers résultats montrent la capacité de notre programme à obtenir la quasi-totalité des modes rencontrés expérimentalement et permettent de dégager certaines grandes règles quant à l'influence des différents paramètres géométriques et matériels sur le type de mode obtenu.
On présente ici uniquement les solutions les plus représentatives, sélectionnées parmi les nombreux calculs effectués. Notons, au passage, que l'absence de prise en compte des phénomènes de contact dans notre programme empêche le prolongement de certains calculs jusqu'à des situations post-critiques encore plus lointaines.

Après quelques tests préalables de convergence vis-à-vis de la finesse du maillage, on décide de se limiter dans tous les calculs qui suivent à des maillages raisonnables de l'ordre de 200 éléments (quadrangles à 8 nœuds isoparamétriques à 5 degrés de liberté par nœud) pour la portion du cylindre modélisée et à 4 points de Gauss pour l'intégration numérique dans l'épaisseur.

### 4.3.1. Comparaison théorique-numérique

On calcule ici le cylindre du sous-paragraphe 3.2 avec les valeurs numériques des tableaux 1 et 2.

<table><caption>**Tableau 1.** Géométrie utilisée à la section 4.3.1</caption><thead><tr><td>Longueur l</td><td>Rayon R</td><td>Epaisseur t</td></tr></thead><tbody><tr><td>2 m</td><td>1 m</td><td>1 cm</td></tr></tbody></table>

<table><caption>**Tableau 2.** Propriétés matérielles utilisées à la section 4.3.1</caption><thead><tr><td>Module d'Young E</td><td>Coefficient de Poisson ν</td><td>Limite d'élasticité σ<sub>0</sub></td><td>Module d'écrouissage isotrope h</td></tr></thead><tbody><tr><td>2.1 10<sup>11</sup> Pa</td><td>0.3</td><td>Bords libres : 1.10<sup>8</sup> Pa<br>Bords appuyés : 3.5 10<sup>8</sup> Pa</td><td>2.10<sup>10</sup> Pa</td></tr></tbody></table>

#### 4.3.1.1. Bords libres

Dans le cas des bords libres, on choisit la contrainte limite élastique $\sigma_0$ suffisamment petite pour bifurquer dans le domaine plastique. Le calcul donne une bifurcation axisymétrique à charge croissante à la valeur critique de 184 MPa qui est égale (à moins de 1 % près) à la contrainte analytique $\lambda_c^l$ obtenue grâce à [44]. Le mode de bifurcation a la forme d'une sinusoïde amortie au centre du cylindre et d'amplitude maximale aux extrémités, conformément aux résultats de Batterman [BAT 68]. La figure 1 montre que la déformée modale numérique est très proche de celle analytique d'un cylindre semi-infini de Batterman.

#### 4.3.1.2. Bords appuyés

On trouve dans le cas des bords simplement appuyés une valeur critique double de la précédente, en accord avec la solution analytique $\lambda_c^a$ obtenue en [47]. De manière qualitative, la bifurcation précède tout juste le maximum de la charge, comme cela a été fréquemment observé dans d'autres applications numériques.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Cylindre élastoplastique à bords libres sous compression axiale : déformées modales axisymétriques analytique et numérique</figcaption></figure>

La figure 2 représente la courbe charge-déplacement radial maximal dans ce cas particulier. La relation [54] nous permet de tracer la tangente théorique au point de bifurcation. La figure 2 montre que cette tangente se raccorde très bien à la courbe numérique.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cylindre élastoplastique à bords appuyés sous compression axiale : évolution de la charge axiale en fonction du déplacement radial maximal</figcaption></figure>
### 4.3.2. Comparaison avec un résultat numérique de la littérature

En restant avec les modes axisymétriques, nous allons comparer nos résultats numériques à ceux de [BRA 97], basés également sur une théorie incrémentale de la plasticité avec un écrouissage isotrope non linéaire. Nous prenons donc la même géométrie et le même matériau que dans [BRA 97] (tableaux 3 et 4).

<table><caption>**Tableau 3.** Géométrie utilisée à la section 4.3.2</caption><thead><tr><td>Longueur l</td><td>Rayon R</td><td>Epaisseur t</td></tr></thead><tbody><tr><td>1.271 10<sup>-1</sup> m</td><td>4.333 10<sup>-2</sup> m</td><td>2.36 10<sup>-3</sup> m</td></tr></tbody></table>

<table><caption>**Tableau 4.** Propriétés matérielles utilisées à la section 4.3.2</caption><thead><tr><td>Module d'Young E</td><td>Coefficient de Poisson ν</td><td>Limite d'élasticité σ<sub>0</sub></td></tr></thead><tbody><tr><td>2.16962 10<sup>11</sup> Pa</td><td>0.274</td><td>1.62722 10<sup>8</sup> Pa</td></tr></tbody></table>

Quant aux coefficients de la loi d'écrouissage isotrope non linéaire, ils sont déterminés par lissage de la courbe donnée dans [BRA 97]. Les calculs sont effectués avec deux conditions aux limites : bords simplement appuyés et bords encastrés. Les courbes d'évolution de l'effort axial en fonction du raccourcissement axial (figure 3) présentent toutes deux un point de bifurcation au voisinage du point limite et sont en très bonne concordance avec les solutions de Brank, Peric et Damjanic. De plus, elles sont très proches l'une de l'autre, ce qui corrobore la remarque faite aux sections 2.1.2 et 3.2.4, selon laquelle les deux types de conditions aux limites conduisent à des résultats très similaires (figure 4). La figure 5 illustre la répartition de la déformation plastique équivalente sur la surface du cylindre (respectivement en moyenne et au maximum dans l'épaisseur de la coque). Cette déformation plastique est maximale au niveau de la localisation axisymétrique aux extrémités du cylindre (zone la plus sombre) et atteint la valeur de 5.6 % à 1.6 % de raccourcissement axial pour le cas des appuis simples et de 6.7 % à 2.3 % de raccourcissement axial pour les bords encastrés.

### 4.3.3. Influence de l'épaisseur, du module d'écrouissage et des conditions aux limites

Après les validations théorique/numérique et numérique/numérique aux sections 4.3.1 et 4.3.2, on utilise dans la suite le programme de calcul pour explorer les modes plus évolués. On fera varier les paramètres géométriques et matériels, ainsi que les conditions aux limites (bords libres, appuyés ou encastrés, contraintes ou déplacements imposés), et on étudiera la réponse post-critique du cylindre sous compression axiale. Il est prévisible que les conditions aux limites aient une forte influence sur l'al-lure de la géométrie déformée, en particulier si elle est non axisymétrique. En effet, on s'aperçoit par exemple dans le cas de bords libres, qu'un déplacement longitudinal imposé aux extrémités conduit plus naturellement à un mode axisymétrique, alors
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Cylindre élastoplastique sous compression axiale : courbes contrainte de compression-raccourcissement axial</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cylindre élastoplastique sous compression axiale : déformées axisymétriques</figcaption></figure>

qu'une force axiale imposée autorise davantage, par la liberté de mouvement de ces extrémités, un mode non axisymétrique.

La gamme des résultats obtenus dans cette section 4.3.3 utilise la géométrie et les caractéristiques mécaniques des tableaux 5 et 6.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cylindre élastoplastique sous compression axiale : déformation plastique équivalente</figcaption></figure>

<table><caption>**Tableau 5.** Géométrie utilisée à la section 4.3.3</caption><thead><tr><td>Longueur l</td><td>Rayon R</td><td>Epaisseur t</td></tr></thead><tbody><tr><td>5 m</td><td>1.5 m</td><td>7 cm, 5 cm, 3 cm, 2 cm, 12 mm, 6 mm selon le cas</td></tr></tbody></table>

Dans l'ensemble des courbes tracées par la suite, on ramène la contrainte de compression $\lambda$ à une grandeur adimensionnelle en la normalisant par la contrainte limite élastique $\sigma_0$, de manière à analyser plus aisément les différents calculs effectués.

### 4.3.3.1. Bords libres

En choisissant une épaisseur relativement importante ($t = 7$ cm), un écrouissage intermédiaire ($h = 2.10^{10}$ Pa) et en forces imposées, on obtient la première solution bifurquée comportant trois ondes circonférencielles représentées sur la figure 6. Cette déformée est conforme à l'observation expérimentale ainsi qu'aux résultats numériques obtenus par Goto et Zhang [GOT 99]. La figure 7 décrit la distribution de la déformation plastique sur l'ensemble du cylindre. Un maximum de 8.4 % est atteint à la jonction de deux ondes circonférencielles successives toujours localisées aux extrémités du cylindre.

En conservant les mêmes conditions aux limites et le même chargement, on diminue l'épaisseur du cylindre à $t = 5$ cm, de façon à observer l'influence du rapport $R/t$ sur le nombre d'ondes circonférencielles. La figure 8 montre qu'il y a bien une aug-
<table><caption>**Tableau 6.** Propriétés matérielles utilisées à la section 4.3.3</caption><thead><tr><th>Module d'Young E</th><th>Coefficient de Poisson ν</th><th>Limite d'élasticité σ<sub>0</sub></th><th>Module d'écrouissage isotrope h</th></tr></thead><tbody><tr><td>2.1 10<sup>11</sup> Pa</td><td>0.3</td><td>4.10<sup>8</sup> Pa</td><td>2.10<sup>10</sup> Pa, 5.10<sup>10</sup> Pa, 5.10<sup>8</sup> Pa selon le cas</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cylindre élastoplastique sous compression axiale : bords libres, forces imposées, $t = 7$ cm</figcaption></figure>

mentation du nombre d'ondes circonférencielles (on est passé de trois ondes à quatre) comme observé dans [GOT 99] et [TVE 83]. Le mode de bifurcation est localisé aux deux extrémités et comporte une série de « décrochements » périodiques vers l'intérieur du cylindre.

La figure 9 représente les courbes de réponse pour $t = 5$ cm et $t = 7$ cm. Dans les deux cas, la plastification (totale) du cylindre précède (tout juste pour l'épaisseur la plus faible) le point de bifurcation qui, lui-même, précède le maximum de la contrainte.

**REMARQUE.** — Concernant le maillage : signalons que pour exhiber trois ondes circonférencielles dans la figure 6, on a maillé un quart du cylindre. Afin de réduire le coût de calcul, on maillera désormais seulement un huitième du cylindre, sachant qu'ainsi on n'obtiendra que des déformées comportant un nombre pair d'ondes circonférencielles, comme c'est le cas de la figure 8. Concernant le type d'écrouissage : à titre de vérification, le dernier calcul a été effectué également pour un écrouissage non plus isotrope mais cinématique linéaire (avec $H = 2.10^{10}$ Pa), et les résultats obtenus sont globalement les mêmes.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Cylindre élastoplastique sous compression axiale : bords libres, forces imposées, $t = 7$ cm, déformation plastique équivalente</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Cylindre élastoplastique sous compression axiale : bords libres, forces imposées, $t = 5$ cm</figcaption></figure>

#### 4.3.3.2. Bords encastrés

On s'intéressera dans la suite aux conditions d'appui simple et d'encastrement, qui sont plus proches des conditions expérimentales. On choisit de prendre un module d'écrouissage isotrope linéaire $h = 5.10^{10}$ Pa, valeur assez forte afin d'obtenir des modes non axisymétriques étendus sur la surface du cylindre (modes « diamant ») ou limités à une région proche des extrémités du cylindre.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Cylindre élastoplastique sous compression axiale (bords libres, forces imposées) : courbes contrainte de compression-raccourcissement axial (○,● : début de plastification ; △,▲ : bifurcation)</figcaption></figure>

En forces imposées, on obtient deux modes « diamant » avec une épaisseur $t = 5 \ cm$ (puis $t = 3 \ cm$), présentés en figure 10, dont les nombres d'ondes aussi bien axiaux que circonférenciels augmentent avec le rapport $R/t$ (même constatation donc que dans le cas des bords libres).

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Cylindre élastoplastique sous compression axiale : bords encastrés, forces imposées, (a) $t = 5 \ cm$ - (b) $t = 3 \ cm$</figcaption></figure>
Les courbes de réponse (figure 11) montrent que le point de bifurcation apparaît à nouveau un peu avant le point limite. On remarque ici que le rapport entre les deux contraintes aux points de bifurcation est très proche du rapport des épaisseurs, ce qui rappelle la proportionnalité entre la contrainte analytique de la relation [47] et l'épaisseur du cylindre.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Cylindre élastoplastique sous compression axiale (bords encastrés, forces imposées) : courbes contrainte de compression-raccourcissement axial (○, ● : début de plastification ; △, ▲ : bifurcation)</figcaption></figure>

**REMARQUE.** — Les déplacements imposés donnent à peu près les mêmes résultats.

Les modes secondaires ne s'obtiennent qu'avec un module d'écrouissage beaucoup plus faible et qu'en déplacements imposés, ce qui favorise la localisation d'un mode axisymétrique. La figure 12 montre deux modes secondaires de type « patte d'éléphant » obtenus en déplacements imposés avec $h = 5.10^8 Pa$. On constate de nouveau que l'ordre circonférenciel dépend essentiellement de l'épaisseur du cylindre.

#### 4.3.3.3. Bords appuyés

Pour la condition d'appuis simples, en forces imposées et avec l'épaisseur $t = 2 cm$, le mode non axisymétrique obtenu ne s'étend pas sur l'ensemble du cylindre (figure 13(a)).

La figure 13(b) montre que dans le cas où l'on impose les déplacements aux extrémités, la déformée est plus étendue que dans le cas des forces imposées et le nombre de « fossettes » est plus important. On peut interpréter ceci par le fait qu'en contraintes imposées les cercles de base ne sont pas astreints à rester plans, la structure est donc
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Cylindre élastoplastique sous compression axiale : bords encastrés, déplacements imposés, (a) $t = 12$ mm - (b) $t = 6$ mm, modes secondaires</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Cylindre élastoplastique sous compression axiale : bords appuyés, (a) forces imposées - (b) déplacements imposés, $t = 2$ cm</figcaption></figure>

moins « rigidifiée » qu’en déplacements imposés et les fossettes peuvent s’estomper plus facilement.
Les deux courbes de réponse correspondantes (figure 14) sont très similaires jus-
qu'au point critique. Contrairement au cas des forces imposées, dans le cas des dépla-
cements imposés, la courbe post-critique présente plusieurs points limites successifs
qui doivent apparaître avec la formation des fossettes supplémentaires.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Cylindre élastoplastique sous compression axiale (bords appuyés, $t = 2 \text{ cm}$) : courbes contrainte de compression-raccourcissement axial (○,● : début de plastification ; △,▲ : bifurcation ; FI : Forces imposées ; DI : Déplacements imposés)</figcaption></figure>

**5. Conclusion**

Un modèle de coque en grandes rotations a été mis en œuvre avec une loi de comportement élastoplastique munie d'un écrouissage mixte. Devant la diversité des formulations en matière de grandes transformations élastoplastiques, nous avons opté pour une formulation lagrangienne totale où la déformation est décomposée de manière additive et les équations constitutives écrites sur le modèle des matériaux standard généralisés dont le cadre thermodynamique a fait ses preuves en plasticité petites transformations.

La formulation élément fini correspondante est associée à des techniques de longueur d'arc et de bifurcation suffisamment robustes et générales pour pouvoir résoudre les problèmes d'instabilité rencontrés. L'ensemble de ces développements numériques a permis d'obtenir un éventail de modes de flambement dans le cas du cylindre sous compression axiale, parmi lesquels des modes secondaires. Une classification de ces différents modes est présentée dans le tableau 7.

Le modèle présenté dans cette étude parvient à décrire l'ensemble de ces modes
dont l'apparition dépend principalement du module d'écrouissage, de l'épaisseur et
des conditions aux limites. En particulier, l'influence majeure de l'épaisseur sur les
<table><caption>**Tableau 7.** Classification des modes de bifurcation du cylindre sous compression axiale selon l'épaisseur et les conditions aux limites</caption><thead><tr><th>Cylindres</th><th>Simplement appuyés</th><th>Encastrés</th><th>Libres</th></tr></thead><tbody><tr><th>Epais</th><td>Axisymétrique<br>(+ forte localisation)<br>(figure 4)</td><td>Axisymétrique<br>(+ localisation)<br>(figure 4)</td><td>Localisation<br>axisymétrique<br>(figure 1)</td></tr><tr><th>Modérément<br>épais</th><td>-</td><td>Mode secondaire<br>« patte d'éléphant »<br>(figure 12)</td><td>-</td></tr><tr><th>Minces</th><td>Diamant<br>(figure 13)</td><td>Diamant<br>(figure 10)</td><td>Localisation non<br>axisymétrique<br>(figures 6 et 8)</td></tr></tbody></table>

nombres d'ondes axiales et circonférencielles a été clairement mis en valeur sur l'ensemble des types de modes rencontrés.

Les résultats de la section 4.3.3 permettent, sans avoir recours aux imperfections géométriques, d'exhiber la majeure partie des modes obtenus par les auteurs préalablement cités [TVE 83],[GOT 99], de manière automatique, en décrivant dans chaque cas le comportement post-critique jusqu'à un état relativement avancé des déformations. Les formes et les courbes de réponse obtenues sont en parfait accord avec la plupart des résultats expérimentaux rapportés dans [AFL 88] et [LIM 91].

Nous avons aussi établi la liaison entre les résultats numériques précédents et la théorie de bifurcation plastique [HUT 74],[NGU 00] et montré la très bonne correspondance entre les résultats numériques et les solutions analytiques envisagés - voir relations [44] et [47] et figures 1 et 2.

## 6. Bibliographie

[AFL 88] AFLAK W. F., « Flambage plastique de coques cylindriques sous compression axiale. Influence des imperfections géométriques et des imperfections de conditions aux limites », Thèse de doctorat, INSA de Lyon, 1988.

[AMA 72] AMAZIGO J. C., BUDIANSKY B., « Asymptotic formulas for the buckling stresses of axially compressed cylinders with localized or random axisymmetric imperfections », J. Applied Mechanics, 1972, p. 179-184.

[ARI 69] ARIARATNAM S. T., DUBEY R. N., « Instability in an elastic-plastic cylindrical shell under axial compression », J. Applied Mechanics, 1969, p. 47-50.

[BAT 65] BATTERMAN S. C., « Plastic buckling of axially compressed cylindrical shells », AIAA Journal, vol. 3, n° 2, 1965, p. 316-325.
[BAT 66] BATTERMAN S. C., LEE L. H. N., « Effect of modes on plastic buckling of compressed cylindrical shells », *AIAA Journal*, vol. 4, n° 12, 1966, p. 2255-2257.

[BAT 68] BATTERMAN S. C., « Free-edge plastic buckling of axially compressed cylindrical shells », *J. Applied Mechanics*, 1968, p. 73-79.

[BES 94] BESSELING J. F., VAN DER GIESSEN E., *Mathematical modelling of inelastic deformation*, Chapman & Hall, 1994.

[BRA 97] BRANK B., PERIC D., DAMJANIC F. B., « On large deformations of thin elasto-plastic shells : implementation of a finite rotation model for quadrilateral shell element », *Int. J. Numerical Methods in Engineering*, vol. 40, 1997, p. 689-726.

[BUD 72] BUDIANSKY B., HUTCHINSON J. W., « Buckling of circular cylindrical shells under axial compression », *Contributions to the Theory of Aircraft Structures*, Delft University Press, 1972, p. 239-259.

[CAS 80] CASEY J., NAGHDI P. M., « A remark on the use of the decomposition $F=FeFp$ in plasticity », *J. Applied Mechanics*, vol. 47, 1980, p. 672-675.

[CRI 83] CRISFIELD M. A., « An arc-length method including line searches and accelerations », *Int. J. Numerical Methods in Engineering*, vol. 19, 1983, p. 1269-1289.

[CRI 84] CRISFIELD M. A., « Accelerating and damping the modified Newton-Raphson method », *Computers and Structures*, vol. 18, n° 3, 1984, p. 395-407.

[DON 76] DONNELL L. H., *Beams, plates and shells*, Mc Graw-Hill, 1976.

[GIE 89a] VAN DER GIESSEN E., « Continuum models of large deformation plasticity, Part I : Large deformation plasticity and the concept of a natural reference state », *Eur. J. Mech. A/Solids*, vol. 8, 1989, p. 15-34.

[GIE 89b] VAN DER GIESSEN E., « Continuum models of large deformation plasticity, Part II : A kinematic hardening model and the concept of a plastically induced orientational structure », *Eur. J. Mech. A/Solids*, vol. 8, 1989, p. 89-108.

[GOT 99] GOTO Y., ZHANG C., « Plastic buckling transition modes in moderately thick cylindrical shells », *J. Engineering Mechanics*, vol. 125, n° 4, 1999, p. 426-434.

[GRE 65] GREEN A. E., NAGHDI P. M., « A general theory of an elastic-plastic continuum », *Archive for Rational Mechanics and Analysis*, vol. 18, 1965, p. 251-281.

[GRE 71] GREEN A. E., NAGHDI P. M., « Some remarks on elastic-plastic deformation at finite strain », *Int. J. Engineering Sciences*, vol. 9, 1971, p. 1219-1229.

[GUP 00] GUPTA N. K., ABBAS H., « Mathematical modeling of axial crushing of cylindrical tubes », *Thin-Walled Structures*, vol. 38, 2000, p. 355-375.

[GUS 00] GUSIC G., COMBESCURE A., JULLIEN J. F., « The influence of circumferential thickness variations on the buckling of cylindrical shells under external pressure », *Computers and Structures*, vol. 74, 2000, p. 461-477.

[HAL 75] HALPHEN B., NGUYEN Q. S., « Sur les matériaux standard généralisés », *Journal de Mécanique*, vol. 14, n° 1, 1975, p. 39-63.

[HUT 72] HUTCHINSON J. W., « On the post-buckling behaviour of imperfection-sensitive structures in the plastic range », *J. Applied Mechanics*, 1972, p. 155-162.

[HUT 74] HUTCHINSON J. W., « Plastic buckling », *Advances in Applied Mechanics*, vol. 14, 1974, p. 67-144.

[LAM 92] LAM W. F., MORLEY C. T., « Arc-length method for passing limit points in structural calculation », *J. Structural Engineering*, vol. 118, n° 1, 1992, p. 169-185.
[LAN 00] LANCASTER E. R., CALLADINE C. R., PALMER S. C., « Paradoxical buckling behaviour of a thin cylindrical shell under axial compression », *Int. J. Mechanical Sciences*, vol. 42, 2000, p. 843-865.

[LEE 67] LEE E. H., LIU D. T., « Finite-strain elastic-plastic theory with application to plane-wave analysis », *J. Applied Physics*, vol. 38, 1967, p. 19-27.

[LEE 69] LEE E. H., « Elastic-plastic deformations at finite strains », *J. Applied Mechanics*, vol. 36, 1969, p. 1-6.

[LIM 91] LIMAM A., « Flambage de coques cylindriques sous combinaison de chargements. Pression interne et compression axiale », *Thèse de doctorat*, INSA de Lyon, 1991.

[MAN 00] MANDAL P., CALLADINE C. R., « Buckling of thin cylindrical shells under axial compression », *Int. J. Solids Structures*, vol. 37, 2000, p. 4509-4525.

[NAG 75] NAGHDI P. M., TRAPP J. A., « Restrictions on constitutive equations for finitely deformed elastic-plastic materials », *Quart. J. Mech. Appl. Math.*, vol. 28, n° 1, 1975, p. 25-46.

[NAG 90] NAGHDI P. M., « A critical review of the state of finite plasticity », *J. Appl. Math. Phys.*, vol. 41, 1990, p. 315-394.

[NEM 79] NEMAT-NASSER S., « Decomposition of strain measures and their rates in finite deformation elastoplasticity », *Int. J. Solids Structures*, vol. 15, 1979, p. 155-166.

[NEM 82] NEMAT-NASSER S., « On finite deformation elasto-plasticity », *Int. J. Solids Structures*, vol. 18, n° 10, 1982, p. 857-872.

[NGU 00] NGUYEN Q. S., *Stabilité et mécanique non-linéaire*, Hermès, 2000.

[PEE 00] PEEK R., « Axisymmetric wrinkling of cylinders with finite strain », *J. Engineering Mechanics*, vol. 126, n° 5, 2000, p. 455-461.

[RIK 79] RIKS E., « An incremental approach to the solution of snapping and buckling problems », *Int. J. Solids Structures*, vol. 15, 1979, p. 524-551.

[RIK 91] RIKS E., « On formulations on path-following techniques for structural stability analysis », *New advances in computational structural mechanics*, Elsevier, 1991, p. 65-80.

[SEY 94] SEYDEL R., *Practical bifurcation and stability analysis. From equilibrium to chaos*, Springer-Verlag, 1994.

[SIM 86] SIMO J. C., TAYLOR R. L., « A return mapping algorithm for plane stress elastoplasticity », *Int. J. Numerical Methods in Engineering*, vol. 22, 1986, p. 649-670.

[SIN 99] SINGACE A. A., « Axial crushing analysis of tubes deforming in the multi-lobe mode », *Int. J. Mechanical Sciences*, vol. 41, 1999, p. 865-890.

[TVE 83] TVERGAARD V., « On the transition from a diamond mode to an axisymmetric mode of collapse in cylindrical shells », *Int. J. Solids Structures*, vol. 19, n° 10, 1983, p. 845-856.