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

HAL
open science

# Modélisation numérique du flambement élastoplastique d'une poutre sous compression axiale

Philippe Le Grognec, Anh Le Van

► To cite this version:

Philippe Le Grognec, Anh Le Van. Modélisation numérique du flambement élastoplastique d'une poutre sous compression axiale. Revue Européenne des Éléments Finis, 2003, 12 (7-8), pp.897-906. 10.3166/reef.12.897-906. hal-01007075

HAL Id: hal-01007075

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

Submitted on 8 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modélisation numérique du flambement élastoplastique d'une poutre 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É.** Cette étude traite de la modélisation numérique par éléments finis du flambement et du post-flambement lointain d'une poutre élastoplastique sous compression axiale en formulation lagrangienne totale. La plasticité est décrite dans le cadre des matériaux standard généralisés avec le critère de von Mises et un écrouissage isotrope linéaire. Les applications numériques montrent l'efficacité des méthodes employées à décrire systématiquement les principales caractéristiques du flambement plastique, connues en théorie, notamment le continuum de points de bifurcation et les modes correspondants.

**ABSTRACT.** In this paper, the problem of buckling and post-buckling of elastoplastic beams under axial compression in finite displacements and finite strains is resolved by the finite element method with a total Lagrangian formulation. The yield criterion is of the von Mises type with linear isotropic hardening in the framework of the Generalized Standard Materials theory. Numerical computations show the ability of the present formulation to describe any buckling mode in a systematic way and to reproduce numerically the main features of the theoretical plastic buckling problem.

**MOTS-CLÉS:** Bifurcation, flambement, élastoplasticité, poutres sous compression axiale.

**KEYWORDS:** Bifurcation, buckling, elastoplasticity, beams under axial compression.
# 1. Introduction

Le problème de la poutre droite comprimée constitue l'application la plus rencon-
trée dans l'histoire du flambement élastique ou plastique.

En élasticité, le calcul de la première charge critique d'une poutre comprimée re-
monte à la première moitié du dix-huitième siècle avec Euler, et repose sur la résolu-
tion d'un problème aux valeurs propres. Il faudra attendre pourtant 1945 et les travaux
de Koiter pour pouvoir calculer le post-flambement élastique par un développement
asymptotique en série entière et prévoir, de manière similaire, la sensibilité aux im-
perfections d'une telle structure.

En revanche, la résolution du problème de flambement plastique d'une poutre com-
primée est beaucoup moins avancée que dans le cas élastique et a fait encore récem-
ment l'objet de nombreuses études tant sur le plan théorique qu'expérimental. Comme
en élasticité, les premiers pas du flambement plastique concernent l'évaluation de la
première charge critique élastoplastique pour une poutre comprimée. Les premiers tra-
vaux datent de la fin du dix-neuvième siècle et sont entrepris par Engesser, Considère
et, plus tard, von Karman. Cependant, les résultats obtenus jusqu'aux années quarante
ne sont pas tout à fait justes, ou s'avèrent tout au moins mal justifiés. Certains ob-
tiennent la valeur critique du module tangent en négligeant la possibilité de décharge
dans la structure, alors que d'autres supposent que le flambement s'opère à charge
constante, trouvant ainsi la valeur critique du module réduit.

Le premier résultat significatif est dû à Shanley en 1947. Il concerne l'étude d'un
modèle discret, introduit par von Karman en élasticité, et connu aujourd'hui sous le
nom de colonne de Shanley. Ce modèle à deux degrés de liberté (deux ressorts élasto-
plastiques à l'extrémité d'une tige rigide sous compression axiale) représente, de ma-
nière simplifiée, la section d'une poutre et fournit des résultats qualitatifs similaires à
ceux d'une structure élastoplastique continue, en matière de flambement.

Shanley apporte ainsi une réponse satisfaisante quant à la valeur et la nature de la
première charge critique. Il s'agit de la charge critique du module tangent à laquelle
le flambement s'effectue à chargement croissant. Il montre, de plus, que la bifurcation
s'accompagne d'une décharge localisée. Ces deux propriétés essentielles sont aisé-
ment généralisables au cas de la poutre continue.

C'est Hill [HIL 58] qui généralise tous ces résultats pour un milieu tridimension-
nel, contribuant à l'élaboration de critères de stabilité et d'unicité en élastoplasticité
basés sur l'utilisation d'un solide élastique de comparaison, et mettant ainsi en valeur
la dissociation entre les notions de stabilité et de bifurcation.

En matière de post-flambement, Hutchinson propose dans sa synthèse [HUT 74]
une méthode de calcul de la branche post-critique au voisinage de la valeur critique
du module tangent par un développement asymptotique. Le modèle de Hutchinson
(ou modèle de Shanley continu), où les deux ressorts précédents sont remplacés par
une distribution continue, constitue l'application la plus simple de cette méthode.
Il l'applique ensuite aux poutres en compression puis à des structures quelconques,
sous l'hypothèse d'un état prébifurqué homogène. Ce développement asymptotique formel, faisant intervenir des puissances fractionnaires, s'obtient difficilement et ne concerne que le post-flambement immédiat. En outre, on ne dispose d'aucun résultat sur la convergence du développement, l'étendue de sa validité et l'existence même des branches.

