# Some new analytical results for plastic buckling and initial post-buckling of plates and cylinders under uniform compression

Philippe Le Grognec, Anh Le Van

► To cite this version:

Philippe Le Grognec, Anh Le Van. Some new analytical results for plastic buckling and initial post-buckling of plates and cylinders under uniform compression. Thin-Walled Structures, 2009, 47 (8-9), pp.879-889. 10.1016/j.tws.2009.02.010. hal-01007313

HAL Id: hal-01007313

https://hal.science/hal-01007313

Submitted on 5 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 Attribution 4.0 International License
# Some new analytical results for plastic buckling and initial post-buckling of plates and cylinders under uniform compression

Philippe Le Grognec<sup>a</sup>, Anh Le van<sup>b</sup>

This paper is devoted to the theoretical aspects of the elastoplastic buckling and initial post-buckling of plates and cylinders under uniform compression. The analysis is based on the 3D plastic bifurcation theory assuming the $J_2$ flow theory of plasticity with the von Mises yield criterion and a linear isotropic hardening. The proposed method is shown to be a systematic and unified way to obtain the critical loads, the buckling modes and the initial slope of the bifurcated branch for rectangular plates under uniaxial or biaxial compression(-tension) and cylinders under axial compression, with various boundary conditions.

## 1. Introduction

Failure of thin structures which is mainly due to the buckling phenomenon implies the analysis of buckling and post-buckling behaviors for their mechanical design, namely the calculation of the critical loads, the bifurcation modes and the post-critical equilibrium branches.

In elasticity, the bifurcation is related to the structural instability as shown in Koiter's theory. The problem is much more difficult with thicker structures where plasticity may occur before buckling. Shanley [1] was the first to give the tangent modulus critical load for a discrete model. He showed that the load increases during the post-bifurcation stage, owing to the elastic unloading in the structure. Hill [2] extended these results to a 3D continuum by using the concept of "comparison elastic solid". He examined the uniqueness and stability criteria, and pointed out the difference between bifurcation and stability. In order to study the plastic post-buckling, Hutchinson [3] extended Koiter's results in elasticity by using a series expansion with fractional exponents, assuming an homogeneous pre-bifurcation state. In the 1980s, Nguyen and Stolz [4] and Akel [5] reformulated the problem within the framework of the generalized standard materials and expanded both the load and the perturbation (i.e. the bifurcation mode amplitude) parameters in power series of a time-like parameter. They showed that their power expansion

yields the same expansion with fractional exponents as obtained by Hutchinson. Much later, Cimetière et al. [6–9] thoroughly solved the bifurcation problem for the Shanley discrete model and a compressed beam, and first provided the necessary theoretical ingredients such as the validity, the convergence of the previous expansion and the existence of the post-critical branches. Another significant result is the existence of continua of bifurcation points in plastic buckling problems, which was discovered by Cimetière [10] when dealing with the case of compressed rectangular plates. These continua enable the structure to bifurcate within intervals of critical loads, by continuously modifying the unloaded zone and the structural stiffness.

The main results were obtained for the plastic bifurcation at the tangent modulus critical value. The load $\lambda$ acts as a bifurcation parameter; the perturbation parameter denoted by $\xi$ ($\xi \ge 0$) can be viewed as a measure of the difference between the fundamental and post-buckling solutions. The asymptotic expansion proposed by Hutchinson [3] for a general structure reads

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

where $\lambda_1$ is positive, thus verifying the hypothesis of an increasing bifurcation load, $\lambda_2$ is negative and the exponent $\beta$ is a rational number such that $0 < \beta < 1$. The $\lambda_2 < 0$ condition leads one to expect a load limit point on the bifurcated branch, i.e. an upper bound for the strength of the elastoplastic structure in the post-critical behavior.

From that time, the main theoretical results available in the literature concern the calculation of the critical loads $\lambda_T$ and the

<footnote><sup>a</sup> Ecole des Mines de Douai, Polymers and Composites Technology and Mechanical Engineering Department, 941 rue Charles Bourseul, BP 10838, 59508 Douai Cedex, France</footnote>

<footnote><sup>b</sup> GeM (Laboratory of Civil and Mechanical Engineering), Faculty of Science, University of Nantes, 2 rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France</footnote>
corresponding buckling modes for different geometries and loading conditions.

The case of a rectangular plate was investigated by a few authors. Among others, Durban and Zuckerman [11] derived semi-analytical solutions for the plastic buckling of rectangular plates under uniform biaxial compression/tension with various boundary conditions. They particularly showed the possible discrepancy between the results provided by the $J_2$ flow and deformation theories, thus illustrating the well-known plastic buckling paradox.

One of the most considered cases in plastic buckling is the axial compression of circular cylindrical shells. In this context, the two following aspects are usually analyzed: the influence of the constitutive law and the effect of the boundary conditions.

Except for very long cylinders which buckle like beams, two types of analytical deformed shapes can occur in the post-buckling behavior of cylinders under axial compression. In the case of thin shells, the buckling mode may usually be non-axisymmetric, whereas for moderately thick shells, when plastic strains appear before buckling, axisymmetric modes are most often obtained.

In the latter case, Batterman [12] derived the critical stress from the rate equilibrium equations, using the flow theory of plasticity:

$$ \sigma_c^p = \frac{Ee}{R} \sqrt{\frac{4H}{3(5 - 4\nu)E + 4H(1 - \nu^2)}} \quad (2) $$

where e denotes the shell thickness, R the radius of the cylinder, and E, ν and H represent the Young's modulus, the Poisson's ratio and the linear hardening modulus, respectively. The critical stress (2) was obtained for a cylinder with simply supported ends. It can be shown that built-in edges almost yield the same critical value and mode shape. In order to study the influence of the boundary conditions, Batterman [13] analytically solved a semi-infinite cylindrical shell with free ends, which should be the most disadvantageous case. With a constant tangent modulus, he obtained a critical stress equal to half of the value for simply supported edges. Batterman also showed that the critical stress is a little bit smaller for finite shells, yet it hardly depends on the cylinder length. Furthermore, the axisymmetric buckling mode looks like an exponential wave, the radial displacement being maximal at the ends of the cylinder and almost zero in the middle.

Ore and Durban [14] also derived semi-analytical values for the critical load of a cylinder under axial compression and various boundary conditions, with a special emphasis on axisymmetric modes. They showed again the discrepancy between the results provided by the flow and deformation theories of plasticity. More recently, Mao and Lu [15] considered the more general case of a cylindrical shell under combined axial and circumferential compressive loading stresses. The transverse shear is taken into account in order to analyze properly the interactive effect of stresses in both directions on the plastic buckling.

Let us finally mention other analytical approaches for this post-buckling problem, involving geometric considerations. From a mathematical point of view, the post-critical deformed shape— whether the buckling is local or global—can be seen as an isometric transformation of the initial cylindrical geometry. Accordingly, it is possible to reach the final deformed shape and evaluate the limit load without describing the whole equilibrium curve. In a slightly different way, some authors—among which Singace [16] and Gupta and Abbas [17]—represented the post-buckling of the cylinder by a cyclic formation of folds giving rise to plastic hinges and investigated the location and the shape of these folds in both axisymmetric and non-axisymmetric cases. Within the scope of the limit analysis, which means an advanced post-buckling state, the obtained solutions were found to be in very good agreement with experimental results [18].

The interested reader can find more details on the theoretical developments and more complete states-of-art in the quoted references.

## 2. Outline of the paper

The aim of this work is to obtain analytical results for the elastoplastic buckling and initial post-buckling of compressed plates and cylinders, namely (i) the critical load, the bifurcation mode and (ii) the initial slope of the bifurcated branch. The proposed method provides a systematic and unified way to obtain the sought quantities since they will all be derived from the 3D elastoplastic bifurcation framework.

Section 3 summarizes the main known results for the 3D plastic bifurcation assuming the $J_2$ flow theory with the von Mises yield criterion and a linear isotropic hardening. In Section 4, the 3D results are particularized to the case of a uniform biaxial stress state with the purpose of treating the plastic bifurcation of plates and cylinders.

Using adequate assumptions for plate and shell models then enables one to apply the previous developments to a rectangular plate under uniaxial or biaxial compression(-tension) in Section 5 and a cylinder under axial compression in Section 6, with various boundary conditions in both geometries.

It should be noted that the initial slope of the bifurcated branch gives information about the load variation in the vicinity of the critical point and is sufficient for the initial post-buckling analysis. The study of advanced post-buckling necessitates higher-order terms like $\lambda_2$ in expansion (1) and is beyond the scope of this work.

## 3. The 3D plastic bifurcation theory

It is essential in the bifurcation analysis to distinguish between the reference and the current configuration so as to take into account the geometric changes and correctly derive the terms responsible for the bifurcation phenomenon. For this purpose, the total Lagrangian formulation is usually adopted since it proves to be a convenient framework in which the bifurcation equation can easily be linearized.

In the 3D context, let us consider a 3D body with reference volume $\Omega_0$; the position of the current particle in the reference configuration will be denoted by $X \in \Omega_0$. The kinematic decomposition is the first step to be defined in a finite plasticity analysis. In general, one has to choose among the different kinematic descriptions proposed in the literature, the most widely used being the Lee multiplicative decomposition of the deformation gradient [19,20]. However, as the structural elements of interest here are plates and cylinders, it is assumed that the pre-critical deformations are small—see assumption (20) below—so that all the kinematic decompositions are equivalent in the present bifurcation analysis and the choice of a kinematic decomposition has no significant consequences. In order to deal with relations analogous to those in small strains, use will be made of the Green and Naghdi additive decomposition [21,22], which splits the Green strain tensor E additively into its elastic and plastic parts:

$$ E = E^e + E^p \quad (3) $$

$E^e$ (resp. $E^p$) is the so-called elastic (resp. plastic) strain tensor. As mentioned by Green and Naghdi [21], here the “elastic” tensor $E^e$ is not the classical elastic tensor, it is named so just by analogy with the case of small deformations. In this approach, the total
strain **E** and the plastic one **E**^p are considered as independent variables, whereas the elastic strain **E**^e is viewed merely as the difference between **E** and **E**^p [23–25]. The plastic strain **E**^p is assumed to be a symmetric and objective tensor; and so is **E**^e.

The description in the sequel is based on the generalized
standard materials theory [26,27]. Let us assume that the free
energy by unit volume w in the reference configuration can
be written as the sum of the elastic energy w^e, depending on the
elastic strain E^e, and the hardening energy w^x, depending on
the hardening parameter α. Since we confine ourselves to small
elastic strains during the whole process, the elastic energy can be
expressed as a quadratic function of the elastic strain. The free
energy by unit volume then writes

$$
\begin{align}
w(\mathbf{E}^e, \alpha) &= w^e(\mathbf{E}^e) + w^x(\alpha) \nonumber \\
&= \frac{1}{2} \mathbf{E}^e : \mathbf{D} : \mathbf{E}^e + w^x(\alpha) \nonumber \\
&= \frac{1}{2} \varepsilon_{ji} D_{ijkl} E_{lk}^x + w^x(\alpha) \tag{4}
\end{align}
$$

where **D** represents the fourth-order elasticity tensor which is assumed to be constant. The material is isotropic so that the elastic coefficients in an orthonormal basis are $D_{ijkl} = \bar{\lambda}\delta_{ij}\delta_{kl} + \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj})$, where $\delta_{ij}$ is the Kronecker symbol, $\bar{\lambda}$ and $\mu$ are the Lamé constants. Use will also be made of Young's modulus $E$ and Poisson's ratio $\nu$ which are related to $\bar{\lambda}$ and $\mu$ by $\bar{\lambda} = \frac{E\nu}{1+(\nu-1)\bar{\lambda}}$.

With the additive decomposition (3), the state laws and the evolution laws are analogous to those in small strains. The state laws express the second Kirchhoff stress tensor Σ (symmetric) and the thermodynamical force A (conjugate of the hardening parameter α) as the derivatives of the free energy with respect to the elastic strain and the hardening variable, respectively:

$$
\Sigma = \frac{\partial w}{\partial \mathbf{E}^e} = \mathbf{D} : \mathbf{E}^e, \quad A = \frac{\partial w}{\partial \alpha} \qquad (5)
$$

We assume that the elastic region is defined by the von Mises
yield function with an isotropic hardening:

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

where $\Sigma^d$ means the deviatoric part of the stress tensor $\Sigma$ and $\sigma_0$ is the yield stress.

In the framework of generalized standard materials [26], the normality rules state that the plastic strain **E**^p as well as the hardening variable α derive from a dissipative potential, which is equal to the yield function:

$$
\dot{\mathbf{E}}^p = \dot{\lambda} \frac{\partial f}{\partial \mathbf{\Sigma}}, \quad \dot{\alpha} = -\dot{\lambda} \frac{\partial f}{\partial A} \qquad (7)
$$

where $\dot{\lambda}$ is the plastic multiplier. Note that the rates of the Lagrangian variables $\mathbf{E}^p$ and $\alpha$ are ordinary derivatives with respect to time. The consistency condition leads to the following expression for the plastic multiplier:

$$
\dot{\lambda} = \frac{\frac{\partial f}{\partial \boldsymbol{\Sigma}} : \mathbf{D} : \dot{\mathbf{E}}}{H + \frac{\partial f}{\partial \boldsymbol{\Sigma}} : \mathbf{D} : \frac{\partial f}{\partial \boldsymbol{\Sigma}}} \quad (8)
$$

where $H = (\partial f / \partial A)^2 \partial^2 w / \partial \alpha^2$ is the hardening modulus. The material tangent elastoplastic tensor is

$$
\mathbf{D}^p = \frac{\partial \boldsymbol{\Sigma}}{\partial \mathbf{E}} = \mathbf{D} - \frac{\mathbf{D} : \frac{\partial f}{\partial \boldsymbol{\Sigma}} \otimes \frac{\partial f}{\partial \boldsymbol{\Sigma}} : \mathbf{D}}{H + \frac{\partial f}{\partial \boldsymbol{\Sigma}} : \mathbf{D} : \frac{\partial f}{\partial \boldsymbol{\Sigma}}} \quad (9)
$$

where the tensor product ⊗ of two second-order tensors **S** and **T** is
defined by (S ⊗ T)_{ijkl} = S_{ij}T_{kl}. Relation (9) can be recast as

$$
\mathbf{D}^p = \mathbf{D} - \mathbf{N} \otimes \mathbf{N} \qquad (10)
$$

where the symmetric tensor **N** is

$$
\mathbf{N} = \frac{\mathbf{D} : \frac{\partial f}{\partial \boldsymbol{\Sigma}}}{\sqrt{H + \frac{\partial f}{\partial \boldsymbol{\Sigma}} : \mathbf{D} : \frac{\partial f}{\partial \boldsymbol{\Sigma}}}} \quad (11)
$$

One then derives the nominal tangent elastoplastic tensor:

$$
\mathbf{K}^p = \frac{\partial \mathbf{\Pi}}{\partial \mathbf{F}} = \mathbf{F} \cdot \frac{\partial \boldsymbol{\Sigma}}{\partial \mathbf{E}} \cdot \mathbf{F}^T + (\| \cdot \| \cdot \boldsymbol{\Sigma})^T \quad (12)
$$

In the previous equation, **F** stands for the deformation gradient, **Π** = **F** · **Σ** the first Kirchhoff (non-symmetric) stress tensor, ‖ represents the fourth-order unit tensor (‖*ijkl* = δᵢⱼδₖₗ) and the superscript T the major transposition of a fourth-order tensor defined by (*D*^T)ᵢₖₗ = *D*ᵢₖₗ. Like (10), Eq. (12) can be rewritten as

$$
\mathbf{K}^p = \mathbf{K}^e - \mathbf{M}^T \otimes \mathbf{M} \quad (13)
$$

with

$$
\mathbf{K}^e = \mathbf{F} \cdot \mathbf{D} \cdot \mathbf{F}^T + (\| \cdot \| \cdot \mathbf{\Sigma})^T, \quad \mathbf{M} = \mathbf{N} \cdot \mathbf{F}^T \quad (14)
$$

Now, let us assume that there exists a fundamental equilibrium
path λ ↦→ **U**<sub>f</sub>(λ), which is the displacement solution of the
elastoplastic problem under an increasing load factor λ (λ̇ > 0),
and that there exists a bifurcated solution **U** at a critical time t<sub>c</sub>
corresponding to the critical load factor λ<sub>c</sub> = λ(t<sub>c</sub>).

The bifurcated branch is described by an asymptotic expansion
with respect to the perturbation parameter ξ ≥ 0:

$$
\left\{
\begin{array}{l}
\lambda = \lambda_c + \lambda_1 \xi + o(\xi) \\
\mathbf{U} = \mathbf{U}_f(\lambda) + \xi \mathbf{X} + o(\xi)
\end{array}
\right.
\qquad (15)
$$

**Assumption 1.** At critical time $t_c$, the plastic zone $\Omega^p(t_c)$ corresponding to the bifurcated solution is equal to that of the fundamental solution $\Omega_f^p(t_c)$, i.e. the bifurcation takes place at the tangent modulus critical load with incipient unloading.

With the above assumption, the critical load $\lambda_c = \lambda_T$ and the bifurcation mode **X** are obtained by solving the following bifurcation equation [28,29]:

$$
\forall \delta \mathbf{U}, \quad \int_{\Omega_0} \nabla^T \delta \mathbf{U} : \mathbf{K}^i (\mathbf{U}_f(\lambda_c)) : \nabla \mathbf{X} d\Omega_0 = 0 \qquad (16)
$$

where $\mathbf{K}^i$ stands for $\mathbf{K}^e$ or $\mathbf{K}^p$ (see Eqs. (13) and (14)), depending on whether the current point belongs to the elastic or plastic zone, respectively.

The coefficient $\lambda_1$ in expansion (15) represents the initial slope of the bifurcated branch at the critical point. Its expression involves some more assumptions, which are listed below.

**Assumption 2.** $\lambda_c$ is a simple eigenvalue of Eq. (16).

**Assumption 3.** The so-called transversality coefficient *Tr* is positive:

$$
\mathrm{Tr} \doteq - \int_{\Omega_0} \nabla^T \mathbf{X} : \frac{d\mathbf{K}^i(\mathbf{U}_f(\lambda))}{d\lambda} \bigg|_{\lambda=\lambda_c} : \nabla \mathbf{X} d\Omega_0 > 0 \quad (17)
$$

**Assumption 4.** The plastic zone corresponding to the fundamental solution at the critical time $\Omega_f^p(t_c)$ is active (i.e. there is loading and no neutral loading).

**Assumption 5.** ∃Δ > 0, ∀X ∈ Ω_f^p(t_c), **M**(X, λ_c) : ∇**U**<sub>f,λ</sub>(X, λ_c) ≥ Δ > 0.

Under these assumptions, one has an angular bifurcation and
the initial slope $\lambda_1$ is given by

$$
\lambda_1 = \max_{X \in \Omega_f^p(t_c)} \left( \frac{\mathbf{M}(X, \lambda_c) : \nabla \mathbf{X}}{\mathbf{M}(X, \lambda_c) : \nabla \mathbf{U}_{f,\lambda}(X, \lambda_c)} \right) \quad (18)
$$
### 4. Biaxial stress state

In this section, we shall derive the expressions for the elastoplastic tangent moduli $D^p$ and $K^p$ (Eqs. (10) and (13)) in the particular case of biaxial stress states. In the next sections, these expressions will be applied to the problem of rectangular plates under uniaxial or biaxial compression(-tension) and the problem of cylinders under axial compression.

Let $(e_1, e_2, e_3)$ be an orthonormal basis and the first Kirchhoff stress tensor $\Pi$ be expressed in this basis as follows:

$$ \Pi = -\lambda \mathbf{e}_1 \otimes \mathbf{e}_1 - \eta \lambda \mathbf{e}_2 \otimes \mathbf{e}_2 = \begin{bmatrix} -\lambda & 0 & 0 \\ 0 & -\eta\lambda & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad (\lambda > 0, \eta \in \mathbb{R}) \qquad (19) $$

where the tensor expression and its matrix representation in the basis $(e_1, e_2, e_3)$ have been put together for brevity. The scalar $\eta$ represents the ratio between the load in the $e_2$-direction and that in the $e_1$-direction. The load is compressive in the $e_1$-direction, it is a tension or compression load in the $e_2$-direction, depending on the sign of $\eta$. In order to get explicit expressions for the tangent elastoplastic tensors, let us make the following assumption, which is satisfied in practice:

**Assumption 6.** The pre-critical deformations are small:

$$ \|\nabla \mathbf{U}\| \ll 1 \qquad (20) $$

Thus, the stress tensor $\Sigma$ writes

$$ \Sigma = \mathbf{F}^{-1} \cdot \Pi \approx \Pi \qquad (21) $$

The tensor $\mathbf{N}$ in Eq. (11) simplifies as follows:

$$ \mathbf{N} = \frac{\mu}{\sqrt{(H + 3\mu)(1 - \eta + \eta^2)}} [\mathbf{I} - 3\mathbf{e}_1 \otimes \mathbf{e}_1 + \eta(\mathbf{I} - 3\mathbf{e}_2 \otimes \mathbf{e}_2)] \\ = \frac{\mu}{\sqrt{(H + 3\mu)(1 - \eta + \eta^2)}} \begin{bmatrix} -2 + \eta & 0 & 0 \\ 0 & 1 - 2\eta & 0 \\ 0 & 0 & 1 + \eta \end{bmatrix} \qquad (22) $$

Hence, the material tangent elastoplastic tensor in Eq. (10) reads

$$ \mathbf{D}^p = \mathbf{D} - \frac{\mu^2}{(H + 3\mu)(1 - \eta + \eta^2)} [\mathbf{I} - 3\mathbf{e}_1 \otimes \mathbf{e}_1 + \eta(\mathbf{I} - 3\mathbf{e}_2 \otimes \mathbf{e}_2)] \\ \otimes [\mathbf{I} - 3\mathbf{e}_1 \otimes \mathbf{e}_1 + \eta(\mathbf{I} - 3\mathbf{e}_2 \otimes \mathbf{e}_2)] \qquad (23) $$

The components of $\mathbf{D}^p$ in the orthonormal basis $(e_1, e_2, e_3)$ are

$$ D_{1111}^p = \bar{\lambda} + 2\mu - \frac{\mu^2 (2 - \eta)^2}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{2222}^p = \bar{\lambda} + 2\mu - \frac{\mu^2 (1 - 2\eta)^2}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{3333}^p = \bar{\lambda} + 2\mu - \frac{\mu^2 (1 + \eta)^2}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{1122}^p = \bar{\lambda} + \frac{\mu^2 (2 - \eta)(1 - 2\eta)}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{1133}^p = \bar{\lambda} + \frac{\mu^2 (2 - \eta)(1 + \eta)}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{2233}^p = \bar{\lambda} - \frac{\mu^2 (1 - 2\eta)(1 + \eta)}{H + 3\mu (1 - \eta + \eta^2)} $$

$$ D_{1212}^p = D_{1313}^p = D_{2323}^p = \mu \qquad (24) $$

The other components are either zero or derived from Eqs. (24) using both major and minor symmetries of tensor $D^p$ ($D_{ijkl}^p = D_{klij}^p = D_{jikl}^p = D_{ijkl}^p$).

The nominal tangent elastoplastic tensor in Eq. (13) becomes

$$ \mathbf{K}^p \approx \frac{\partial \Sigma}{\partial \mathbf{E}} + (\mathbf{I} \cdot \mathbf{\Sigma})^T = \mathbf{D}^p - \lambda \mathbf{e}_i \otimes \mathbf{e}_1 \otimes \mathbf{e}_1 \otimes \mathbf{e}_i - \eta \lambda \mathbf{e}_i \otimes \mathbf{e}_2 \otimes \mathbf{e}_2 \otimes \mathbf{e}_i \qquad (25) $$

which is independent of the spatial coordinates. Furthermore, when dealing with two-dimensional models like plates or shells, an ad hoc assumption has to be added in order to enforce some specific stress state in the body. Namely, the normal material stress is assumed to be zero: $\Sigma_{33} = 0$ (although this contradicts the equality $E_{33} = 0$ induced by the classical kinematics). Taking into account this assumption leads one to recast the 3D constitutive law and to replace tensor $D^p$ in (25) with a suitable one denoted by $C^p$, which will be computed in the next section.

Eventually, let us make the following assumption which is satisfied by the plates and cylinders considered in this work:

**Assumption 7.** The whole solid is plastified on the fundamental branch:

$$ \Omega_f^p(t_c) = \Omega_0 \qquad (26) $$

which implies that $\mathbf{K}^p = \mathbf{K}^p$ throughout the solid.

Then, the bifurcation equation (16) writes in the biaxial stress case:

$$ \forall \delta \mathbf{U}, \quad \int_{\Omega_0} \nabla^T \delta \mathbf{U} : (\mathbf{C}^p - \lambda \mathbf{e}_i \otimes \mathbf{e}_1 \otimes \mathbf{e}_1 \otimes \mathbf{e}_i - \eta \lambda \mathbf{e}_i \otimes \mathbf{e}_2 \otimes \mathbf{e}_2 \otimes \mathbf{e}_i) : \nabla \mathbf{X} d\Omega_0 = 0 \qquad (27) $$

Furthermore, from Eq. (14) and assumption (20), one has $\mathbf{M} = \mathbf{N} \cdot \mathbf{F}^T \approx \mathbf{N}$. Hence:

$$ \mathbf{M} : \nabla \mathbf{X} = \frac{\mu}{\sqrt{(H + 3\mu)(1 - \eta + \eta^2)}} [\mathbf{I} - 3\mathbf{e}_1 \otimes \mathbf{e}_1 + \eta(\mathbf{I} - 3\mathbf{e}_2 \otimes \mathbf{e}_2)] : \nabla \mathbf{X} = \frac{\mu}{\sqrt{(H + 3\mu)(1 - \eta + \eta^2)}} (\mathrm{div} \, \mathbf{X} - 3(\nabla \mathbf{X})_{11} + \eta(\mathrm{div} \, \mathbf{X} - 3(\nabla \mathbf{X})_{22})) \qquad (28) $$

Accordingly, the initial slope of the bifurcated branch (18) becomes

$$ \lambda_1 = \max_{\mathbf{X} \in \Omega_0} \left( \frac{3(\nabla \mathbf{X})_{11} - \mathrm{div} \, \mathbf{X} + \eta(3(\nabla \mathbf{X})_{22} - \mathrm{div} \, \mathbf{X})}{3(\nabla \mathbf{U}_{f,\lambda}^c)_{11} - \mathrm{div} \, (\mathbf{U}_{f,\lambda}^c) + \eta(3(\nabla \mathbf{U}_{f,\lambda}^c)_{22} - \mathrm{div} \, (\mathbf{U}_{f,\lambda}^c))} \right) \qquad (29) $$

where $\mathbf{U}_{f,\lambda}^c$ denotes the differentiation $\mathbf{U}_{f,\lambda}(X, \lambda)$ evaluated at $\lambda = \lambda_c$.

Eventually, let us check that the transversality condition (17) is satisfied in the biaxial stress state. From expression (25) for the nominal tangent elastoplastic tensor and assumption (26), it can be readily shown that the transversality coefficient $\mathrm{Tr}$ in (17) is

$$ \mathrm{Tr} = \int_{\Omega_0} (|\mathbf{X}_1|^2 + \eta |\mathbf{X}_2|^2) d\Omega_0 \qquad (30) $$

It will be seen in the next sections that the eigenmode $\mathbf{X}$ does depend on abscissa $X_1$, so that $\mathrm{Tr}$ is strictly positive as long as $\eta$ is not too negative and the transversality condition (17) is always satisfied.

**Remark 1.** (i) All the previous developments apply to the special case of a uniaxial stress state by taking $\eta = 0$.

(ii) Whereas the above results have been obtained using the $J_2$ flow theory of plasticity, similar expressions can straightforwardly be derived for the deformation theory. The same methodology is valid, provided the elastoplastic tangent moduli in Eq. (24) are replaced with those corresponding to the deformation theory.
## 5. Plastic bifurcation of a rectangular plate under uniaxial or biaxial compression(-tension)

Let us consider a rectangular plate with reference dimensions a along **x** = **e**₁, b along **y** = **e**₂ (a ≥ b or a < b) and thickness e along **z** = **e**₃ (e ≪ a, b). The edges are subjected to a compressive nominal stress Πₓₓ = -λ < 0 normal to the sides of length b and another stress Πᵧᵧ = -ηλ normal to the sides of length a (which is compressive when η > 0). The biaxial stress state can be dealt with by applying the results established in Section 4.

The plate is made of an elastoplastic material with a linear isotropic hardening (make A = Hα and H = constant in (5) and (8)). Use will also be made of the tangent modulus Eₜ related to the Young's modulus E and the isotropic hardening modulus H by 1/Eₜ = 1/E + 1/H. The yield stress σ₀ is assumed small enough for the plastic strains to appear before buckling.

The kinematics is described by the Love-Kirchhoff model, hence the displacement field is computed from the three displacement components U(X, Y, t), V(X, Y, t) and W(X, Y, t) of the middle surface:

$$ \mathbf{U} = \begin{cases} U - ZW_{,X} \\ V - ZW_{,Y} \\ W \end{cases} \qquad (31) $$

### 5.1. Critical load and bifurcation mode

The bifurcation mode **X** and the displacement variation δ**U** are expressed according to the Love-Kirchhoff kinematics:

$$ \mathbf{X} = \begin{cases} \mathcal{U} - Z\mathcal{W}_{,X}, \\ \mathcal{V} - Z\mathcal{W}_{,Y}, \\ \mathcal{W} \end{cases}, \quad \delta\mathbf{U} = \begin{cases} \delta U - Z\delta W_{,X} \\ \delta V - Z\delta W_{,Y} \\ \delta W \end{cases} \qquad (32) $$

Whence the bifurcation mode gradient:

$$ \nabla \mathbf{X} = \begin{bmatrix} \mathcal{U}_{,X} - Z\mathcal{W}_{,XX} & \mathcal{U}_{,Y} - Z\mathcal{W}_{,XY} & -\mathcal{W}_{,X} \\ \mathcal{V}_{,X} - Z\mathcal{W}_{,YY} & \mathcal{V}_{,Y} - Z\mathcal{W}_{,YY} & -\mathcal{W}_{,Y} \\ \mathcal{W}_{,X} & \mathcal{W}_{,Y} & 0 \end{bmatrix} \qquad (33) $$

In the case of a thin plate, the normal stress Σ<sub>zz</sub> is assumed to be zero (the subscripts Z and 3 are used alike). Then, the reduced tensor **C**<sup>p</sup> involved in (27) can be recast as

$$ C_{ijkl}^{p} = D_{ijkl}^{p} - \frac{D_{ij33}^{p} D_{33kl}^{p}}{D_{3333}^{p}} \quad (i,j) \neq (3,3), \quad (k,l) \neq (3,3) \qquad (34) $$

The tensor **C**<sup>p</sup> has the major and both minor symmetries. After integrating the bifurcation equation (27) through the thickness, integrating by parts with respect to X and Y, and eliminating negligible higher-order terms, one obtains the following local equation for the out-of-plane component **W** of the eigenmode:

$$ e^2[\alpha\mathcal{W}_{,XXXX} + \gamma\mathcal{W}_{,YYYY} + (2\beta + 4\mu)\mathcal{W}_{,XYY}] + 12\lambda_c(\mathcal{W}_{,XX} + \eta\mathcal{W}_{,YY}) = 0 \qquad (35) $$

where the notations (α, β, γ) defined as follows have been used for brevity:

$$ \alpha = C_{1111}^{p} = E \frac{1 + 4\eta(\eta - 1) + 3\frac{E_T}{E}}{(5 - 4\nu)(1 + \eta^2) - 2(4 - 5\nu)\eta - (1 - 2\nu)[(1 - 2\nu)(1 + \eta^2) - 2(2 - \nu)\eta]\frac{E_T}{E}} $$

$$ \begin{aligned} \beta &= C_{1122}^{p} \\ &= E \frac{2 - 5\eta + 2\eta^2 - [2(1 - 2\nu)(1 + \eta^2) - (5 - 4\nu)\eta]\frac{E_T}{E}}{(5 - 4\nu)(1 + \eta^2) - 2(4 - 5\nu)\eta - (1 - 2\nu)[(1 - 2\nu)(1 + \eta^2) - 2(2 - \nu)\eta]\frac{E_T}{E}} \\ \gamma &= C_{2222}^{p} \\ &= E \frac{(2 - \eta)^2 + 3\eta^2\frac{E_T}{E}}{(5 - 4\nu)(1 + \eta^2) - 2(4 - 5\nu)\eta - (1 - 2\nu)[(1 - 2\nu)(1 + \eta^2) - 2(2 - \nu)\eta]\frac{E_T}{E}} \end{aligned} \qquad (36) $$

Eq. (35) can be analytically solved in the classical case where the four edges are simply supported on the Oxy plane. The mode is assumed to be sinusoidal along the **x** and **y**-axes in order to automatically satisfy the boundary conditions:

$$ \left\{ \begin{array}{l} \mathcal{U} = 0 \\ \mathcal{V} = 0 \\ \mathcal{W} = \frac{e}{2} \sin \frac{p\pi X}{a} \sin \frac{q\pi Y}{b} \quad (p, q \in \mathbb{N}^*) \end{array} \right. \qquad (37) $$

By substituting for **W** from Eq. (37) into Eq. (35), one gets the expression for the critical stress:

$$ \lambda_c = \frac{\pi^2 e^2}{12} \frac{\alpha \frac{p^4}{a^4} + \gamma \frac{q^4}{b^4} + (2\beta + 4\mu) \frac{p^2 q^2}{a^2 b^2}}{\frac{p^2}{a^2} + \eta \frac{q^2}{b^2}} \qquad (38) $$

The integers p and q minimizing λ<sub>c</sub> are obtained after solving the following equations:

$$ \left\{ \begin{array}{l} \alpha \frac{p^4}{a^4} + 2\alpha\eta \frac{p^2 q^2}{a^2 b^2} + ((2\beta + 4\mu)\eta - \gamma) \frac{q^4}{b^4} = 0 \\ (2\beta + 4\mu - \alpha\eta) \frac{p^4}{a^4} + 2\gamma \frac{p^2 q^2}{a^2 b^2} + \gamma\eta \frac{q^4}{b^4} = 0 \end{array} \right. \qquad (39) $$

Since the expressions for p and q are intricate in the case of arbitrary a/b, let us confine ourselves to the case of a square plate (a = b) in the equi-biaxial stress state (η = 1). The problem is symmetric with respect to **x** and **y**, so that integers p and q in (37) must be equal and Eq. (38) becomes

$$ \lambda_c = \frac{\pi^2 e^2}{24a^2}(\alpha + 2\beta + \gamma + 4\mu)p^2 \qquad (40) $$

Clearly, the integer p minimizing λ<sub>c</sub> is p = 1. Relations (36) and (40) then yield the critical stress:

$$ \lambda_c = E \frac{\pi^2 e^2}{12a^2} \frac{1 + \frac{3E_T}{E}}{(1 + \nu) \left[ 1 + (1 - 2\nu) \frac{E_T}{E} \right]} \qquad (41) $$

In the case of a square plate, with η ∈ [-0.1, 1], the sought mode is obtained by taking p = q = 1. Table 1 displays the numerical values which are used throughout. Fig. 1 plots the plastic buckling stress of a square plate under biaxial compression(-tension) with simply supported edges versus the load ratio η, including the particular case of a uniaxial loading (η = 0) which will be studied in details in the sequel. As indicated in the figure, adding a small tensile stress in the **y**-direction has a stabilizing effect on the behavior of the plate under uniaxial compression. On the contrary,

<table><caption>**Table 1**

Geometric and material parameters.</caption><thead><tr><th>Dimension b</th><th>Thickness e</th><th>Young's modulus E</th><th>Poisson's ratio v</th><th>Tangent modulus E<sub>T</sub></th></tr></thead><tbody><tr><td>100 mm</td><td>3 mm</td><td>70 000 MPa</td><td>0.3</td><td>30 000 MPa</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Plastic buckling stress for a square plate under biaxial compression(-tension) with simply supported edges.</figcaption></figure>

adding a compressive stress in the y-direction lowers the critical load. The minimum value for the critical load corresponds to the equi-biaxial compression case ($\eta = 1$) and is given by Eq. (41).

### 5.1.1. Particular case of a uniaxial loading

The particular case of a uniaxial stress state is readily obtained from the biaxial case by taking $\eta = 0$. The entailed simplifications will enable one to discuss further the critical value and consider another boundary condition type.

#### 5.1.1.1. Simply supported edge solution.
The critical stress is derived from Eq. (38):

$$ \lambda_c = \frac{\pi^2 e^2}{12} \left[ \alpha \frac{p^2}{a^2} + \gamma \frac{q^4 a^2}{p^2 b^4} + (2\beta + 4\mu) \frac{q^2}{b^2} \right] \quad (42) $$

It can readily be checked that the wave numbers corresponding to the smallest critical load are

$$ \left\{ \begin{array}{l} p = \text{Int}\left(\sqrt[4]{\frac{\gamma}{\alpha}}\frac{a}{b}\right) = \text{Int}\left(\sqrt[4]{\frac{4}{1 + 3\frac{E_T}{E}\frac{b}{E}}}\right) \\ q = 1 \end{array} \right. \quad (43) $$

where $Int(r)$ with $r \in \mathbb{R}$ means the one of the two integers closest to $r$ and minimizing Eq. (42) with $q = 1$.

With the so-defined integer $p$, the critical stress reads

$$ \lambda_c = E \frac{\pi^2 e^2}{12b^2} \left[ \frac{2}{1+v} + \frac{4 \left[ 1 - (1-2v) \frac{E_T}{E} \right] + \left( 1 + 3 \frac{E_T}{E} \right) \left( \frac{pb}{a} \right)^2 + 4 \left( \frac{a}{pb} \right)^2}{5 - 4v - (1-2v)^2 \frac{E_T}{E}} \right] \quad (44) $$

In the particular case of a square plate ($a = b$), since $E_T/E$ lies between 0 and 1 in practice, Eq. (43) leads to choosing between $p = 1$ and $p = 2$, i.e. choosing the greater number of $\alpha + \gamma$ and $4\alpha + \gamma/4$. It can readily be checked that for all $E_T/E \in [0, 1]$, $\alpha + \gamma < 4\alpha + \gamma/4$, so that $p = 1$ is the proper value. One eventually finds the critical stress obtained in [10] by a different way:

$$ \lambda_c = E \frac{\pi^2 e^2}{12a^2} \left[ \frac{2}{1+v} + \frac{9 + (8v - 1) \frac{E_T}{E}}{5 - 4v - (1-2v)^2 \frac{E_T}{E}} \right] \quad (45) $$

#### 5.1.1.2. Simply supported and built-in edge solution.
Let us consider now the little-known problem of a rectangular plate with simply supported and built-in edges. The plate is the same as before, except that the simply supported edge at $X = 0$ is now replaced with a clamped edge. The edges $X = 0$ and $X = a$ remain subjected to a compressive nominal stress $\Pi_{XX} = -\lambda < 0$ along the x-axis.

The mode is still assumed to be sinusoidal along the y-axis in order to automatically satisfy the boundary conditions on edges $Y = 0$ and $Y = b$:

$$ \mathcal{W} = \mathcal{F}(X) \sin \frac{q\pi Y}{b} \quad (q \in \mathbb{N}^*) \quad (46) $$

where function $\mathcal{F}$ acting as a magnitude varying with abscissa $X$ is to be determined. By substituting for $\mathcal{W}$ from Eq. (46) into Eq. (35), it leads to the following differential equation in $\mathcal{F}$:

$$ e^2 \left[ \alpha \mathcal{F}_{,XXXX} - (2\beta + 4\mu) \left(\frac{q\pi}{b}\right)^2 \mathcal{F}_{,XX} + \gamma \left(\frac{q\pi}{b}\right)^4 \mathcal{F} \right] + 12 \lambda_c \mathcal{F}_{,XX} = 0 \quad (47) $$

Taking into account the order of magnitude of the geometric and material parameters, one finds that the general solution of Eq. (47) can be written as

$$ \mathcal{F}(X) = k_1 \cos(rX) + k_2 \sin(rX) + k_3 \cos(sX) + k_4 \sin(sX) \quad (48) $$

with

$$ \left\{ \begin{array}{l} r = \sqrt{\frac{12\lambda - (2\beta + 4\mu)\left(\frac{eq\pi}{b}\right)^2 - \sqrt{\left[12\lambda - (2\beta + 4\mu)\left(\frac{eq\pi}{b}\right)^2\right]^2 - 4\alpha\gamma\left(\frac{eq\pi}{b}\right)^4}}{2e^{2\alpha}}} \\ s = \sqrt{\frac{12\lambda - (2\beta + 4\mu)\left(\frac{eq\pi}{b}\right)^2 + \sqrt{\left[12\lambda - (2\beta + 4\mu)\left(\frac{eq\pi}{b}\right)^2\right]^2 - 4\alpha\gamma\left(\frac{eq\pi}{b}\right)^4}}{2e^{2\alpha}}} \end{array} \right. \quad (49) $$

The boundary conditions at the clamped edge $X = 0$ ($\mathcal{W} = 0$ and $\mathcal{W}_{,X} = 0$) enable one to eliminate $k_3$ and $k_4$ according to

$$ \left\{ \begin{array}{l} k_3 = -k_1 \\ k_4 = -k_2 \frac{r}{s} \end{array} \right. \quad (50) $$
With expression (46), it is not a trivial matter to satisfy the remaining boundary conditions at X = a:

$$ \begin{cases} \mathcal{W} = 0 \\ \alpha \mathcal{W}^{\prime},_{XX} + \beta \mathcal{W}^{\prime},_{YY} = 0 \end{cases} \qquad (51) $$

These boundary conditions lead to the following linear equation system:

$$ \begin{cases} k_1(\cos(\mathrm{ra}) - \cos(\mathrm{sa})) + k_2(\sin(\mathrm{ra}) - \frac{r}{s}\sin(\mathrm{sa})) = 0 \\ k_1(s^2 \cos(\mathrm{sa}) - r^2 \cos(\mathrm{ra})) + k_2(rs \sin(\mathrm{sa}) - r^2 \sin(\mathrm{ra})) = 0 \end{cases} \qquad (52) $$

The critical loads $\lambda_c$ are such that the determinant of the previous system vanishes, in order to ensure a non-zero solution for $(k_1, k_2)$. Hence the equation for the critical loads $\lambda_c$:

$$ r(s^2 - r^2) \cos(\mathrm{ra}) \sin(\mathrm{sa}) + s(r^2 - s^2) \cos(\mathrm{sa}) \sin(\mathrm{ra}) = 0 \qquad (53) $$

The minimum critical load is then obtained by numerically solving Eq. (53). A single half-wave mode along the y-axis is obtained ($q = 1$) when minimizing the critical load with respect to parameter $q$. The critical load for a large range of the aspect ratio $a/b$ is depicted in Fig. 2 together with the one given by Eq. (44) corresponding to the previous case of simply supported edges, for comparison purposes.

As expected, the two boundary condition types considered yield the same critical load for large $a/b$ values. On the other hand, the critical load $\lambda_c$ is identical for the two boundary condition types, when the aspect ratio $a/b$ takes some particular values, e.g. $a/b \approx 1.25, 2.2, 3$, as shown in Fig. 2. This fact is noteworthy, although the authors have no explanations for it.

## 5.2. Initial slope of the bifurcated branch

The initial slope $\lambda_1$ of the bifurcated branch is given by Eq. (18) and depends on the fundamental solution. The fundamental solution for the plate under biaxial loading is obtained by solving the elastoplasticity equation set, Eqs. (3)–(7). Without giving

details, one finds the fundamental displacement field:

$$ \mathbf{U}_f = \left[ \begin{array}{l} -\frac{\lambda_c}{E}(1 - \eta v) + \frac{1}{2H}\left(\lambda - \frac{\sigma_0}{\sqrt{1 - \eta + \eta^2}}\right)(-2 + \eta) \mathbf{X}\mathbf{x} \\ + \left[\frac{\lambda_c}{E}(\eta - v) + \frac{1}{2H}\lambda - \frac{\sigma_0}{\sqrt{1 - \eta + \eta^2}}\right)(-2\eta + 1) \mathbf{Y}\mathbf{y} \end{array} \right] \qquad (54) $$

In the case of simply supported edges, substituting Eqs. (37) and (54) into Eq. (29) yields the initial slope of the bifurcated branch:

$$ \lambda_1 = E_T \frac{\pi^2 e^2}{2} \frac{\frac{p^2}{a^2}(2 - \eta) + \frac{q^2}{b^2}(2\eta - 1)}{5 - 8\eta + 5\eta^2 - (1 - 2v)(1 - 4\eta + \eta^2)\frac{E_T}{E}} \qquad (55) $$

where $p$ and $q$ are the integers minimizing $\lambda_c$ in (38), as explained in (39). In the case of a square plate ($a = b$) with $\eta \in [-0.1, 1]$, Eq. (55) yields

$$ \lambda_1 = E_T \frac{\pi^2 e^2}{2a^2} \frac{1 + \eta}{5 - 8\eta + 5\eta^2 - (1 - 2v)(1 - 4\eta + \eta^2)\frac{E_T}{E}} \qquad (56) $$

The initial slope $\lambda_1$ is plotted versus the load ratio $\eta$ in Fig. 3, which shows that its maximum value is reached at $\eta \approx 0.88$.

In the case of a square plate in the equi-biaxial stress state, the initial slope is

$$ \lambda_1 = E_T \frac{\pi^2 e^2}{2a^2} \frac{1}{1 + (1 - 2v)\frac{E_T}{E}} \qquad (57) $$

In the case of a uniaxial stress state, the initial slope of the bifurcated branch is

$$ \lambda_1 = E_T \frac{\pi^2 e^2 \left( \frac{p^2}{a^2} - \frac{1}{2b^2} \right)}{5 - (1 - 2v) \frac{E_T}{E}} \qquad (58) $$

where $p$ is given by (43).

# 6. Plastic bifurcation of a cylinder under axial compression

The last section of this paper is devoted to thin cylinders under axial compression. Let us consider a cylindrical tube defined in the reference configuration by the length $l$ along the x-axis, the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Plastic buckling stress for a rectangular plate under uniaxial compression.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Initial slope of the bifurcated branch for a square plate under biaxial compression(-tension) with simply supported edges.</figcaption></figure>

average radius R and the thickness e ≪ R. One end of the cylinder lies in the plane X = 0 where the radial displacements and the rotations are free; whereas the other end at X = l is loaded by a uniformly distributed nominal axial compressive stress Π$_{XX}$ = −λ < 0. The material properties are the same as those of the plate in the previous section. To simplify the discussion, we still assume that the yield stress σ$_{0}$ is chosen small enough so that the plastification occurs before the tangent modulus critical load is reached.

The kinematics is described by the Koiter model, equivalent to the Love–Kirchhoff plate model. The displacement field of the middle surface is defined by the three components: U(X, θ, t) (axial), V(X, θ, t) (circumferential) and W(X, θ, t) (radial). Consider a current particle in the cylinder, whose position in the reference configuration is defined by (X, θ) and the thickness parameter ρ (−e/2 ≪ ρ ≪ e/2). Its displacement can be expressed in the local cylindrical basis (x, e$_{θ}$, e$_{r}$) as:

$$
\mathbf{U} = \begin{bmatrix} U - \rho W_x \\ V - \frac{\rho}{R}(W_{,\theta} - V) \\ W \end{bmatrix} \tag{59}
$$

6.1. Critical load and bifurcation mode

The bifurcation mode **X** and the displacement variation δ**U** are
expressed according to the kinematics (59):

$$
\mathbf{X} = \begin{bmatrix} \mathcal{U} - \rho \mathcal{W}_{,x} \\ \mathcal{V} - \frac{\rho}{R}(\mathcal{W}_{,\theta} - \mathcal{V}) \\ \mathcal{W} \end{bmatrix}, \quad \delta \mathbf{U} = \begin{bmatrix} \delta U - \rho \delta W_{,x} \\ \delta V - \frac{\rho}{R}(\delta W_{,\theta} - \delta V) \\ \delta W \end{bmatrix} \quad (60)
$$

Whence the bifurcation mode gradient:

$$
\nabla \mathbf{X} = \begin{bmatrix} \mathcal{U}_{,x} - \rho \mathcal{W}_{,xx} & \frac{1}{R}(\mathcal{U}_{,\theta} - \rho \mathcal{W}_{,x\theta}) & -\mathcal{W}_{,x} \\ \mathcal{V}_{,x} - \frac{\rho}{R}(\mathcal{W}_{,\theta x} - \mathcal{V}_{,x}) & \frac{1}{R}(\mathcal{V}_{,\theta} - \frac{\rho}{R}(\mathcal{W}_{,\theta\theta} - \mathcal{V}_{,\theta}) + \mathcal{W}_{,\theta}) & \frac{1}{R}(\mathcal{V}_{,\theta} - \mathcal{W}_{,\theta}) \\ \mathcal{W}_{,x} & \frac{1}{R}(\mathcal{W}_{,\theta} - \mathcal{V}_{,x} + \frac{\rho}{R}(\mathcal{W}_{,\theta} - \mathcal{V}_{,x})) & 0 \end{bmatrix} \quad (61)
$$

As in the case of plates, the normal stress Σrr is assumed to be
zero and the reduced tensor C^p in Eq. (34) can be used again. One
has to integrate the bifurcation equation (27) in the shell
thickness and then adequately integrate by parts with respect to

X and θ. Without going into details, we give the resulting local
equations:

$$
\left\{
\begin{aligned}
& \left(\alpha - \lambda_c\right) \mathcal{U}_{,XX} + \frac{\mu}{R^2} \mathcal{U}_{,00} + \frac{\beta + \mu}{R} \mathcal{V}_{,X\theta} + \frac{\beta}{R} \mathcal{W}_{,X} = 0 \\
& + \frac{\beta + \mu}{R} \mathcal{U}_{,X\theta} - \mu \frac{e^2}{12R^2} \mathcal{V} + (\mu - \lambda_c) \left(1 + \frac{e^2}{12R^2}\right) \mathcal{V}_{,XX} \\
& + \frac{\gamma}{R^2} \left(1 + \frac{e^2}{12R^2}\right) \mathcal{V}_{,00} \\
& + \mu \frac{e^2}{12R^2} \mathcal{W}_{,\theta} - (\beta + 2\mu - \lambda_c) \frac{e^2}{12R^2} \mathcal{W}_{,XX\theta} - \gamma \frac{e^2}{12R^4} \mathcal{W}_{,000} = 0 \\
& \frac{\beta}{R} \mathcal{U}_{,X} + \left(\frac{\gamma}{R^2} + \mu \frac{e^2}{12R^4}\right) \mathcal{V}_{,\theta} - (\beta + 2\mu - \lambda_c) \frac{e^2}{12R^2} \mathcal{V}_{,XX\theta} \\
& - \gamma \frac{e^2}{12R^4} \mathcal{V}_{,000} + \frac{\gamma}{R^2} \mathcal{W} + \lambda_c \mathcal{W}_{,XX} - \mu \frac{e^2}{12R^4} \mathcal{W}_{,00} \\
& + (\alpha - \lambda_c) \frac{e^2}{12} \mathcal{W}_{,XXXX} + (2\beta + 4\mu - \lambda_c) \frac{e^2}{12R^2} \mathcal{W}_{,XX00} \\
& + \gamma \frac{e^2}{12R^4} \mathcal{W}_{,0000} = 0
\end{aligned}
\right.
\tag{62}
$$

where coefficients α, β and γ are defined as in (36) by making
η = 0.

In the axisymmetric case, the circumferential displacements V and
**V** are zero and all fields are independent of the circumfer-
ential parameter θ. Only two equations remain among Eqs. (62),
which after omitting negligible terms can be written as

$$
\left\{
\begin{aligned}
& \alpha \mathcal{U}_{,X} + \frac{\beta}{R} \mathcal{W} = 0 \\
& \alpha \frac{e^2}{12} \mathcal{W}_{,XXXX} + \lambda_c \mathcal{W}_{,XX} + \frac{\gamma}{R^2} \mathcal{W} + \frac{\beta}{R} \mathcal{W}_{,X} = 0
\end{aligned}
\right.
\tag{63}
$$

Combining these two relations leads to a single equation
involving the radial component of the eigenmode:

$$
\alpha \frac{e^2}{12} \mathcal{W}_{,XXXX} + \lambda_c \mathcal{W}_{,XX} + \frac{1}{R^2} \left( \gamma - \frac{\beta^2}{\alpha} \right) \mathcal{W} = 0 \quad (64)
$$

Expression (64), valid for plastic buckling, is similar to the
simplified form (in the axisymmetric case) of the Donnell
equation [30] for the well-known elastic buckling of a cylinder
under axial compression (see also [31]).
### 6.1.1. Free edge solution

Since it is difficult to determine the bifurcation mode in the case of a free edge cylinder with finite length, we limit consideration to a cylinder of infinite length. The symmetry makes it possible to study one half of the cylinder only, with the free edge at $X = 0$ and the plane of symmetry put at infinity.

The axisymmetric buckling mode is assumed to be a sinusoidal function weighted by a decreasing exponential function:

$$ \mathcal{W}' = (k_1 \cos(\xi X) + k_2 \sin(\xi X)) \exp(-\xi X) \quad (65) $$

where $k_1, k_2, \xi$ and $\zeta$ are to be determined. Relation (65) expresses the fact that the radial displacement is significant at the end of the cylinder and is zero in the plane of symmetry.

By inserting Eq. (65) in (64), one obtains

$$ \left\{ \begin{aligned} & \alpha \frac{e^2}{12} (k_1 \xi^4 + 4k_2 \xi^3 \zeta - 6k_1 \xi^2 \zeta^2 - 4k_2 \xi \zeta^3 + k_1 \zeta^4) \\ & + \lambda_c (-k_1 \zeta^2 - 2k_2 \xi \zeta + k_1 \zeta^2) + \left(\gamma - \frac{\beta^2}{\alpha}\right) \frac{k_1}{R^2} = 0 \\ & \alpha \frac{e^2}{12} (k_2 \xi^4 - 4k_1 \xi^3 \zeta - 6k_2 \xi^2 \zeta^2 + 4k_1 \xi \zeta^3 + k_2 \zeta^4) \\ & + \lambda_c (-k_2 \xi^2 + 2k_1 \xi \zeta + k_2 \zeta^2) + \left(\gamma - \frac{\beta^2}{\alpha}\right) \frac{k_2}{R^2} = 0 \end{aligned} \right. \quad (66) $$

The boundary conditions at infinity—namely the displacement $\mathcal{W}'$ and its derivatives with respect to $X$ are zero—result in the single condition $\zeta > 0$. On the edge $X = 0$, the boundary conditions are less trivial:

$$ \left\{ \begin{aligned} \mathcal{W}'_{,XX} &= 0 \\ \alpha \frac{e^2}{12} \mathcal{W}'_{,XXX} + \lambda_c \mathcal{W}'_{,X} &= 0 \end{aligned} \right. \quad (67) $$

that is to say

$$ \left\{ \begin{aligned} & -k_1 \xi^2 - 2k_2 \xi \zeta + k_1 \zeta^2 = 0 \\ & \alpha \frac{e^2}{12} (-k_2 \xi^3 + 3k_1 \xi^2 \zeta + 3k_2 \xi \zeta^2 - k_1 \zeta^3) + \lambda_c (k_2 \xi - k_1 \zeta) = 0 \end{aligned} \right. \quad (68) $$

The bifurcation mode is defined within a factor which can be determined by enforcing $k_2 = -e$. The four Eqs. (66)–(68) then give

$$ \begin{aligned} k_1 &= e\sqrt{3} \\ \xi &= \frac{p}{\sqrt{Re}}, \quad p \doteq \sqrt[4]{\frac{27\left[5 - 4\nu - (1 - 2\nu)^2 \frac{E_T}{E}\right] \frac{E_T}{E}}{\left(1 + 3\frac{E_T}{E}\right)^2}} \\ \zeta &= \frac{p}{\sqrt{3Re}} \end{aligned} \quad (69) $$

and the critical stress:

$$ \lambda_c^{\text{free edges}} = \frac{Ee}{R} \sqrt{\frac{\frac{E_T}{E}}{3\left[5 - 4\nu - (1 - 2\nu)^2 \frac{E_T}{E}\right]}} \quad (70) $$

This relation compares very well with Batterman's ones [13], which show that the critical stress is a little bit smaller for a finite cylinder and it hardly depends on the cylinder length.

The complete eigenmode is obtained by solving Eqs. (63) with the boundary conditions $\mathcal{W}'(0) = 0$ and (67):

$$ \left\{ \begin{aligned} \mathcal{W}' &= -e \sqrt{\frac{3e q}{R}} \frac{\sin\left(\frac{pX}{\sqrt{Re}}\right)}{\sin\left(\frac{pX}{\sqrt{3Re}}\right)} \exp\left(-\frac{pX}{\sqrt{3Re}}\right) \\ \mathcal{W}' &= 0 \\ \mathcal{W}'' &= e \left[ \sqrt{3} \cos\left(\frac{pX}{\sqrt{Re}}\right) - \sin\left(\frac{pX}{\sqrt{Re}}\right) \right] \exp\left(-\frac{pX}{\sqrt{3Re}}\right) \end{aligned} \right. \quad (71) $$

where

$$ q = 2 \frac{1 - (1 - 2\nu) \frac{E_T}{E}}{1 + 3 \frac{E_T}{E}} \quad (72) $$

In the above, the mode has been normalized so that $\mathcal{W}'(0) = e\sqrt{3}$.

### 6.1.2. Simply supported edge solution

Consider now the case of simply supported edges, by assuming that the radial displacements are fixed and the rotations free on both edges. This case is simpler than the free edge case and does not require the infinite length assumption. The same bifurcation equation (64) can be used in order to get the critical load and the bifurcation mode, providing that the free edges fundamental solution is used. This approximation is justified as long as the pre-critical deformations are small.

The bifurcation mode is assumed to be sinusoidal so as to identically satisfy the kinematic boundary conditions:

$$ \mathcal{W}' = e \sin \frac{n\pi X}{l} \quad (n \in \mathbb{N}^*) \quad (73) $$

By inserting Eq. (73) in Eq. (64), we obtain the critical stress as a function of integer $n$:

$$ \lambda_c(n) = \alpha n^2 \pi^2 \frac{e^2}{12l^2} + \left(\gamma - \frac{\beta^2}{\alpha}\right) \frac{1}{n^2 \pi^2 R^2} \quad (74) $$

The $n$ value minimizing $\lambda_c(n)$ is

$$ n = \frac{2lp}{\pi\sqrt{3Re}} \quad (75) $$

In general, the previous relation provides a real value for $n$. However, if the cylinder is long enough this real value is large enough to be a good approximation of the sought integer $n$. Inserting Eq. (75) in Eq. (74) then yields the critical value:

$$ \lambda_c^{\text{supported edges}} = \frac{2Ee}{R} \sqrt{\frac{\frac{E_T}{E}}{3\left[5 - 4\nu - (1 - 2\nu)^2 \frac{E_T}{E}\right]}} = 2\lambda_c^{\text{free edges}} \quad (76) $$

Again, this result is in good agreement with Batterman's solution [12] obtained for the same boundary conditions, yet in a different way.

The eigenmode is obtained by solving Eqs. (63):

$$ \left\{ \begin{aligned} \mathcal{W}' &= \frac{qle}{n\pi R} \left( \cos \frac{n\pi X}{l} - 1 \right) \\ \mathcal{W}' &= 0 \\ \mathcal{W}'' &= e \sin \frac{n\pi X}{l} \end{aligned} \right. \quad (77) $$

where $n$ is given by Eq. (75). The mode has been normalized so that $\mathcal{W}'(0) = e$.
## 6.2. Initial slope of the bifurcated branch

### 6.2.1. Free edge solution

The initial slope $\lambda_1$ is given by (29). Here, the symbol $\Omega_0$ stands for the reference volume of the cylinder and the fundamental solution is

$$ \mathbf{U}_f = \begin{vmatrix} U_f(X) \\ 0 \\ W_f \end{vmatrix} = \begin{vmatrix} \left[ \frac{\lambda}{E_T} + \sigma_0 \left( \frac{1}{E_T} - \frac{1}{E} \right) \right] X \\ 0 \\ \frac{R}{2} \left[ \lambda \left( \frac{1}{E_T} - \frac{1-2\nu}{E} \right) - \sigma_0 \left( \frac{1}{E_T} - \frac{1}{E} \right) \right] \end{vmatrix} \quad (78) $$

In the case of an infinite cylinder, Eq. (60) entails $\text{div} \mathbf{X} = (\mathcal{W}^-/r) + \mathcal{W}_{\times} - \rho \mathcal{W}_{\times\times}$, where $r = R + \rho$ is the radius at a current position. From Eq. (78), one has

$$ \lambda_1 = \max_{r \in [R-e/2, R+e/2]} \frac{\exp\left(-\frac{pX}{\sqrt{3Re}}\right) \left\{ \left(2q_R^e + \frac{e}{r}\right) \left[\sqrt{3} \cos\left(\frac{pX}{\sqrt{Re}}\right) - \sin\left(\frac{pX}{\sqrt{Re}}\right)\right] + \frac{16}{3} p^2 \rho^2 \sin\left(\frac{pX}{\sqrt{Re}}\right) \right\}}{\frac{1}{2E_T} \left[ 4 + \frac{R}{r} - \frac{R}{r} (1-2\nu) \frac{E_T}{E} \right]} \quad (79) $$

Since the cylinder is a thin shell, $1/r$ can be replaced by $1/R$ in Eq. (79). Hence:

$$ \lambda_1^{\text{free edges}} = \frac{2E_T e}{R} \frac{5 + (8\nu - 1)\frac{E_T}{E}}{\left(1 + 3\frac{E_T}{E}\right) \left[5 - (1 - 2\nu)\frac{E_T}{E}\right]} \max_{\tau \ge 0} [-f(\tau, \rho)] \quad (80) $$

where

$$ \begin{gathered} f(\tau, \rho) = \exp\left(-\frac{\tau}{\sqrt{3}}\right) \left[\sqrt{3} \cos \tau - (1 - C\rho) \sin \tau\right] \\ C = \frac{16\sqrt{3} \left[5 - 4\nu - (1 - 2\nu)^2 \frac{E_T}{E} \right] \frac{E_T}{E}}{5 + (8\nu - 1) \frac{E_T}{E}} > 0 \end{gathered} \quad (81) $$

and $\rho$ stands for $\rho/e$ for brevity.

Assume for definiteness that the ratio $E_T/E$ is small enough to have $C<4$. A thorough analysis of coefficient C defined by (81) shows that $E_T/E$ must be bounded by $E_T/E<5/49 \approx 0.10204$ for

the inequality $C<4$ to be fulfilled. One then finds

$$ \begin{aligned} \max(-f) = & -\min(f) = \frac{1}{4}\sqrt{3(16 + 4C + C^2)} \\ & \times \exp\left[-\frac{1}{\sqrt{3}}\tan^{-1}\left(-\frac{\sqrt{3}4 + C}{4 - C}\right)\right] \end{aligned} \quad (82) $$

The extremum value of $f$ is reached at $\rho = -\frac{1}{2}$ and the smaller $\tau > 0$ value satisfying $\tan \tau = -\frac{\sqrt{3}(4 + C)}{4 - C}$.

So far, the initial slope $\lambda_1$ has been computed using the eigenmode X given by Eqs. (71), such that the radial displacement $W(0)$ at $X=0$ is positive, i.e. the free edge is bent outwards. Now let us change the signs of the modal displacements $\mathcal{W}$ and $\mathcal{W}^-$ in (71), which implies that the radial displacement $W(0)$ at $X=0$ is negative and the free edge is bent inwards. In that case, the same calculation as above leads to another expression for the initial slope:

$$ \lambda_1^{\text{inward free edges}} = \frac{2E_T e}{R} \frac{5 + (8\nu - 1)\frac{E_T}{E}}{\left(1 + 3\frac{E_T}{E}\right) \left[5 - (1 - 2\nu)\frac{E_T}{E}\right]} \max_{\tau \ge 0} f(\tau, \rho) \quad (83) $$

where

$$ \max(f) = \sqrt{3} \neq \max(-f) \quad (84) $$

The extremum value for $f$ is reached at $\tau = 0$, i.e. on the whole interval $\rho \in [-1/2, 1/2]$.

In the case of compressed plates, Section 5, one has to find the maximum value of a sinusoidal function $f(\tau)$, thus satisfying $\max(-f) = \max(f)$. One thereby gets the same value for $\lambda_1$ whatever the sign chosen for the bifurcation mode X. For the free edge cylinder in hand, the extremal values of the damped sinusoidal function $f(\tau, \rho)$ are such that $\max(-f) \neq \max(f)$. Depending on the sign retained for the eigenmode X, i.e. depending on whether the free edge moves outwards or inwards, one obtains two different values for the initial slope $\lambda_1$. More precisely, the $\lambda_1$ value corresponding to $W(0)<0$ is greater than that corresponding

<figure><img src="image_5.png" /><figcaption>Fig. 4. Initial slope of bifurcated branches for a cylinder under axial compression.</figcaption></figure>
to $W(0)>0$. As regards the load-displacement curves, this fact gives rise to two bifurcating (half-)branches with different initial slopes, stemming at the same critical load $\lambda_c$.

### 6.2.2. Simply supported edge solution

Again, the case of simply supported edges is simpler than the free edge case and does not require the infinite length assumption. From the fundamental solution (78) and Eq. (29), one obtains

$$ \lambda_1^{\text{supported edges}} = \frac{2E_T e}{R} \frac{5 + (8\nu - 1)\frac{E_T}{E} + 4\sqrt{3\left[5 - 4\nu - (1 - 2\nu)\frac{2E_T}{E}\right]\frac{E_T}{E}}}{\left(1 + 3\frac{E_T}{E}\right)\left[5 - (1 - 2\nu)\frac{E_T}{E}\right]} \quad (85) $$

Contrary to the free edge case, there is only one single value for the initial slope, whatever the sign of the radial displacements. The three values of initial slope (two for the free edges and one for the simply supported edges) are compared in Fig. 4 versus the moduli ratio $E_T/E$. The numerical values are the same as in the plate problem (Section 5), in addition the mean radius of the cylinder is $R = 100$ mm.

Fig. 4 shows that the initial slope when the free edge is bent inwards is always greater than in the case when the free edge is bent outwards. This stiffening effect on the post-critical behavior is mainly due to the compressive hoop stress in the inward free edge case, giving rise to incipient elastic unloading at the edge $X = 0$ where the exponential displacement field is maximal.

One particular cylinder geometry has been chosen to be numerically computed in [32] and the numerical initial slope obtained compares very well with the analytical one given by Eq. (85). Eventually, it should be mentioned that the case of built-in ends, which is not studied here, provides results similar to those with simply supported ends.

## 7. Conclusion

An analytical approach has been proposed to investigate the elastoplastic buckling and initial post-buckling of plates under uni- or bi-axial loadings and cylinders under axial compression. First the analysis has been carried out in the unified framework of the 3D elastoplastic bifurcation theory, and then the resulting equations have been solved under ad hoc assumptions for plate and shell models.

The following geometries and boundary conditions have been considered:

(a) Simply supported plates under biaxial compression.

(b) Simply supported (and possibly built-in) plates under uniaxial compression.

(c) Free edge cylinders under axial compression.

(d) Simply supported edge cylinders under axial compression.

The following results have been obtained:

(1) The critical loads and the bifurcation modes. Some of the results are in good agreement with those existing in the literature. The results obtained for square simply supported plates under biaxial loading and rectangular simply supported (and possibly built-in) plates under uniaxial loading are to the best of the authors' knowledge not known.

(2) The initial slopes of the bifurcated branches. These are the most significant results of the proposed elastoplastic bifurcation analysis.

## References

[1] Shanley FR. Inelastic column theory. International Journal of Aeronautical Sciences 1947;14:261-7.

[2] Hill R. A general theory of uniqueness and stability in elastic-plastic solids. Journal of the Mechanics and Physics of Solids 1958;6:236-49.

[3] Hutchinson JW. Plastic buckling. Advances in Applied Mechanics 1974;14: 67-144.

[4] Nguyen QS, Stolz C. Sur la méthode de développement asymptotique en flambage plastique (On the asymptotic development method in plastic buckling). Comptes Rendus de l'Académie des Sciences 1985;Paris t. 300(7): 235-8.

[5] Akel S. Sur le flambage des structures élastoplastiques (On the buckling of elastoplastic structures). PhD thesis, Ecole Nationale des Ponts et Chaussées, France; 1987.

[6] Cimetière A, El Koulani A, Léger A. Flambage naissant et post-flambage d'un modèle simple en élastoplasticité (Buckling and post-buckling of an elastic-plastic simple model). Comptes Rendus de l'Académie des Sciences 1994;Paris t. 319(11):1263-9.

[7] El Koulani A. Continuation dans les problèmes à frontières libres de type bifurcations plastiques (Continuation in moving boundary problems of the plastic bifurcation type). PhD thesis, University of Metz, France; 1996.

[8] El Koulani A, Léger A. Solutions bifurquées du problème en vitesses initiales pour une poutre élastoplastique (Bifurcated solutions of the initial velocity problem for an elastic-plastic beam). Comptes Rendus de l'Académie des Sciences 1996;Paris t. 322(10):1007-13.

[9] El Koulani A, Léger A. Comportement post-critique des poutres élastoplastiques : existence et régularité des branches bifurquées (Post-critical behavior of elastoplastic beams: existence and regularity of the bifurcated branches). Comptes Rendus de l'Académie des Sciences 1997;Paris t. 324(11):1307-13.

[10] Cimetière A. Sur la modélisation et le flambage des plaques élastoplastiques (On the modelling and buckling of elastoplastic plates). PhD thesis, University of Poitiers, France; 1987.

[11] Durban D, Zuckerman Z. Elastoplastic buckling of rectangular plates in biaxial compression/tension. International Journal of Mechanical Sciences 1999;41: 751-65.

[12] Batterman SC. Plastic buckling of axially compressed cylindrical shells. AIAA Journal 1965;3(2):316-25.

[13] Batterman SC. Free-edge plastic buckling of axially compressed cylindrical shells. Journal of Applied Mechanics 1968;73-9.

[14] Ore E, Durban D. Elastoplastic buckling of axially compressed circular cylindrical shells. International Journal of Mechanical Sciences 1992;34(9): 727-42.

[15] Mao R, Lu G. Plastic buckling of circular cylindrical shells under combined in-plane loads. International Journal of Solids and Structures 2001;38:741-57.

[16] Singace AA. Axial crushing analysis of tubes deforming in the multi-lobe mode. International Journal of Mechanical Sciences 1999;41:865-90.

[17] Gupta NK, Abbas H. Mathematical modeling of axial crushing of cylindrical tubes. Thin-Walled Structures 2000;38:355-75.

[18] Jones N. Structural impact. Cambridge: Cambridge University Press; 1997.

[19] Lee EH, Liu DT. Finite-strain elastic-plastic theory with application to plane-wave analysis. Journal of Applied Physics 1967;38:19-27.

[20] Lee EH. Elastic-plastic deformations at finite strains. Journal of Applied Mechanics 1969;36:1-6.

[21] Green AE, Naghdi PM. A general theory of an elastic-plastic continuum. Archive for Rational Mechanics and Analysis 1965;18:251-81.

[22] Green AE, Naghdi PM. Some remarks on elastic-plastic deformation at finite strain. International Journal of Engineering Science 1971;9:1219-29.

[23] Naghdi PM, Trapp JA. Restrictions on constitutive equations for finitely deformed elastic-plastic materials. Quarterly Journal of Mechanics and Applied Mathematics 1975;28(1):25-46.

[24] Casey J, Naghdi PM. A remark on the use of the decomposition $F = FeFp$ in plasticity. Journal of Applied Mechanics 1980;47:672-5.

[25] Naghdi PM. A critical review of the state of finite plasticity. Journal of Applied Mathematics and Physics 1990;41:315-94.

[26] Halphen B, Nguyen QS. Sur les matériaux standard généralisés (On the generalized standard materials). Journal de Mécanique 1975;14(1):39-63.

[27] Germain P, Nguyen QS, Suquet P. Continuum thermodynamics. Journal of Applied Mechanics 1983;50:1010-20.

[28] Hutchinson JW. On the post-buckling behavior of imperfection-sensitive structures in the plastic range. Journal of Applied Mechanics 1972;155-62.

[29] Nguyen QS. Stability and non-linear solid mechanics. New York: Wiley; 2000.

[30] Donnell LH. Beams, plates and shells. New York: McGraw-Hill; 1976.

[31] Timoshenko SP, Gere JM. Theory of elastic stability. New York: McGraw-Hill; 1961.

[32] Le Grognec P, Le van A. Elastoplastic bifurcation and collapse of axially loaded cylindrical shells. International Journal of Solids and Structures 2008; 45(1):64-86.