Cimetière, El Koulani et Léger [CIM 94], [ELK 96a], [ELK 96b], [ELK 97] sont les premiers à apporter une justification théorique de la validité du développement et de l'existence des branches, en fournissant des résultats de bifurcation globale pour le modèle de Shanley puis pour une poutre comprimée.

Un résultat antérieur essentiel concerne la mise en évidence de continua de points de bifurcation, propres à la plasticité, par Cimetière dans sa thèse [CIM 87] dans le cas de plaques élastoplastiques comprimées. Ce phénomène de spectre continu témoigne de la faculté d'une structure de pouvoir bifurquer à une multitude de valeurs critiques, modifiant en conséquence la zone de décharge mise en jeu au point critique et la rigidité correspondante. Dans l'évolution sur la branche bifurquée, tout au moins à la charge critique du module tangent, la croissance de la zone de décharge augmente progressivement la rigidité globale de la structure et permet à la bifurcation de s'effectuer à charge croissante.

Dans le cas de la poutre comprimée, le développement asymptotique effectué par Hutchinson exprimant la charge en fonction de l'amplitude du mode de bifurcation (c'est-à-dire la projection de la solution sur le mode) est limité au premier terme non linéaire du fait de la complexité de sa mise en œuvre et des coordonnées utilisées liées à la description de l'évolution de la zone de décharge. A la fin des années quatre-vingts, Nguyen et Stolz [NGU 85] et Akel [AKE 87] utilisent une méthode de calcul s'appuyant sur le formalisme des matériaux standard généralisés en développant les deux quantités en séries entières d'un paramètre « temps ».

Plus récemment, Cheng, Lu et Fang [CHE 95] proposent une approche plus simple et applicable également aux plaques et aux coques, leur permettant de calculer le second terme non linéaire dans le cas d'une poutre à section rectangulaire. S'affranchissant des difficultés directement liées à la présence d'une zone de décharge, ils retrouvent les trois premiers termes classiques, déterminent le quatrième et montrent qu'on ne peut aller plus loin sous peine de briser la continuité de la solution, ne donnant plus de sens au développement ainsi défini. L'écriture du principe des travaux virtuels à différents ordres d'un certain paramètre cinématique leur permet de construire une nouvelle branche bifurquée à la charge critique du module tangent avec en particulier une nouvelle charge maximale bien supérieure à celle obtenue dans les travaux antérieurs. Ces auteurs enrichissent plus tard dans [CHE 98] leur développement asymptotique exact précédent en y ajoutant l'effet d'un défaut géométrique. Ils mettent l'accent sur l'importance, tout au moins dans un premier temps, de la décharge élastique comparée aux non-linéarités géométriques. Un défaut croissant modifie la valeur puis provoque la disparition complète de la charge maximale rencontrée dans le cas de la poutre parfaite.
Parallèlement à tous ces développements théoriques, l'absence relative de travaux numériques sur le flambement élastoplastique des poutres s'explique aisément par la quantité importante de résultats analytiques disponibles à ce jour pour un tel problème.

Dans ce travail, on choisit d'aborder le sujet par une approche purement numérique. L'objectif est de pouvoir dégager des méthodes de calcul facilement générables à d'autres types de structure de géométrie plus complexe. En même temps, la modélisation numérique de ce problème doit permettre d'élargir la connaissance et la compréhension du phénomène de flambement d'une poutre élastoplastique en tenant compte davantage des non-linéarités géométriques dans la description du post-flambement lointain.

La première section décrit les différentes notations et hypothèses retenues dans le cadre de notre étude et rassemble les principaux résultats théoriques fondamentaux issus de la littérature, qui serviront de comparaison avec nos propres résultats numériques. La seconde section rappelle les caractéristiques principales de la formulation élément fini retenue. La troisième section présente l'essentiel des méthodes numériques employées et propose une première application du programme éléments finis qui en découle dans le cas élastique. Dans la dernière section, l'utilisation de cet outil numérique sur le même problème en plasticité permet de retrouver, entre autres, les caractéristiques précédentes spécifiques au flambement plastique des structures.

## 2. Cadre théorique

On considère une poutre droite de section constante $S$ et de longueur $l$, bi-appuyée et soumise à une force de compression $\lambda$. Cette poutre est constituée d'un matériau élastoplastique avec une loi d'écrouissage isotrope linéaire. On note respectivement $E$ et $E_T$ le module d'Young et le module tangent supposés tous deux constants. On envisage la flexion de la poutre dans un plan de symétrie et on note $I$ le moment quadratique associé.

L'état prébifurqué, ou la solution fondamentale de ce problème, correspond à une compression homogène élastique puis élastoplastique. Chaque continuum de points de bifurcation (un par mode), observé en plasticité, s'étend de la valeur critique du module tangent $\lambda_T$ à la valeur critique élastique $\lambda_E$ correspondante en passant par une valeur intermédiaire $\lambda_R$ appelée charge critique du module réduit et définie comme étant la valeur à laquelle la bifurcation se fait à charge constante. On s'intéressera ici principalement au premier mode, les résultats concernant les autres modes s'en déduisant facilement.

Afin de simplifier la discussion, on fera les choix suivants sur les propriétés matérielles : on choisit un rapport $E/E_T$ suffisamment faible pour que les différents continua ne se chevauchent pas. De plus, on choisit la contrainte limite élastique $\sigma_0$ suffisamment faible pour que la plastification survienne avant la charge critique du module tangent.
De manière générale, on distingue trois cas possibles quant à l'ordre d'apparition des phénomènes de plasticité et de flambement, dépendant de la position relative de la charge limite élastique $\lambda_0 = \sigma_0 S$ et des diverses charges critiques. Si la limite élastique précède le premier continuum ($\lambda_0 \le \lambda_T$), le flambement (plastique) survient dès la charge critique du module tangent. Si, à l'inverse, $\lambda_0 \ge \lambda_E$, le flambement est purement élastique. Enfin, dans le cas intermédiaire où $\lambda_T \le \lambda_0 \le \lambda_E$, on parle de claquage plastique [AKE 87], et le flambement (plastique) se produit à la charge même de plastification $\lambda_0$. On envisagera seulement le premier cas, caractéristique d'une poutre relativement épaisse, les autres n'apportant rien de nouveau à la description du problème.

Les résultats principaux concernent la bifurcation à la valeur critique du module tangent. La charge $\lambda$ sert de paramètre de bifurcation et on définit le paramètre de perturbation $\xi$ ($\xi \ge 0$) comme la projection de la solution sur le mode de flambement (l'infiniment petit $\xi$ mesure en quelque sorte l'écart entre les solutions fondamentale et bifurquée). Le développement asymptotique général proposé par Hutchinson [HUT 74] pour une structure quelconque est de la forme :

$$ \lambda = \lambda_T + \lambda_1\xi + \lambda_2\xi^{1+\beta} + \dots \quad [1] $$

où $\lambda_1 > 0$, satisfaisant ainsi l'hypothèse de bifurcation à charge croissante, $\lambda_2 < 0$, et l'exposant $\beta$ est rationnel et vérifie $0 < \beta < 1$. La condition $\lambda_2 < 0$ laisse prévoir l'existence d'un maximum pour la charge sur la branche bifurquée, qui limite la résistance de la structure élastoplastique dans son comportement post-critique, d'autant plus que cette charge maximale est proche de la valeur critique du module tangent.

Pour une poutre à section circulaire, la décharge s'initie en un seul point et $\beta = \frac{1}{3}$. Dans le cas d'une section rectangulaire, auquel on s'intéressera uniquement par la suite, c'est toute une ligne qui décharge instantanément en surface et $\beta = \frac{2}{5}$. Notons que dans le cas d'une section symétrique, le développement précédent est encore valable avec $\xi < 0$ en changeant $\xi$ en $-\xi$ dans [1] (prolongeant ainsi $\lambda(\xi)$ par parité). Il n'en est plus de même pour une section non symétrique, comme l'ont souligné Tvergaard et Needleman [TVE 75].

Le développement de Cheng, Lu et Fang [CHE 95], établi pour une section rectangulaire, comprend deux termes non linéaires et s'écrit, en normalisant la charge $\lambda$ par $\lambda_T = \frac{\pi^2 E_T I}{l^2}$ :

$$ \lambda^* = \frac{\lambda}{\lambda_T} = 1 + 3\xi + \lambda_2^* \xi^{\frac{7}{5}} + \lambda_3^* \xi^{\frac{9}{5}} + o(\xi^{\frac{9}{5}}) \quad [2] $$

avec :

$$ \lambda_2^* = -\frac{15}{7} \left( \frac{15\pi E_T}{4\sqrt{2}(E-E_T)} \right)^{\frac{2}{5}} \quad [3] \\ \lambda_3^* = \frac{17}{42} \left( \frac{15\pi E_T}{4\sqrt{2}(E-E_T)} \right)^{\frac{4}{5}} $$

La flèche correspondante est donnée par :

$$ w(x) = W(x)\xi + o(\xi^{\frac{9}{5}}) \quad [4] $$
$W(x)$ étant la fonction propre (demi-sinusoïde s'annulant aux extrémités) normalisée.

Si, dans la pratique, la bifurcation se produit à la charge critique du module tangent, elle peut très bien, en théorie, avoir lieu à partir de n'importe quel point du continuum. Les branches bifurquées partant de ces points sont moins connues et seuls quelques résultats qualitatifs sont disponibles.

La bifurcation s'effectue à charge croissante entre $\lambda_T$ et $\lambda_R = \frac{\pi^2 E_R I}{l^2}$ (avec le module réduit $E_R = \frac{4E E_T}{(\sqrt{E}+\sqrt{E_T})^2}$ dans le cas d'une section rectangulaire) et à charge décroissante entre $\lambda_R$ et $\lambda_E$. Nguyen [NGU 84] démontre, au moyen d'un critère de stabilité basé sur un potentiel plastique, que les branches bifurquées à charge croissante ($\lambda_c \in [\lambda_T, \lambda_R]$) sont stables. Enfin, le comportement des branches bifurquées loin de la solution fondamentale dépendra du modèle choisi et particulièrement des non-linéarités géométriques retenues.

### 3. Formulation numérique

Les calculs numériques suivants s'appuient sur un modèle de coque élastoplastique en grandes rotations développé dans [LEG 03].

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, tous deux indépendants 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 et $\vec{T}(\xi, \eta)$ et $\vec{t}(\xi, \eta)$ les vecteurs directeurs dans les configurations initiale et déformée respectivement.

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}.\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}{\theta} \vec{\Theta} + \frac{1 - \cos \theta}{\theta^2} \vec{\theta} \otimes \vec{\theta} \quad [5] $$

où $\theta = \|\vec{\theta}\|$ désigne l'amplitude de la rotation.
Le tenseur des déformations de Green-Lagrange $\bar{\bar{E}}$ s'écrit :

$$ \bar{\bar{E}} = \frac{1}{2}(\bar{\bar{F}}^T \cdot \bar{\bar{F}} - \bar{\bar{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 [6] $$

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{aligned} 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}) \\ 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{aligned} \quad [7] $$

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

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

$$ \bar{\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 [8] $$

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 $\Sigma$ et l'opérateur tangent consistant $\partial\Sigma/\partial E$ en fonction de l'incrément de déformation, en contraintes planes. La discrétisation en temps des relations en vitesses 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 effectuées par un procédé d'intégration numérique dans l'épaisseur. La programmation a été réalisé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 [9] $$

avec :

$$ \begin{aligned} \eta_{n+1} &= \frac{1}{1 + \frac{2}{3} \Lambda H} \Gamma(\Lambda) C^{-1} \eta^E \\ p_{n+1} &= p_n + \Lambda \sqrt{\frac{2}{3} \eta_{n+1}^T P \eta_{n+1}} \end{aligned} \quad [10] $$

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 [11] $$

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 [12] $$
### 4. Tests numériques préliminaires en élasticité

Avant d'aborder la plasticité, nous allons présenter dans ce paragraphe les résultats concernant la compression axiale d'une poutre bi-appuyée en élasticité. Ces résultats doivent être considérés comme des tests préliminaires en vue de s'assurer de la qualité de l'élément fini utilisé ainsi que du bon fonctionnement des techniques numériques développées, lorsque la poutre travaille en régime élastique.

On considère une poutre bi-appuyée en flexion dans le plan $(Oxz)$, de longueur $l = 1 \text{ m}$, de section rectangulaire, de largeur $b = 5 \text{ cm}$ (selon $y$) et d'épaisseur $t = 1 \text{ mm}$ (selon $z$), figure 1. Une extrémité de la poutre est fixée selon $x$ et l'autre chargée par une force de compression linéique uniforme. La valeur très faible de l'épaisseur permet de mener les calculs numériques à l'aide d'un élément fini coque isoparamétrique à 8 nœuds dont la formulation est décrite dans la section précédente. Un seul élément fini dans la direction $y$ suffit, alors qu'on prend 50 éléments suivant la longueur de façon à rendre compte correctement des grandes rotations. Seul le nœud milieu de chaque extrémité est fixé dans la direction $y$, et non pas l'ensemble des nœuds aux extrémités, pour permettre la dilatation homogène de la poutre. Dans cette configuration, on choisit de prendre 2 points de Gauss dans chacune des directions de la surface de l'élément fini et d'intégrer numériquement dans l'épaisseur avec 4 points de Gauss dans cette direction.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Géométrie, maillage et conditions aux limites</figcaption></figure>

Les propriétés élastiques du matériau sont : module d'Young $E = 2.1 \times 10^{11} \text{ Pa}$ et coefficient de Poisson $\nu = 0.3$.

Le suivi des courbes est effectué avec un pilotage par longueur d'arc sphérique, le facteur variable étant la force de compression. La correction en longueur d'arc né-
cessitant la résolution d'une équation du second degré, une méthode spécifique a été mise en place afin de traiter le cas des racines complexes [LAM 92], [LEG 03]. Le programme de calcul est prévu pour détecter automatiquement les points critiques et pour pouvoir bifurquer correctement sur une branche choisie. Sur la courbe fondamentale, les deux premières charges critiques d'Euler ont été obtenues à $\lambda_{E1} = 8.64 N$ et $\lambda_{E2} = 34.5 N$, soit à moins de 0.05 % des valeurs théoriques $\lambda_{En} = \frac{n^2\pi^2EI}{l^2}$ ($I = \frac{bt^3}{12}$) pour le mode $n$.

On tracera pour chaque courbe bifurquée la force $\lambda$ en fonction du déplacement vertical $w$ en un point d'abscisse $x = 0.3 m$, figure 2 (on évite le milieu de la poutre qui est un point fixe pour le second mode). On observe que la branche bifurquée correspondant au mode 1 n'admet pas de point de bifurcation jusqu'à un post-flambement lointain. La branche bifurquée correspondant au mode 2 admet, quant à elle, des branches secondaires qui correspondent à une déformée hors plan dont l'existence est liée aux conditions aux limites transversales retenues.

Ces premiers résultats, obtenus en élasticité, montrent en particulier que l'élément fini coque utilisé supporte très bien les grandes rotations.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Courbes de réponse et déformées en élasticité</figcaption></figure>
# 5. Poutre élastoplastique bi-appuyée sous compression axiale

## 5.1. Courbes de réponse

On reprend la même poutre que précédemment en considérant cette fois-ci un matériau élastoplastique. On se place dans le cadre des matériaux standard généra- lisés [HAL 75], [GRE 65] avec le critère de von Mises. La contrainte limite élas- tique est $\sigma_0 = 5.10^4 Pa$ et l'écrouissage est supposé linéaire isotrope de module $h = 1.5 10^{11} Pa$. Le module tangent élastoplastique s'en déduit par la relation :

$$E_T = \frac{Eh}{E + h} = 8.75 10^{10} Pa \quad [13]$$

Le rapport $E/E_T$ est choisi dans l'intention d'obtenir des continua de points de bifurcation disjoints et de simplifier l'analyse.

Pour une section rectangulaire, on peut définir le module réduit en fonction du module élastique et du module tangent :

$$E_R = \frac{4EE_T}{(\sqrt{E} + \sqrt{E_T})^2} = 1.29 10^{11} Pa \quad [14]$$

A partir de ces trois différents modules, on peut définir les charges critiques cor- respondantes (respectivement du module tangent, du module réduit et élastique, dans l'ordre croissant) :

$$\lambda_{Tn} = \frac{n^2 \pi^2 E_T I}{l^2}, \quad \lambda_{Rn} = \frac{n^2 \pi^2 E_R I}{l^2}, \quad \lambda_{En} = \frac{n^2 \pi^2 E I}{l^2} \quad [15]$$

ce qui donne, pour le premier mode :

$$\lambda_{T1} = 3.60 N, \quad \lambda_{R1} = 5.30 N, \quad \lambda_{E1} = 8.64 N \quad [16]$$

Si chaque point du continuum est un point de bifurcation potentiel, c'est qu'il existe en chacun d'eux une solution en vitesses non triviale, qui correspond non plus à la poutre entièrement plastifiée, mais à la poutre partiellement déchargée, la forme de la zone de décharge variant selon la valeur atteinte de la force de compression. A la charge critique du module tangent, la notion du solide élastique de comparaison s'applique et tout se déroule de la même manière qu'en élasticité : le programme de calcul détecte la charge critique $\lambda_T$ lorsque la matrice raideur devient singulière et la prédiction de la branche bifurquée associée (qui s'initie avec un seul point de dé- charge) est faite avec le vecteur propre correspondant. Pour bifurquer à un autre point du continuum, nous prenons comme prédiction le vecteur propre associé à la valeur propre de plus petit module, à défaut de connaître le bon champ des vitesses corres- pondant à la forme exacte de la zone de décharge. Les résultats numériques montrent que ce procédé permet de bifurquer sur l'ensemble du continuum. En concordance
avec la théorie, les branches trouvées sont comprises entre celles issues des charges critiques du module tangent $\lambda_T$ et élastique $\lambda_E$. Les courbes de la figure 3 représentent la flèche au milieu de la poutre en fonction de l'effort de compression. Les résultats numériques obtenus montrent que la bifurcation s'effectue à charge croissante pour $\lambda_c \in [\lambda_T, \lambda_R]$, à charge décroissante pour $\lambda_c \in [\lambda_R, \lambda_E]$ et à charge constante en $\lambda_R$. Ces comportements sont en accord avec la théorie. En poussant les calculs jusqu'aux déformées assez importantes, on obtient les courbes de la figure 4. On peut y distinguer trois zones :

– la zone post-critique immédiate (petites flèches) où c'est la non-linéarité matérielle qui l'emporte sur celle géométrique;

– la zone post-critique lointaine (flèches relativement grandes) où c'est l'inverse qui se produit;

– et une zone intermédiaire où les deux non-linéarités ont des effets équivalents.

Si on se restreint à la zone post-critique immédiate, correspondant à $w < 0.0005$ sur la figure 3, on peut considérer que les branches bifurquées tendent toutes vers une asymptote verticale d'équation $\lambda = \lambda_R$, ce qui est conforme à la solution analytique d'El Koulani [ELK 96a]. Par contre, la figure 4 montre que dans la zone post-critique lointaine, les branches bifurquées tendent vers une autre asymptote verticale. Il est normal que la solution théorique précitée, établie dans l'hypothèse de petites perturbations, ne prévoit pas cette nouvelle asymptote dans cette zone où ce sont les non-linéarités géométriques qui sont prépondérantes.

Analysons plus en détail la branche bifurquée issue de la charge critique du module tangent qui correspond à une bifurcation avec charge plastique totale de la poutre. La figure 5 regroupe quatre courbes post-critiques :

– la courbe analytique ($\lambda, w$) correspondant au développement asymptotique [2] en ne conservant que le seul terme linéaire;

– la même courbe avec le terme non linéaire à la puissance $\frac{7}{5}$, obtenu par Hutchinson [HUT 74];

– la même courbe avec les deux termes non linéaires – à la puissance $\frac{7}{5}$ et $\frac{9}{5}$ – obtenus par Cheng, Lu et Fang [CHE 95];

– et enfin la courbe issue des résultats numériques.

La comparaison de ces différentes courbes permet de juger de l'importance des termes successifs du développement asymptotique. La courbe avec le seul terme linéaire n'admet pas de maximum en charge. La prise en compte d'un terme non linéaire dans le développement [2] permet de rendre compte de la charge maximale supportable par la structure, mais la valeur de cette charge s'avère trop conservative. Les résultats numériques sont très proches de ceux théoriques avec deux termes non linéaires, notamment au niveau de la charge maximale. Les derniers points s'écartent de la courbe théorique, ce qui peut s'expliquer par le fait que d'une part la théorie utilise un développement asymptotique tronqué et que d'autre part elle ne tient pas compte
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Courbes de réponse en plasticité – post-flambement immédiat</figcaption></figure>

des non-linéarités géométriques qui deviennent prépondérantes lorsque la flèche *w* est assez grande. En tous cas, lorsque les non-linéarités géométriques n'ont encore que très peu d'effet, les résultats numériques sont pratiquement confondus avec ceux théoriques de Hutchinson et de Cheng, Lu et Fang, et se rapprochent très nettement de la courbe analytique monotone d'El Koulani [ELK 96a].

**5.2. Zones de décharge**

Intéressons-nous maintenant à la topologie des zones de décharge dans la poutre aux différents points de bifurcation du continuum et à leur évolution sur la branche bifurquée correspondante.

On étudie toujours une poutre bi-appuyée de longueur *l* = 1 m, mais afin de pouvoir visualiser les zones de décharge on prend cette fois-ci la largeur *b* = 10 cm selon *y*, l'épaisseur *t* = 1 cm selon *z*, et on considère la flexion dans le plan (Oxy), figure 6. Ainsi, on fixe tous les degrés de liberté de rotation et ceux de translation suivant l'axe (Oz) orthogonal au plan de flexion.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Courbes de réponse en plasticité – post-flambement lointain</figcaption></figure>

Le matériau utilisé reste inchangé, excepté $\sigma_0$ qui passe à la valeur de $5.10^7$ Pa, de manière à assurer la plastification avant la bifurcation au module tangent. Le maillage est enfin adapté à la nouvelle géométrie (figure 6) avec 30 éléments selon $x$ et 4 éléments selon $y$. On utilise 7 points de Gauss par élément selon $x$ et $y$. Le nombre de points de Gauss servant à l’intégration numérique selon $z$ est fixé à 2 : il n’a en effet nul besoin d’être important, compte tenu du comportement uniforme suivant cette direction.

### 5.2.1. Bifurcation au module réduit

Les résultats théoriques montrent que la zone de décharge observée à la bifurcation passe d'un domaine convexe à un domaine concave lorsque la charge critique décrit le continuum $[\lambda_T, \lambda_E]$ en passant par la valeur de transition égale à la charge critique du module réduit $\lambda_R$ qui est par conséquent caractérisée par une zone de décharge rec-tangulaire (ou plusieurs morceaux rectangulaires s'il ne s'agit pas du premier mode). On calcule ici les branches bifurquées aux deux premières charges critiques du mo-dule réduit, $\lambda_{R1}$ et $\lambda_{R2}$. Connaissant les valeurs théoriques $\lambda_{R1} = 1.06 \ 10^6$ N et $\lambda_{R2} = 4.24 \ 10^6$ N, on impose manuellement la bifurcation à ces valeurs.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Bifurcation au module tangent : comparaison développements asymptotiques et résultats numériques</figcaption></figure>

D'après la théorie, l'ordonnée de la frontière horizontale séparant les zones de charge plastique et de décharge élastique, comptée à partir de la ligne moyenne, s'écrit $d_R = Y_R \frac{b}{2}$, où $Y_R$ vaut pour une poutre à section rectangulaire $Y_R = \pm \frac{1-\sqrt{\gamma}}{1+\sqrt{\gamma}}$ avec $\gamma = E_T/E$. La distance $d_R$ (ou $Y_R$) est en fait racine d'un polynôme du second degré et son signe dépend du sens vers lequel la poutre va flamber. En considérant des flèches négatives, la zone de décharge est dans la partie inférieure de la poutre et $d_R = -0.01075$.

En théorie, on distingue les points en charge plastique de ceux en décharge en considérant la vitesse de déformation plastique instantanée. Ici dans le calcul numérique, la détection des points en charge plastique est basée sur la comparaison des déformations plastiques équivalentes entre deux incréments successifs. Dans la mesure où la taille de la zone de décharge est croissante (du moins en post-critique proche), ce procédé, appliqué à deux incréments $n$ et $n + 1$, permet d'obtenir l'information voulue à l'incrément $n$.

Les résultats numériques montrent que la frontière entre zones de charge plastique et de décharge élastique est bien horizontale sur la majeure partie de la longueur de la poutre (figure 7). On relève l'ordonnée de cette frontière pour les deux premiers modes
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Géométrie, maillage et conditions aux limites pour l'étude des zones de décharge</figcaption></figure>

et on obtient respectivement $d_{R1} = -0.00995$ et $d_{R2} = \pm 0.0053$ ($\pm$ par symétrie), soit une erreur respective vis-à-vis de la solution analytique de l'ordre de 1 % et 5 % de l'épaisseur $b$ de la poutre.

### 5.2.2. Bifurcation au module tangent

On calcule ici la branche bifurquée à la première charge critique du module tangent $\lambda_{T1}$. Contrairement aux bifurcations aux modules réduits, celle à $\lambda_{T1}$ se fait de manière automatique.

A la charge critique du module tangent, la zone de décharge est réduite à un seul point à partir duquel se développe une zone de taille croissante le long de la branche bifurquée. Notons $a$ la demi-longueur de la zone de décharge (symétrique par rapport au plan médian de la poutre) le long de la couche supérieure (figure 8). On définit $\tau$ le rapport entre l'étendue en surface de la zone de décharge et la longueur de la poutre, soit $\tau = \frac{2a}{l}$. Le scalaire $d$ représente l'ordonnée de la frontière élastique-plastique au milieu de la poutre. La forme de la zone de décharge est décrite par l'équation suivante, paramétrée par $a$ [HUT 74], [NGU 00] :

$$ 
\begin{aligned}
\lambda &= \lambda_T + \lambda_5 \frac{\tau^5}{5!} \\
d &= \frac{b}{2} \left( 1 - \frac{\pi^2}{8} \tau^2 \right)
\end{aligned} 
\qquad [17] $$

$$ \text{avec } \lambda_5 = \frac{3\pi^6(E-E_T)I}{4l^2} \quad (I = \frac{b^3t}{12}). $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Zones de décharge à la bifurcation au module réduit (modes 1 et 2)</figcaption></figure>

La figure 9 visualise les zones de décharge analytique et numérique pour une même valeur de chargement. On peut constater que la zone de décharge théorique est très bien représentée par les calculs numériques.

L'évolution théorique de l'ordonnée *d* en fonction de la charge λ sur la branche post-critique est représentée en trait continu sur la figure 10. La courbe de réponse numérique est quant à elle obtenue en relevant à différents incréments de la branche bifurquée les valeurs numériques de λ et *d* (pour le relevé des valeurs de l'ordonnée de la frontière, la précision dépend bien sûr du nombre de points de Gauss retenu dans la direction *y*). Pour une valeur de *d* donnée, l'écart entre les deux courbes est de 3 % en moyenne par rapport à la force de compression théorique.

**6. Conclusion**

On s'est intéressé dans ce travail à la résolution numérique du problème de flam-
bement d'une poutre élastoplastique sous compression axiale. On a retenu le critère
de von Mises avec un écrouissage isotrope linéaire dans le cadre des matériaux stan-
dard généralisés en formulation lagrangienne totale. Le modèle élément fini associé,
le pilotage par longueur d'arc et les méthodes de bifurcation mis en œuvre ont permis
de mettre en valeur les spécificités suivantes du flambement plastique relativement
connues en théorie :
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Caractéristiques géométriques de la zone de décharge sur la branche issue de la bifurcation au module tangent (mode 1)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Zones de décharge sur la branche issue de la bifurcation au module tangent (mode 1)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evolution de la profondeur de la zone de décharge sur la branche issue de la bifurcation au module tangent (mode 1)</figcaption></figure>

– On a mis en évidence par le calcul numérique l’existence d’un continuum de points de bifurcation à chaque mode allant de la charge critique du module tangent à la charge critique élastique avec une bifurcation à charge soit croissante, soit décroissante (constante pour le cas limite de la charge critique du module réduit).

Il convient de remarquer que pour simplifier l’analyse des résultats, nous avons voilà :
- lontairement choisi d’une part, un rapport *E*/*E*<sub>T</sub> suffisamment faible pour que les différents continua ne se chevauchent pas et, d’autre part, une contrainte limite élastique *σ*<sub>0</sub> suffisamment faible pour que la plastification survienne avant la charge critique du module tangent.

Concernant la bifurcation au module tangent, la courbe de réponse post-critique calculée numériquement s'accorde très bien avec le développement asymptotique correspondant de Cheng, Lu et Fang [CHE 95], avec en particulier une charge maximum très proche entre les deux cas. Pour l'ensemble des branches bifurquées, nous avons obtenu – en petits déplacements – la même asymptote analytique d'El Koulanî [ELK 96a]. Pour les déplacements plus importants, nous avons mis en évidence une autre asymptote du fait de la prise en compte de l'ensemble des non-linéarités géométriques dans notre modèle.

– On a observé enfin l’initiation et la progression des zones de décharge dans la poutre le long des branches bifurquées, propres au flambement plastique des struc-
tures. La comparaison de la forme et de la taille des zones de décharge obtenues numériquement, en différents points du continuum et sur la branche issue de la bifurcation au module tangent, avec les résultats théoriques de Nguyen [NGU 00] s'est avérée satisfaisante.

Les résultats obtenus nous permettent d’envisager l’utilisation de ce programme dans des cas plus complexes (le problème du cylindre sous compression axiale a déjà fait l’objet d’une étude numérique approfondie [LEG 03]).

Une autre perspective envisageable sur la poutre sous compression axiale concerne l’analyse numérique de la sensibilité aux imperfections, qui viendra compléter la présente étude à l’instar de ce qui est fait dans les travaux théoriques de Cheng, Fang et Lu [CHE 98]. D’un point de vue numérique, l’ajout d’un défaut – géométrique par exemple – ne pose pas de difficulté technique et doit permettre d’évaluer la chute non négligeable de la charge limite applicable à la structure. Si les imperfections sont évitables en pratique, l’étude que nous avons menée sur une poutre parfaite n’est pour autant pas dénuée d’intérêt dans la mesure où elle fournit le squelette des courbes de réponse de la structure avec défaut.

## 7. Bibliographie

[AKE 87] AKEL S., « Sur le flambage des structures élastoplastiques », *Thèse de doctorat*, Ecole Nationale des Ponts et Chaussées, 1987.

[CHE 95] CHENG Y. S., LU W. D., FANG H., « Plastic post-buckling of a simply supported column with a solid rectangular cross-section », *Applied Mathematics and Mechanics*, vol. 16, n° 8, 1995, p. 713-722.

[CHE 98] CHENG Y. S., FANG H., LU W. D., « Imperfection sensitivity analysis of a rectangular column compressed into the plastic range », *Applied Mathematics and Mechanics*, vol. 19, n° 1, 1998, p. 9-14.

[CIM 87] CIMETIÈRE A., « Sur la modélisation et le flambage des plaques élastoplastiques », *Thèse de doctorat*, Université de Poitiers, 1987.

[CIM 94] CIMETIÈRE A., EL KOULANI A., LÉGER A., « Flambage naissant et post-flambage d'un modèle simple en élastoplasticité », Paris, 1994, C.R. Acad. Sc., t. 319.

[ELK 96a] EL KOULANI A., « Continuation dans les problèmes à frontières libres de type bifurcations plastiques », *Thèse de doctorat*, Université de Metz, 1996.

[ELK 96b] EL KOULANI A., LÉGER A., « Solutions bifurquées du problème en vitesses initiales pour une poutre élastoplastique », Paris, 1996, C.R. Acad. Sc., t. 322.

[ELK 97] EL KOULANI A., LÉGER A., « Comportement post-critique des poutres élastoplastiques : existence et régularité des branches bifurquées », Paris, 1997, C.R. Acad. Sc.

[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.

[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.
[HIL 58] HILL R., « A general theory of uniqueness and stability in elastic-plastic solids », *J. Mech. Phys. Solids*, vol. 6, 1958, p. 236-249.

[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.

[LEG 03] LE GROGNEC P., LE VAN A., « Bifurcations primaires et secondaires d'un cylindre sous compression axiale », *Revue Européenne des Eléments Finis*, vol. 12, n° 1, 2003, p. 7-41.

[NGU 84] NGUYEN Q. S., « Bifurcation et stabilité des systèmes irréversibles obéissant au principe de dissipation maximale », *Journal de Mécanique théorique et appliquée*, vol. 3, n° 1, 1984, p. 41-61.

[NGU 85] NGUYEN Q. S., STOLZ C., « Sur la méthode de développement asymptotique en flambage plastique », Paris, 1985, C.R. Acad. Sc., t. 300.

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

[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.

[TVE 75] TVERGAARD V., NEEDLEMAN A., « On the buckling of elastic-plastic columns with asymmetric cross-sections », *Int. J. Mechanical Sciences*, vol. 17, 1975, p. 419-424.