# Elastoplastic bifurcation and collapse of axially loaded cylindrical shells

Philippe Le Grognec, Anh Le Van

► To cite this version:

Philippe Le Grognec, Anh Le Van. Elastoplastic bifurcation and collapse of axially loaded cylindrical shells. International Journal of Solids and Structures, 2008, 45 (1), pp.64 - 86. 10.1016/j.ijsolstr.2007.07.017. hal-01006791

HAL Id: hal-01006791

https://hal.science/hal-01006791

Submitted on 14 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
# Elastoplastic bifurcation and collapse of axially loaded cylindrical shells

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

<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

<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

In this paper, a shell finite element is designed within the total Lagrangian formulation framework to deal with the plastic buckling and post-buckling of thin structures, such as cylindrical shells. First, the numerical formulation is validated using available analytical results. Then it is shown to be able to provide the bifurcation modes—possibly the secondary ones—and describe the complex advanced post-critical state of a cylinder under axial compression, where the theory is no longer operative.

*Keywords:* Bifurcation; Buckling; Elastoplasticity; Finite element method; Axially compressed cylinders

## 1. Introduction

Failure of slender or 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.

The problem is particularly difficult with thicker structures where plasticity may occur before buckling. Shanley (1947) 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 (1958) extended these results to a three-dimensional 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.

More precisely, the critical problem of a cylinder under axial compression was the subject of a lot of experimental and numerical investigations (Bushnell, 1982, 1985). A wide variety of experimental results were
obtained in terms of the boundary conditions, the geometries and the material parameters. On the other hand, different numerical models were devised in order to predict the critical loads as well as the advanced post-buckling behavior. Some pioneering experimental, theoretical and numerical results in the literature will be briefly presented in the sequel.

Large discrepancies have always been observed between theoretical and experimental results. The most important feature is that the experimental critical loads sometimes fall down to less than half of the analytical or numerical values. These differences, due to unavoidable imperfections in experiments, were explained by many authors. Among others, in his asymptotic analysis of plastic post-buckling, Hutchinson (1974) examined the influence of various geometric imperfection types on the buckling and post-critical behavior of several structures. It allows one to determine safety factors for the critical load versus the imperfection size and assess the stability of the bifurcated branch. Despite the significant influence of the imperfections on the buckling behavior of a cylinder under axial compression, we shall not account for any type of imperfection in this paper and will only analyze the critical and post-critical behavior of a perfect cylinder.

From a material point of view, different critical values are obtained in the literature, depending on the plasticity theory considered. Indeed, when numerically (or theoretically) solving the problem of plastic buckling, one can choose between the deformation and the incremental theories, which have each their own advantages and drawbacks. The critical stresses obtained by the deformation theory are fairly far from the ones corresponding to the flow theory. As a matter of fact, the solutions derived from the deformation theory turn out to compare best with the experimental results, although such a theory does not take into account the elastic unloading possibility. On the contrary, the flow theory overpredicts the critical loads. However, the discrepancies can be accounted for and quantified through an analysis of imperfection sensitivity.

More than 50 years ago, Batdorf (1949) and Bijlaard (1949) already tried to explain these differences between the two theories. Later, Christoffersen and Hutchinson (1979) proposed a new "corner yield surface" theory which realized a compromise between both theories and was used in number of numerical developments in the literature. Dealing with axially compressed cylinders, Gellin (1979) analyzed the effects of an axisymmetric modal imperfection on the plastic buckling, using the deformation theory. Tvergaard analyzed the effect of geometric and material parameters on the buckling behavior of a cylinder under axial compression (Tvergaard, 1983a,b). He numerically solved the compressed shells using the phenomenological theory of Christoffersen and Hutchinson, and validated his numerical approach by the experimental results obtained by Batterman, Lee, Horton and Johnson in the sixties. He considered an initial sinusoidal axisymmetric defect and looked for a possible bifurcation towards a non-axisymmetric mode. For relatively thin shells, he found a bifurcation point just before the maximum load. For thicker shells, however, the maximal load occurs simultaneously as the axisymmetric deformed shape is localized at both ends in the form of outward axial half-waves. Consequently, the bifurcation is delayed or may merely disappear (very thick shells display an axisymmetric failure). When the bifurcation precedes the localization, a "diamond" mode is observed, displaying both axial and circumferential waves throughout the structure. Otherwise, the localization is followed by the loss of axisymmetry and gives rise to an "elephant foot" mode with regularly distributed "toes" at both ends of the cylinder. According to Tvergaard, the bifurcation is delayed because of the localization which stiffens the structure. Moreover, Tvergaard showed that the circumferential order of the solution decreases as the radius-to-thickness ratio decreases. Also, he noticed that a weak hardening makes it easier for either an axisymmetric failure or an "elephant foot" mode to occur.

Recently, the analysis of this transition between the axisymmetric mode and the "diamond" mode was conducted on a perfect cylinder by Goto and Zhang (1999) using the flow plasticity theory. It was shown therein that the boundary conditions are a deciding factor for obtaining the previous types of buckling modes. In the case of a moderately thick cylinder with built-in ends, Goto and Zhang found two successive bifurcation points: the first one (on the fundamental branch) near the limit load leading to an axisymmetric mode, the second one (on the decreasing primary bifurcated branch) leading to a localized non-axisymmetric secondary mode, named "elephant foot". The thinner the shell, the earlier the secondary mode occurs. Thus, if the shell is thin enough, the second bifurcation point supersedes the first one, prevents from any localization and gives rise to a "diamond" mode. Conversely, if the shell is thick enough, the second bifurcation point disappears—at least in the computed range—and the failure is very likely axisymmetric. The circumferential wave number decreases with the thickness. Less obvious is the observed fact that an increase of the length-to-radius
ratio makes the second bifurcation occur earlier too, without changing the circumferential order of the deformed shape. In the case of simply supported ends, the wave is very near the ends, so that the localization—if it exists—never leads to a second bifurcation. Eventually, the free end condition produces one single non-axisymmetric bifurcation, which is immediately localized. Again, the circumferential order of the solution decreases with the thickness.

Comparisons between the predictions of both flow and deformation theories and experimental results were then performed in Kyriakides et al. (2005), Bardi and Kyriakides (2006) and Bardi et al. (2006), in the case of axially loaded cylindrical shells, including the effect of anisotropy with the use of the Hill's yield criterion. The flow theory significantly overpredicts both the critical stresses and strains, whereas the deformation theory gives rise to predictions much closer to the experimental values. Furthermore, better predictions were reached with the deformation theory, when including imperfections whose amplitude was only a small fraction of the thickness instead of the full thickness used by Gellin and Tvergaard.

Finally, Peek (2000) derived analytical solutions for the immediate post-critical behavior of a cylinder under axial compression in the axisymmetric case only, considering finite strains.

All the above-mentioned facts show that it is necessary to use a numerical tool to analyze the influence of the geometries, the material parameters and the boundary conditions on the complete non-linear buckling and post-buckling response of compressed shells. All plastic bifurcations are strongly influenced by the constitutive model used. Usually, the deformation theory of plasticity or the corner theory are preferred as they predict more accurate critical loads than the flow theory, although some recent works have tried to show that the flow theory is also able to give good results. In this paper, our aim is not to compare the merits of the different theories. Rather, we adopt the $J_2$ flow theory and aim to test the ability of numerical tools developed in predicting plastic bifurcations and post-bifurcation paths. It will be shown that (i) the proposed formulation allows to obtain primary and secondary buckling modes and to exhibit the localization phenomenon, as observed in experiments; (ii) and that the predicted bifurcation stresses may indeed be higher than the experimental ones. Nothing will be done to compare between numerical and experimental results, since this is a tremendous task which is out of the scope of the paper.

## 2. Outline of the paper

This paper is devoted to the numerical computations of the critical loads, the bifurcation modes and the advanced post-critical behaviors of compressed shells in the elastoplastic range.

The principle of virtual work and the total Lagrangian formulation will be used to build a shell finite element encompassing the geometric non-linearities (finite rotations and finite strains) and the material non-linearities (elastoplastic behavior with non-linear isotropic hardening). The shell element is designed to deal with buckling and advanced post-buckling of thin structures, it will be shown to correctly exhibit a wide variety of bifurcation modes and possibly the secondary ones.

First, the numerical formulation will be validated using available analytical results. Then it will be seen that the numerical computations can provide more complete results on the advanced post-critical state of a cylinder under axial compression, where the theory is no longer operative.

## 3. Shell finite element formulation

### 3.1. Shell kinematics with finite rotations

Let us consider a three-dimensional shell-like body with reference volume $\Omega_0$ at the initial time, subjected to finite quasi-static elastoplastic transformations. In this work, the shell-like body is modeled by the Mindlin–Reissner theory, it is represented by a middle surface and an inextensible director field which is not constrained to be normal to the middle surface, except at the reference configuration.

The reference and current positions of a particle on the middle surface are denoted by $\vec{P}_0(\xi, \eta)$ and $\vec{P}(\xi, \eta)$, respectively, where $(\xi, \eta)$ are the surface parameters. The displacement of the middle surface particle is denoted $\vec{U}(\xi, \eta) = \vec{P}(\xi, \eta) - \vec{P}_0(\xi, \eta)$.

The reference and current positions of any particle in the shell are given respectively by:
$$
\begin{align}
\vec{X} &= \vec{P}_0(\xi, \eta) + \zeta \vec{T}(\xi, \eta) \nonumber \\
\vec{x} &= \vec{P}(\xi, \eta) + \zeta \vec{t}(\xi, \eta) \tag{1}
\end{align}
$$

where $\zeta$ is the through-thickness parameter, $\vec{T}(\xi, \eta)$ and $\vec{t}(\xi, \eta)$ are the reference and current director vectors, respectively. The displacement of any particle is:

$$
\vec{D}(\xi, \eta, \zeta) = \vec{U}(\xi, \eta) + \zeta(\vec{t}(\xi, \eta) - \vec{T}(\xi, \eta)). \quad (2)
$$

All the vectors will be resolved in a fixed Cartesian basis:

$$
\vec{X} = X^i \vec{e}_i, \quad \vec{T} = T^i \vec{e}_i, \quad \vec{x} = x^i \vec{e}_i, \quad \vec{t} = t^i \vec{e}_i \tag{3}
$$

where implicit summations are made on repeated indices.
The deformation gradient is:

$$
\mathbf{F} = \vec{g}_i \otimes \vec{G}^i \tag{4}
$$

where {$\vec{G}^i$}$_{i=1,2,3}$ is the contravariant natural basis in the reference configuration, dual to {$\vec{G}_i$}$_{i=1,2,3}$, and
{$\vec{g}_i$}$_{i=1,2,3}$ is the covariant natural basis in the current configuration:

$$
\begin{align}
\vec{G}_1 &= \frac{\partial \vec{X}}{\partial \xi}, & \vec{G}_2 &= \frac{\partial \vec{X}}{\partial \eta}, & \vec{G}_3 &= \frac{\partial \vec{X}}{\partial \zeta} \notag \\
\vec{g}_1 &= \frac{\partial \vec{x}}{\partial \xi}, & \vec{g}_2 &= \frac{\partial \vec{x}}{\partial \eta}, & \vec{g}_3 &= \frac{\partial \vec{x}}{\partial \zeta}. \tag{5}
\end{align}
$$

Here, the variation of the metric through the shell thickness is neglected, so that the bases {$\vec{G}_i$}$_{i=1,2,3}$ and {$\vec{G}^i$}$_{i=1,2,3}$
are equal to those defined on the middle surface and denoted by {$\vec{A}_i$}$_{i=1,2,3}$ and {$\vec{A}^i$}$_{i=1,2,3}$.
Thus, the metric tensor
is defined, respectively, in its covariant and contravariant components as $A_{ij} = \vec{A}_i \cdot \vec{A}_j$ and $A^{ij} = \vec{A}^i \cdot \vec{A}^j$.

It follows from the definition of the deformation gradient:

$$
\mathbf{F} = (\vec{a}_{\alpha} \otimes \vec{A}^{\alpha} + \vec{t} \otimes \vec{T}) + \zeta \vec{t}_{,\alpha} \otimes \vec{A}^{\alpha} \quad (6)
$$

where the Greek indices take the values 1 or 2.

The motion of the inextensible director is described by:

$$
\vec{t} = \mathbf{R} \cdot \vec{T} \tag{7}
$$

where the rotation tensor $\mathbf{R}$ is orthogonal ($\mathbf{R}^{-1} = \mathbf{R}^T$) and can be associated with the skew-symmetric tensor $\boldsymbol{\Theta}$ through the exponential mapping:

$$
\mathbf{R} = \exp \Theta = \sum_{k=0}^{+\infty} \frac{\Theta^k}{k!}. \qquad (8)
$$

Let us define the axial vector $\vec{\theta}$ by $\Theta \cdot \vec{\theta} = \vec{0}$ and $\Theta \cdot \vec{v} = \vec{\theta} \wedge \vec{v}$, $\forall \vec{v} \in \mathbb{R}^3$. From the equality $\mathbf{R} \cdot \vec{\theta} = \vec{\theta}$, it follows that $\mathbf{R}$ is the rotation about $\vec{\theta}$. Then an explicit expression for Eq. (8) is given by the Rodrigues formula:

$$
\mathbf{R} = \cos \theta \mathbf{I} + \frac{\sin \theta}{\theta} \mathbf{\Theta} + \frac{1 - \cos \theta}{\theta^2} \vec{\theta} \otimes \vec{\theta} \tag{9}
$$

where $\theta = \|\vec{\theta}\|$ is the magnitude of the rotation.

Since $\vec{\theta}$ is perpendicular to $\vec{T}$, the director field $\vec{t}$ can be recast in the simplified form:

$$
\vec{t} = \mathbf{R} \cdot \vec{T} = \cos \theta \vec{T} + \frac{\sin \theta}{\theta} \vec{\theta} \wedge \vec{T}. \quad (10)
$$

Denoting $\vec{\theta} = \theta^i \vec{e}_i$, one gets the components of $\vec{t}$ in the fixed Cartesian basis $\{\vec{e}_i\}_{i=1,2,3}$:
$$ \begin{Bmatrix} t^1 \\ t^2 \\ t^3 \end{Bmatrix} = \cos \theta \begin{Bmatrix} T^1 \\ T^2 \\ T^3 \end{Bmatrix} + \frac{\sin \theta}{\theta} \begin{Bmatrix} \theta^2 T^3 - \theta^3 T^2 \\ \theta^3 T^1 - \theta^1 T^3 \\ \theta^1 T^2 - \theta^2 T^1 \end{Bmatrix}. \quad (11) $$

By defining a new local orthonormal basis {$\vec{E}_i$}$_{i=1,2,3}$ at each point of the middle surface as $\vec{E}_3 = \vec{T}$, $\vec{E}_1 = \frac{\vec{A}_1}{\|\vec{A}_1\|}$ and $\vec{E}_2 = \vec{E}_3 \wedge \vec{E}_1$, and by denoting $\vec{\theta} = \hat{\theta}^1 \vec{E}_1 + \hat{\theta}^2 \vec{E}_2$, Eq. (11) becomes:

$$ \begin{Bmatrix} \hat{t}^1 \\ \hat{t}^2 \\ \hat{t}^3 \end{Bmatrix} = \cos \theta \begin{Bmatrix} 0 \\ 0 \\ 1 \end{Bmatrix} + \frac{\sin \theta}{\theta} \begin{Bmatrix} \hat{\theta}^2 \\ -\hat{\theta}^1 \\ 0 \end{Bmatrix}. \quad (12) $$

## 3.2. Variations of kinematic variables

The linearization of the kinematic variables at a given configuration is performed by using the directional derivative with respect to a small parameter $\epsilon$. The following relations hold for small values of $\epsilon$:

$$ \begin{aligned} \vec{P}_{\epsilon} &= \vec{P} + \epsilon \delta \vec{U} \\ \vec{t}_{\epsilon} &= \mathbf{R}_{\epsilon} \cdot \vec{T} = \mathbf{R}(\vec{\theta} + \epsilon \delta \vec{\theta}) \cdot \vec{T} \end{aligned} \quad (13) $$

which yield:

$$ \begin{aligned} \delta \vec{P} &= \frac{d}{d\epsilon} (\vec{P}_{\epsilon})_{|_{\epsilon=0}} = \delta \vec{U} \\ \delta \vec{t} &= \frac{d}{d\epsilon} (\vec{t}_{\epsilon})_{|_{\epsilon=0}} = \frac{\sin \theta}{\theta} \left[ -(\vec{\theta} \cdot \delta \vec{\theta}) \vec{T} + \delta \vec{\theta} \wedge \vec{T} \right] = \frac{\theta \cos \theta - \sin \theta}{\theta^3} (\vec{\theta} \cdot \delta \vec{\theta}) (\vec{\theta} \wedge \vec{T}). \end{aligned} \quad (14) $$

The latter relation can be written in the local basis {$\vec{E}_i$}$_{i=1,2,3}$:

$$ \begin{aligned} \begin{Bmatrix} \delta \hat{t}^1 \\ \delta \hat{t}^2 \\ \delta \hat{t}^3 \end{Bmatrix} &= [\Phi] \begin{Bmatrix} \delta \hat{\theta}^1 \\ \delta \hat{\theta}^2 \end{Bmatrix} \\ [\Phi] &= \begin{bmatrix} -\frac{\sin \theta}{\theta} \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} + \frac{\sin \theta - \theta \cos \theta}{\theta^3} \begin{bmatrix} -\hat{\theta}^2 \hat{\theta}^1 & -\hat{\theta}^2 \hat{\theta}^2 \\ \hat{\theta}^1 \hat{\theta}^1 & \hat{\theta}^1 \hat{\theta}^2 \end{bmatrix} \\ & \qquad + \frac{\sin \theta - \theta \cos \theta}{\theta^3} \begin{bmatrix} 0 & 0 \end{bmatrix} \end{bmatrix} \end{aligned} \quad (15) $$

and then transformed in the Cartesian basis {$\vec{e}_i$}$_{i=1,2,3}$ for computational purposes.

## 3.3. Strains

The Green strain tensor $\mathbf{E}$ is written in the contravariant basis {$\vec{A}^i$}$_{i=1,2,3}$ as:

$$ \mathbf{E} = \frac{1}{2}(\mathbf{F}^T \cdot \mathbf{F} - \mathbf{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 = E_{ij}\vec{A}^i \otimes \vec{A}^j \quad (16) $$

or:

$$ \mathbf{E} = E_{\alpha\beta}\vec{A}^{\alpha} \otimes \vec{A}^{\beta} + E_{\alpha 3}\vec{A}^{\alpha} \otimes \vec{T} + E_{3\alpha}\vec{T} \otimes \vec{A}^{\alpha} + E_{33}\vec{T} \otimes \vec{T}. \quad (17) $$

Since $\vec{a}_\alpha = \vec{A}_\alpha + \vec{U}_{,\alpha}$, the components of $\mathbf{E}$ can be written as functions of the through-thickness parameter $\zeta$:

$$ E_{ij} = E_{ij}^{(0)} + \zeta E_{ij}^{(1)} + \zeta^2 E_{ij}^{(2)} \quad (18) $$

where:
$$
\begin{align}
E_{\alpha\beta}^{(0)} &= \frac{1}{2} (\vec{A}_{\alpha} \cdot \vec{U}_{,\beta} + \vec{A}_{\beta} \cdot \vec{U}_{,\alpha} + \vec{U}_{,\alpha} \cdot \vec{U}_{,\beta}) \nonumber \\
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}) \nonumber \\
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}) \tag{19} \\
E_{\alpha\beta}^{(2)} &= \frac{1}{2} (\vec{t}_{,\alpha} \cdot \vec{t}_{,\beta} - \vec{T}_{,\alpha} \cdot \vec{T}_{,\beta}). \nonumber
\end{align}
$$

All other components, in particular $E_{33}^{(0)}$, $E_{33}^{(1)}$ and $E_{33}^{(2)}$, are zero since the director field is assumed to be inex-
tensible. However, the expression for $E_{33}$ will rather be derived from the plane stress condition, as seen further.
The variation of the Green tensor is computed by:

$$
\delta \mathbf{E} = \frac{d\mathbf{E}(\vec{P}_v, \vec{t}_v)}{d\epsilon} \Big|_{\epsilon=0} \qquad (20)
$$

which yields:

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

3.4. Stresses

The second Piola–Kirchhoff stress tensor $\Sigma$ is resolved in basis $\{\vec{A}_i\}_{i=1,2,3}$:

$$
\Sigma = \Sigma^{ij} \vec{A}_i \otimes \vec{A}_j = \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 (22)
$$

In the sequel, use will be made of the stress resultant tensors defined as:

$$
\begin{align}
\mathbf{n} &= n^{ij} \vec{A}_i \otimes \vec{A}_j \qquad n^{ij} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \Sigma^{ij} d\zeta \notag \\
\mathbf{m}^{(1)} &= m^{ij(1)} \vec{A}_i \otimes \vec{A}_j \qquad m^{ij(1)} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \Sigma^{ij} \zeta d\zeta \tag{23} \\
\mathbf{m}^{(2)} &= m^{ij(2)} \vec{A}_i \otimes \vec{A}_j \qquad m^{ij(2)} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \Sigma^{ij} \zeta^2 d\zeta. \notag
\end{align}
$$

3.5. Elastoplasticity formulation

The equation set for finite elastoplasticity is established on the basis of the total Lagrangian formulation
and the generalized standard materials theory introduced by Halphen and Nguyen (1975). Use is made of
the von Mises yield criterion with a non-linear isotropic hardening:

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

where $\Sigma^d$ denotes the deviatoric part of stress tensor $\Sigma$, $\sigma_0$ the initial yield stress and $a$ the thermodynamic force characterizing the isotropic hardening.
Following Green and Naghdi (Green and Naghdi, 1965, 1971), the Green strain tensor E is split additively into its elastic and plastic parts:

$$ \mathbf{E} = \mathbf{E}^e + \mathbf{E}^p. \qquad (25) $$

The elastic strain is related to the second Piola–Kirchhoff stress by the three-dimensional Saint–Venant–Kirchhoff law:

$$ \Sigma = \mathbf{D} : \mathbf{E}^e = \lambda \operatorname{tr}(\mathbf{E}^e) \mathbf{I} + 2\mu \mathbf{E}^e \qquad (26) $$

where $\lambda$ and $\mu$ are the Lamé parameters, related to Young's modulus E and Poisson's ratio $\nu$ by:

$$ \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)} \qquad \mu = \frac{E}{2(1 + \nu)}. \qquad (27) $$

For computational purposes, all the equations are written in matrix representations using the covariant or contravariant bases. For instance, the relationship (26) reads in the covariant and contravariant bases:

$$ \Sigma^{ij} = D^{ijkl} E_{lk}^{e} = (\lambda A^{ij} A^{kl} + \mu (A^{ik} A^{jl} + A^{il} A^{kj})) E_{lk}^{e}. \qquad (28) $$

The thermodynamic force $a$ in (24) is derived from the so-called hardening energy $w^\alpha(\alpha)$, where $\alpha$ is the hardening variable:

$$ a = \frac{\partial w^\alpha(\alpha)}{\partial \alpha}. \qquad (29) $$

The flow rule and the hardening law are derived from the yield function:

$$ \dot{\mathbf{E}}^p = \dot{\lambda} \frac{\partial f}{\partial \Sigma} = \dot{\lambda} \frac{\Sigma^d}{\sqrt{\frac{2}{3} \Sigma^d : \Sigma^d}} \qquad \dot{\alpha} = -\dot{\lambda} \frac{\partial f}{\partial a} = \dot{\lambda} \qquad (30) $$

where $\dot{\lambda}$ is the plastic multiplier. Use will also be made of the equivalent plastic strain $P$ defined by its rate:

$$ \dot{P} = \sqrt{\frac{2}{3} \dot{\mathbf{E}}^p : \dot{\mathbf{E}}^p} = \dot{\lambda} = \dot{\alpha}. \qquad (31) $$

Since the normal material stress is assumed to be zero, $\Sigma^{33} = 0$, the elastic fourth-order tensor $\mathbf{D}$ must be replaced with the reduced one $\mathbf{C}$ for the shell model:

$$ C^{ijkl} = D^{ijkl} - \frac{D^{ij33} D^{33kl}}{D^{3333}} \quad (i, j) \neq (3, 3), (k, l) \neq (3, 3). \qquad (32) $$

Substituting Eq. (18) into (23) gives the expressions for the stress resultant tensors (23):

$$ \begin{align} n^{ij} &= \int_{-\frac{h}{2}}^{\frac{h}{2}} C^{ijkl} (E_{lk}^{e(0)} + \zeta E_{lk}^{e(1)} + \zeta^2 E_{lk}^{e(2)}) d\zeta \notag \\ m^{ij(1)} &= \int_{-\frac{h}{2}}^{\frac{h}{2}} C^{ijkl} (E_{lk}^{e(0)} + \zeta E_{lk}^{e(1)} + \zeta^2 E_{lk}^{e(2)}) \zeta d\zeta \tag{33} \\ m^{ij(2)} &= \int_{-\frac{h}{2}}^{\frac{h}{2}} C^{ijkl} (E_{lk}^{e(0)} + \zeta E_{lk}^{e(1)} + \zeta^2 E_{lk}^{e(2)}) \zeta^2 d\zeta \notag \\ (i, j) \neq (3, 3), (k, l) \neq (3, 3). \notag \end{align} $$

According to the thin shell assumptions, the strain components $E_{ij}^{e(2)}$ in Eq. (18) will be neglected in the sequel.

### 3.6. The discretized non-linear equations

Use is made of the principle of virtual work in order to obtain the discretized non-linear equations:

$$ \delta W_{\text{int}} + \delta W_{\text{ext}} = 0 \qquad (34) $$
where:

$$ \delta W_{\text{int}} = - \int_{\Omega_0} \Sigma : \delta \mathbf{E} \, dV = - \int_{\mathcal{M}} \left( \int_{-\frac{h}{2}}^{\frac{h}{2}} \Sigma^{ij} \left( \delta E_{ji}^{(0)} + \zeta \delta E_{ji}^{(1)} \right) d\zeta \right) dA \quad (35) $$

and:

$$ \delta W_{\text{ext}} = \int_{\Omega_0} \vec{f}^P \cdot \delta \vec{P} \, dV + \int_{\Omega_0} \vec{f}^t \cdot \delta \vec{t} \, dV + \int_{\partial \Omega_0} \vec{\tau}^P \cdot \delta \vec{P} \, dS + \int_{\partial \Omega_0} \vec{\tau}^t \cdot \delta \vec{t} \, dS \quad (36) $$

where $\mathcal{M}$ represents the shell middle surface, and $\vec{f}$ and $\vec{\tau}$ are the volume and surface distributed loads, respectively.

Substituting Eqs. (35) and (36) into (34) gives:

$$ - \int_{\mathcal{M}} \left( n^{ij} \delta E_{ji}^{(0)} + m^{ij(1)} \delta E_{ji}^{(1)} \right) dA + \delta W_{\text{ext}} = 0 \quad (37) $$

where the components $n^{ij}$ and $m^{ij(1)}$ are defined in (23).

The middle surface is discretized using isoparametric shell finite elements. There are five degrees of freedom at each node: three translations and two rotations in the tangent space. The reference position $\vec{P}_0$, the displacement $\vec{U}$ on the middle surface, and the director vectors $\vec{T}$ and $\vec{t}$ are interpolated in the same way:

$$ \begin{align} \vec{P}_0(\xi, \eta) &= \sum_{i=1}^{N} N^i(\xi, \eta) \vec{P}_{0i} \notag \\ \vec{T}(\xi, \eta) &= \sum_{i=1}^{N} N^i(\xi, \eta) \vec{T}_i \notag \\ \vec{U} &= \sum_{i=1}^{N} N^i(\xi, \eta) \vec{U}_i \notag \\ \vec{t} &= \sum_{i=1}^{N} N^i(\xi, \eta) \vec{t}_i \tag{38} \end{align} $$

where $N$ is the element node number and $(N^i)_{i \in [1, N]}$ are the shape functions.

Attention should be drawn on two types of errors due to the geometric discretization:

- The interpolation of the reference position $\vec{P}_0$ of a particle on the middle surface leads to several normals at each node, one per contiguous element.

- The interpolated reference director vector $\vec{T}$ does not coincide with the normal at the surface defined by the interpolation of $\vec{P}_0$.

Let us mention that it is necessary to define a local orthonormal basis at each node in correspondence with the Cartesian fixed basis in order to describe $\vec{\theta}$ in every adjoined element in a coherent way.

The tangent vectors on the middle surface are computed from the interpolation of $\vec{P}_0$:

$$ \begin{align} \vec{A}_1 &= \vec{P}_{0,1} = \sum_{i=1}^{N} N_{,1}^i \vec{P}_{0i} \notag \\ \vec{A}_2 &= \vec{P}_{0,2} = \sum_{i=1}^{N} N_{,2}^i \vec{P}_{0i}. \tag{39} \end{align} $$

The following variations and derivatives of these fields will be required:
$$
\begin{align}
\vec{U}_{,\alpha} &= \sum_{i=1}^{N} N_{,\alpha}^{i} \vec{U}_{i} \nonumber \\
\delta \vec{U}_{,\alpha} &= \sum_{i=1}^{N} N_{,\alpha}^{i} \delta \vec{U}_{i} \nonumber \\
\vec{t}_{,\alpha} &= \sum_{i=1}^{N} N_{,\alpha}^{i} \vec{t}_{i} \tag{40} \\
\delta \vec{t} &= \sum_{i=1}^{N} N^{i} \delta \vec{t}_{i} \nonumber \\
\delta \vec{t}_{,\alpha} &= \sum_{i=1}^{N} N_{,\alpha}^{i} \delta \vec{t}_{i}. \nonumber
\end{align}
$$

The second derivatives will also be required:

$$
\begin{align}
\Delta(\delta\vec{t}) &= \sum_{i=1}^{N} N^i \Delta(\delta\vec{t}_i) \tag{41} \\
\Delta(\delta\vec{t}_{,\alpha}) &= \sum_{i=1}^{N} N_{,\alpha}^i \Delta(\delta\vec{t}_i). \nonumber
\end{align}
$$

It should be noted that $\Delta(\delta\vec{U})$ is always zero. The integrations over the middle surface and through the thickness are performed using the Gaussian quadrature. Eventually, the virtual work principle (37) yields the non-linear discretized equation system:

$$
\{\mathbf{R}(\mathbf{U})\} = \{\Psi(\mathbf{U})\} - \{\Phi\} = \{0\} \quad (42)
$$

where {$\mathbf{U}$} denotes the nodal displacement vector of the whole structure, {$\mathbf{R}$} the residual vector, {$\Psi$} and {$\Phi$}
are the internal and external force vectors, respectively.

**4. Solution procedure**

4.1. *The Newton–Raphson procedure*

The non-linear equation system (42) is solved by the iterative Newton–Raphson procedure, which requires
the computation of the structural tangent stiffness matrix:

$$
[\mathbf{K}_T] = \frac{\partial\{\mathbf{R}(\mathbf{U})\}}{\partial\{\mathbf{U}\}} = \frac{\partial\{\Psi(\mathbf{U})\}}{\partial\{\mathbf{U}\}}. \quad (43)
$$

The geometric part of [K<sub>T</sub>] requires the computation of Δ(δ⃗t<sub>i</sub>) at the middle surface nodes and entails a com-
plicated development. The interested reader is referred to Brank et al. (1997) for more details.

At each iteration, the local integration consists in solving the following problem: assuming the strain **E** is
known, compute the plastic strain **E**<sup>P</sup>, the stress **Σ** and the reduced consistent elastoplastic tangent modulus **C**<sup>P</sup>
as in (32). The implicit Euler scheme and an algorithm of radial return type are employed to integrate the state
and evolution laws as described in Simo and Hughes (1998). In particular, the von Mises criterion (24) is
rewritten as:

$$
f = \frac{1}{2} \Sigma^{\mathrm{T}} \mathbb{P} \Sigma - \frac{1}{3} (\sigma_{0} + a)^{2} \quad (44)
$$

where $\mathbb{P}$ is the deviatoric operator ($\Sigma^d = \mathbb{P}\Sigma$). At each step $n$, after computing the prediction $\Sigma^E$ and the consistency parameter $\lambda$, the consistent elastoplastic tangent operator for computing the structural tangent stiffness matrix is obtained by:
$$C_n^p = \frac{\partial \Sigma_n}{\partial E_n} = \Gamma(\lambda) - \frac{(\Gamma(\lambda)\mathbb{P}\Sigma_n)(\Gamma(\lambda)\mathbb{P}\Sigma_n)^T}{\Sigma_n^T \mathbb{P}^T \Gamma(\lambda)\mathbb{P}\Sigma_n + \mu} \quad (45)$$

where:

$$\begin{align*}
\Gamma(\lambda) &= (C^{-1} + \lambda \mathbb{P})^{-1} \\
\Sigma_n &= \Gamma(\lambda) C^{-1} \Sigma^E \\
P_n &= P_{n-1} + \lambda \sqrt{\frac{2}{3} \Sigma_n^T \mathbb{P} \Sigma_n} \\
\mu &= \frac{2a'(P_n)}{3\gamma} \Sigma_n^T \mathbb{P} \Sigma_n \quad (\gamma = 1 - \frac{2}{3}\lambda a'(P_n)).
\end{align*} \quad (46)$$

In the above, C denotes the matrix of tensor C in Eq. (32).

It should be noticed that the plane stress condition $\Sigma^{33} = 0$ is a non-trivial constraint for the non-linear equation set which is strain driven.

## 4.2. Arc-length method

The nodal displacement vector {$\tilde{\mathbf{U}}$} is split into two parts as shown in Eq. (47) below: one denoted by {$\tilde{\mathbf{U}}$} contains the unknown degrees of freedom, the other denoted by {$\tilde{\mathbf{U}}$} contains the prescribed degrees of freedom. The external force and the residual vector in Eq. (42) are split in a similar way. For instance, the external force {$\tilde{\mathbf{\Phi}}$} is seen as the union of two parts: {$\tilde{\mathbf{\Phi}}$} corresponding to {$\tilde{\mathbf{U}}$} contains the prescribed force components, and {$\tilde{\mathbf{\Phi}}$} corresponding to {$\tilde{\mathbf{U}}$} contains the unknown reaction force components. Likewise, the tangent stiffness matrix {$\tilde{\mathbf{K}}_T$} is partitioned into four sub-matrices:

$$\{\tilde{\mathbf{U}}\} = \begin{Bmatrix} \tilde{\mathbf{U}} \\ \tilde{\mathbf{U}} \end{Bmatrix}, \quad \{\tilde{\mathbf{R}}\} = \begin{Bmatrix} \tilde{\mathbf{R}} \\ \tilde{\mathbf{R}} \end{Bmatrix}, \quad \{\tilde{\mathbf{\Phi}}\} = \begin{Bmatrix} \tilde{\mathbf{\Phi}} \\ \tilde{\mathbf{\Phi}} \end{Bmatrix}, \quad [\tilde{\mathbf{K}}_T] = \begin{bmatrix} \tilde{\tilde{\mathbf{K}}}_T & \tilde{\tilde{\mathbf{K}}}_T \\ \tilde{\tilde{\mathbf{K}}}_T & \tilde{\tilde{\mathbf{K}}}_T \end{bmatrix}. \quad (47)$$

Either the prescribed displacement or the external loading is assumed to be proportional:

$$\{\tilde{\mathbf{U}}\} = \{\tilde{\mathbf{U}}\}^0 + \lambda \{\tilde{\mathbf{U}}\}^{\text{ref}} \quad (48a)$$

$$\{\tilde{\mathbf{\Phi}}\} = \lambda \{\tilde{\mathbf{\Phi}}\}^{\text{ref}} \quad (48b)$$

where $\lambda$ is the control parameter, {$\tilde{\mathbf{U}}$}^{ref} and {$\tilde{\mathbf{\Phi}}$}^{ref} denote reference prescribed quantities. Vector {$\tilde{\mathbf{U}}$}^0 related to zero prescribed displacements does not change the value of {$\tilde{\mathbf{U}}$}. At each increment, the quadratic arc-length method is used in order to proceed on the solution branches given a specified arc-length $\Delta l$, see Wempner (1971), Riks (1979) and Crisfield (1991). The constraint equation is either of the following relations, depending on whether one has a proportional prescribed displacement or loading:

$$\|\Delta\tilde{\mathbf{U}}\|^2 + \Delta\lambda^2 \|\tilde{\mathbf{U}}^{\text{ref}}\|^2 = \Delta l^2 \quad (49a)$$

$$\|\Delta\tilde{\mathbf{U}}\|^2 + \Delta\lambda^2 C_{\text{ref}}^2 = \Delta l^2 \quad (49b)$$

where $\Delta\tilde{\mathbf{U}} = \tilde{\mathbf{U}} - \tilde{\mathbf{U}}_{n-1}$ is the difference between the displacement $\tilde{\mathbf{U}}$ at the current increment (n) and that at the previous increment (n-1). Similarly, $\Delta\lambda = \lambda - \lambda_{n-1}$. In Eq. (49b), the scalar $C_{\text{ref}}$ is a scale factor which makes the relation consistent dimensionally.

Combining Eq. (49a,b) with the equilibrium equation leads to a quadratic equation. When solving this quadratic equation, one may encounter severe computational difficulties due to complex roots which occur repeatedly. It is found that an efficient way to cope with these complex roots is to modify the standard solution scheme according to Lam and Morley (1992). Let us explain the procedure in the case of prescribed loadings (48b). The main idea is to project the residual force onto the external load vector. At a current iteration where complex roots occur, the residual force is split into one component in the load direction and another
component orthogonal to this load. The last component is mainly responsible for the complex roots and should be eliminated. The standard iteration loop is modified as follows.

Iteration loop

[... ]

In the standard arc-length method, solve the quadratic equation.

If the roots are complex, then:

(i) Compute $g = \tilde{\mathbf{R}}^T \tilde{\Phi}^{\text{ref}} / \tilde{\Phi}^{\text{ref}T} \tilde{\Phi}^{\text{ref}}$ and $\mathbf{H} = \tilde{\mathbf{R}} - g \tilde{\Phi}^{\text{ref}}$.

(ii) Compute:

$$ \delta \tilde{\mathbf{U}}_H = - \tilde{\mathbf{K}}_T^{-1} \mathbf{H} \quad \text{and} \quad \tilde{\mathbf{U}}_{\phi \text{ref}} = \tilde{\mathbf{K}}_T^{-1} \tilde{\Phi}_{\text{ref}}. \qquad (50) $$

(iii) Solve a quadratic equation in $\eta$: $\alpha\eta^2 + 2\beta\eta + \gamma = 0$, where the coefficients $\alpha, \beta$ and $\gamma$ depend on $g, C_{\text{ref}}^2$, $\tilde{\mathbf{U}}_{\phi \text{ref}}$ and $\delta \tilde{\mathbf{U}}_H$. Hopefully, this equation in $\eta$ gives real roots $\eta_1$ and $\eta_2$. Choose a value for $\eta$ at 5% of $|\eta_2 - \eta_1|$ from $\eta_1$ or $\eta_2$.

(iv) Solve a quadratic equation in $\xi$: $A\xi^2 + 2B(\eta)\xi + C(\eta) = 0$, where the coefficient $A$ depends on $C_{\text{ref}}^2$ and $\tilde{\mathbf{U}}_{\phi \text{ref}}, B(\eta)$ and $C(\eta)$ are functions of $\eta, g, C_{\text{ref}}^2, \tilde{\mathbf{U}}_{\phi \text{ref}}$ and $\delta \tilde{\mathbf{U}}_H$. This equation in $\xi$ gives real roots. Choose that root $\xi$ which enables the solution point to advance in the desired direction.

(v) Make the corrections for the ith-iteration: $\delta\lambda = g + \xi$ and $\delta\tilde{\mathbf{U}} = \eta\delta\tilde{\mathbf{U}}_H + \xi\tilde{\mathbf{U}}_{\phi \text{ref}}$.

End if

[... ]

End of iteration loop.

The modified arc-length method described above is also applicable to the case of prescribed displacements (48a), providing that $\tilde{\Phi}_{\text{ref}}$ is replaced with $\tilde{\mathbf{K}}_T \tilde{\mathbf{U}}_{\text{ref}}, C_{\text{ref}}^2$ with $\|\tilde{\mathbf{U}}_{\text{ref}}\|^2$ and $\tilde{\mathbf{U}}_{\phi \text{ref}}$ with:

$$ \tilde{\mathbf{U}}_{\text{ref}} = - \tilde{\mathbf{K}}_T^{-1} \tilde{\mathbf{K}}_T \tilde{\mathbf{U}}_{\text{ref}}. \qquad (51) $$

For numerical purposes, the computer program is developed so as to handle both prescribed loads and prescribed displacements in a unified way.

### 4.3. Branching method

The branch switching techniques are included in the numerical procedure in order to detect the bifurcation points and bifurcate onto a given branch. Specific methods are implemented following Riks (Riks, 1979, 1991) and Seydel (1994). The four fundamental steps are given below.

1. At the end of each increment, it must be checked whether one has gone across one or several critical points. The detection of critical points is based on the singularity of the tangent stiffness matrix, which is factorized following the Crout formula $\mathbf{K}_T = \mathbf{L}\mathbf{D}\mathbf{L}^T$, where $\mathbf{L}$ is a lower triangular matrix with unit diagonal elements and $\mathbf{D}$ is a diagonal matrix. Since the number of negative eigenvalues of $\mathbf{K}_T$ is equal to the number of negative diagonal elements (pivots) of $\mathbf{D}$, the critical points are determined by counting the negative pivot number.

2. Each critical point detected has to be isolated in order to determine its nature: limit point or bifurcation point. To do this, the current arc-length $\Delta l$ is re-estimated several times using a dichotomy-like method. In the case of prescribed loadings, a simple way to distinguish a limit point from a bifurcation point is to determine the sign of the current stiffness parameter introduced by Bergan et al. (1978):

$$ k = \frac{\{\tilde{\Phi}^{\text{ref}}\}^T \{\tilde{\mathbf{U}}_{\phi \text{ref}}\}}{\{\tilde{\mathbf{U}}_{\phi \text{ref}}\}^T \{\tilde{\mathbf{U}}_{\phi \text{ref}}\}} \qquad (52) $$
where {$\tilde{\Phi}^{\text{ref}}$} is defined in (48b) and {$\tilde{\mathbf{U}}_{\text{ref}}$} in (50). The sign of parameter $k$ changes when passing a limit point, whereas it remains unchanged when passing a bifurcation point. In the case of prescribed displacements, parameter $k$ is computed by:

$$k = - \frac{\{\tilde{\mathbf{K}}_{\text{T}} \tilde{\mathbf{U}}^{\text{ref}}\}^{\text{T}} \{\tilde{\mathbf{U}}_{\text{ref}}\}}{\{\tilde{\mathbf{U}}_{\text{ref}}\}^{\text{T}} \{\tilde{\mathbf{U}}_{\text{ref}}\}} \qquad (53)$$

where {$\tilde{\mathbf{U}}^{\text{ref}}$} is defined in (48a) and {$\tilde{\mathbf{U}}_{\text{ref}}$} in (51).

3. If it is a bifurcation point, the step increment is renewed so as to reach a point just behind the bifurcation point.

4. Finally, the switching on a bifurcated branch is performed using the mode injection method (Wagner and Wriggers, 1988; Seydel, 1994): at the first step of a bifurcating branch, the eigenvector $\tilde{\mathbf{Z}}$, solution of $\tilde{\mathbf{K}}_{\text{T}}\tilde{\mathbf{Z}} = \mathbf{0}$, is computed and the following predictions are used:

$$\delta\lambda = 0 \qquad \delta\tilde{\mathbf{U}} = \pm \frac{\tilde{\mathbf{Z}}}{\|\tilde{\mathbf{Z}}\|} \qquad (54)$$

## 5. Numerical results

Let us consider a cylindrical tube defined in the reference configuration by the length $l$ along the x-axis, the average radius $R$ and the thickness $t \ll R$. One end of the cylinder lies in the plane $X=0$ (more details will be given according to the boundary condition type) whereas the other end at $X=l$ is loaded by a uniformly distributed nominal axial compressive stress $\Pi_{XX} = -\lambda < 0$.

First, comparisons will be made between the numerical results and available analytical solutions, in the case of axisymmetric modes with a linear isotropic hardening. Next, still considering axisymmetric modes, we will compare our numerical results with those obtained by Brank et al. (1997) using the flow plasticity theory with a non-linear isotropic hardening. Eventually, we will go further in the description of the post-buckling behavior and obtain both axisymmetric and non-axisymmetric deformed shapes, especially the so-called “diamond” modes or the “elephant foot” secondary modes. It will be shown how the geometric and material parameters, together with the boundary conditions, affect the buckling mode type. For simplicity sake, the numerical results will be presented assuming a linear isotropic hardening. The obtained results show that the proposed formulation is able to provide complex experimental modes which are hard to predict in theory.

In all examples, the Newton-Raphson method is used for solving the matrix non-linear equations of the problem. The path-following is carried out either by displacement (49a) or by force control (49b), use is also made of the extended version of the arc-length method in Lam and Morley (1992) as explained in Section 4 to deal with the complex roots in the solution scheme.

A few preliminary mesh convergence tests lead us to use about 200 quadratic eight-node elements for the considered portion of the cylinder and four Gaussian points for the numerical quadrature through the shell thickness.

### 5.1. Comparison between analytical and numerical results

The comparison between analytical and numerical results is possible within the context of a linear isotropic hardening and axisymmetric bifurcation modes. In order to simplify the discussion, let us take the yield stress $\sigma_0$ small enough so that the plasticization occurs before the tangent modulus critical load is reached. Furthermore, it is assumed that the hardening energy $w^\alpha(\alpha)$ in Eq. (29) is a quadratic function, i.e. $w^\alpha = \frac{1}{2}h\alpha^2$, where constant $h$ is the hardening modulus. The numerical values are shown in Tables 1 and 2. The hardening modulus $h$ is given a very high value just to make sure of obtaining an axisymmetric mode, so that the numerical results can be compared with the analytical solution.
<table><caption>Table 1
Geometry used in Section 5.1</caption><thead><tr><td>Length l</td><td>Radius R</td><td>Thickness t</td></tr></thead><tbody><tr><td>2 m</td><td>1 m</td><td>1 cm</td></tr></tbody></table>

<table><caption>Table 2
Material properties used in Section 5.1</caption><thead><tr><td>Young's modulus E</td><td>Poisson's ratio v</td><td>Yield stress σ₀</td><td>Isotropic hardening modulus h</td></tr></thead><tbody><tr><td>2.1 × 10¹¹ Pa</td><td>0.3</td><td>Free edges: 1 × 10⁸ Pa<br>Simply supported edges: 3.5 × 10⁸ Pa</td><td>2 × 10¹⁰ Pa</td></tr></tbody></table>

### 5.1.1. *Free edges*

In the case of free edge cylinders, the yield stress $σ₀$ is supposed to be small enough for the buckling to occur in the plastic range.

The analytical results for a cylinder of infinite length can be obtained by describing the shell kinematics with the Koiter model, which is equivalent to the Love–Kirchhoff plate model. The critical stress was given by Batterman (1968):

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

where $E_T$ is the tangent modulus, related to the Young's modulus $E$ and the isotropic hardening modulus $h$ by $\frac{1}{E_T} = \frac{1}{E} + \frac{1}{h}$. Batterman showed that the critical stress is a little bit smaller for a finite shell and it hardly depends on the length. The eigenmode is:

$$ \left\{ \begin{array}{l} U = -2t \sqrt{\frac{3t}{R}} \frac{1}{p} \frac{1-(1-2v)^{\frac{E_T}{E}}}{1+3^{\frac{E_T}{E}}} \sin\left(\frac{pX}{\sqrt{Rt}}\right) \exp\left(-\frac{pX}{\sqrt{3Rt}}\right) \\ V = 0 \\ W = t\left(\sqrt{3}\cos\left(\frac{pX}{\sqrt{Rt}}\right) - \sin\left(\frac{pX}{\sqrt{Rt}}\right)\right) \exp\left(-\frac{pX}{\sqrt{3Rt}}\right) \end{array} \right. \quad (56) $$

where $U$, $V$ and $W$ are axial, circumferential and radial displacements, respectively, and:

$$ p = 4 \sqrt{\frac{27 \left[ 5 - 4v - (1 - 2v)^2 \frac{E_T}{E} \right] \frac{E_T}{E}}{\left( 1 + 3 \frac{E_T}{E} \right)^2}} \quad (57) $$

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

The numerical computations show that an axisymmetric bifurcation under increasing load is obtained at the critical value of 184 MPa, which is identical with the analytical critical stress $\lambda_c^{\text{free edges}}$ derived from Eq. (55), within 1%. The bifurcation shape is a sinusoid weighted by a decreasing exponential function, so that its magnitude is minimal at the middle of the cylinder and maximal at the ends. This result is in full accordance with Batterman's results (Batterman, 1968). Fig. 1 shows that the numerical modal shape is very close to the Batterman's analytical one corresponding to the semi-infinite cylinder.

### 5.1.2. *Simply supported edges*

The case of simply supported edges is simpler and does not require the infinite length assumption. Let the radial displacements be fixed and the rotations be free on both edges.

Batterman (1965) showed that, if the cylinder is long enough, a good approximate expression for the critical value is:
<figure><img src="image_2.png" /><figcaption>Fig. 1. Analytical and numerical axisymmetric bifurcation modes of an elastoplastic cylinder under axial compression with free edges.</figcaption></figure>

$$ \lambda_c^{\text{supported edges}} = \frac{2Et}{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 (58) $$

In the case of simply supported edges, it is also possible to obtain the analytical slope to the post-critical curve at the bifurcation point. Since the calculation is somewhat intricate, we give the final result without proof:

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

In the case of simply supported edges, the numerical critical value is found to be twice as large as in the free edge case, which is in agreement with the analytical solution $\lambda_c^{\text{supported edges}}$ obtained from Eq. (58). From a qualitative point of view, the bifurcation comes just before the maximum load value, as commonly observed in other numerical applications in the literature.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Elastoplastic cylinder under axial compression with simply supported edges: compressive stress versus axial shortening.</figcaption></figure>
Fig. 2 shows the load versus the axial shortening in the simply supported edge case. The analytical tangent to the post-critical curve at the bifurcation point is also plotted using Eq. (59) and fits well with the numerical curve.

## 5.2. Comparison with a numerical result of the literature

We still consider axisymmetric modes and now compare the numerical results with those obtained by Brank et al. (1997) using the flow plasticity theory with a non-linear isotropic hardening. The computations are carried out with the same geometry and material as in Brank et al. (1997), see Tables 3 and 4. The coefficients in the non-linear isotropic hardening law are determined by fitting the data given in Brank et al. (1997), see Table 5.

The computations are performed for two boundary conditions: simply supported edges and built-in edges. The bifurcation loads are 280 MPa in the case of simply supported edges and 300 MPa in the case of built-in edges, which agree very well with the solutions of Brank et al. (1997). As shown in Fig. 3, the bifurcation loads are close to the maximum loads. Both boundary conditions lead to similar bifurcation modes, see Fig. 4. In addition, Fig. 3 shows that although the bifurcation loads are similar, the corresponding axial shortenings are significantly distinct.

Fig. 5 shows the maximum equivalent plastic strain *P* in the shell thickness. Fig. 5a is related to the case of simply supported edges when the axial shortening equals 1.6%; Fig. 5b to the case of built-in edges when the axial shortening equals 2.3%. As regards the maximum plastic strain in the whole shell, it is 5.6% in the case of simply supported edges and 6.7% in the case of built-in edges. The points where the plastic strain *P* reaches the maximum are also indicated in Fig. 5. For both boundary condition types, it is found that *P* reaches its maximum value near the ends of the cylinder where the axisymmetric localization takes place.

## 5.3. Influence of thickness, hardening modulus and boundary conditions

In the previous Sections 5.1 and 5.2, the numerical formulation has been validated by comparing analytical versus numerical results and numerical versus numerical results. Attention is next focused on more complex modes. Several thicknesses, isotropic hardening moduli and boundary conditions (free, simply supported or built-in edges, prescribed displacements or stresses) are considered in order to exhibit the variety of post-critical responses of the cylinder under axial compression. It is expected that the boundary conditions have a strong influence on the deformed shapes, especially for non-axisymmetric modes. Indeed, in the particular case of free edges, prescribing the axial displacements at the ends of the cylinder leads more naturally to an axisym-

<table><caption>Table 3
Geometry used in Section 5.2</caption><thead><tr><td>Length l</td><td>Radius R</td><td>Thickness t</td></tr></thead><tbody><tr><td>1.271 × 10<sup>-1</sup> m</td><td>4.333 × 10<sup>-2</sup> m</td><td>2.36 × 10<sup>-3</sup> m</td></tr></tbody></table>

<table><caption>Table 4
Material properties used in Section 5.2</caption><thead><tr><td>Young's modulus E</td><td>Poisson's ratio v</td><td>Yield stress σ<sub>0</sub></td><td>Isotropic hardening coefficients</td></tr></thead><tbody><tr><td>2.16962 × 10<sup>11</sup> Pa</td><td>0.274</td><td>1.62722 × 10<sup>8</sup> Pa</td><td>Obtained by fitting values in Brank et al. (1997)</td></tr></tbody></table>

<table><caption>Table 5
Non-linear isotropic hardening curve</caption><tbody><tr><td>P (10<sup>-3</sup>)</td><td>0</td><td>1</td><td>2.55</td><td>4.25</td><td>6.75</td><td>9.25</td><td>14.25</td><td>24.25</td><td>54.25</td></tr><tr><td>K(P) (MPa)</td><td>0</td><td>68.95</td><td>86.19</td><td>106.2</td><td>123.4</td><td>133.1</td><td>151</td><td>166.9</td><td>211</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Elastoplastic cylinder under axial compression: compressive stress versus axial shortening. Case of non-linear isotropic hardening.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Axisymmetric bifurcation modes of an elastoplastic cylinder under axial compression. Case of non-linear isotropic hardening.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Maximum equivalent plastic strain in the thickness of the cylinder under axial compression. Case of non-linear isotropic hardening: (a) simply supported edges, axial shortening = 1.6%; (b) built-in edges, axial shortening = 2.3%.</figcaption></figure>
metric mode, whereas prescribing the stresses leads to non-axisymmetric modes as the ends are free to move now.

The numerical results in this section are obtained using the geometry and the material properties shown in Tables 6 and 7. The length-to-radius ratio is constant and equal to $l/R = 3.33$, while several radius-to-thickness ratios will be used.

**Remark.** Why the unrealistically high values of the hardening modulus $h$ in Table 7? It is found that with usual $h$ values, the numerical computations in plasticity give axisymmetric modes. Only very high moduli give rise to non-axisymmetric modes, similar to those encountered in elasticity. The very high moduli are admittedly unrealistic, but they are useful to check how and when non-axisymmetric modes can be obtained using the proposed plasticity formulation.

### 5.3.1. Free edges

Consider a rather large thickness $t = 7$ cm, a medium hardening modulus $h = 2 \times 10^{10}$ Pa, and prescribed stresses at the ends of the cylinder. Fig. 6 shows that the first bifurcated solution has three circumferential waves, which is in agreement with experimental and numerical results obtained by Goto and Zhang (1999). Fig. 7 depicts the distribution of the equivalent plastic strain $P$ in the cylinder when the nominal axial compressive stress is $\lambda = 1.617 \times 10^8$ Pa. The maximum value of 8.4% is reached at the junction of two circumferential waves which are localized at the ends of the cylinder.

Let us maintain the boundary conditions and reduce the thickness to $t = 5$ cm so as to see the influence of the radius-to-thickness ratio $R/t$ on the circumferential wave number. Fig. 8 shows that the wave number increases with the ratio $R/t$—now there are four circumferential waves with the new thickness—as was already observed in Goto and Zhang (1999) and Tvergaard (1983a). The bifurcation mode is localized at both ends and displays a series of periodic inward dimples.

Fig. 9 describes the load-displacement curves for $t = 5$ cm and $t = 7$ cm. In both cases, the (complete) plastication of the cylinder just precedes the bifurcation, which itself precedes the maximum stress value.

<table><caption>Table 6
Geometries used in Section 5.3</caption><thead><tr><td>Length l</td><td>Radius R</td><td>Thickness t</td></tr></thead><tbody><tr><td>5 m</td><td>1.5 m</td><td>7 cm, 5 cm, 3 cm, 2 cm, 12 mm, 6 mm according to the considered case</td></tr></tbody></table>

<table><caption>Table 7
Material properties used in Section 5.3</caption><thead><tr><td>Young's modulus E</td><td>Poisson's ratio v</td><td>Yield stress σ₀</td><td>Isotropic hardening modulus h</td></tr></thead><tbody><tr><td>2.1 × 10¹¹ Pa</td><td>0.3</td><td>4 × 10⁸ Pa</td><td>2 × 10¹⁰ Pa, 5 × 10¹⁰ Pa, 5 × 10⁸ Pa according to the considered case</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Non-axisymmetric bifurcation mode of an elastoplastic cylinder under axial compression (free edges, prescribed stresses, $t = 7$ cm).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Maximum equivalent plastic strain in the thickness of the cylinder under axial compression (free edges, prescribed stresses, $t = 7$ cm, nominal axial stress $\lambda = 1.617 \times 10^8$ Pa).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Non-axisymmetric bifurcation mode of an elastoplastic cylinder under axial compression (free edges, prescribed stresses, $t = 5$ cm).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Elastoplastic cylinder under axial compression (free edges, prescribed stresses): compressive stress versus axial shortening.</figcaption></figure>

**Remark.** About the mesh: it is necessary to mesh a quarter of the cylinder in order to get three circumferential waves as shown in Fig. 6. From now on only one eighth of the cylinder will be meshed so as to reduce the computational cost. Consequently, one necessarily obtains even circumferential wave numbers, as shown in Fig. 8.
### 5.3.2. Built-in edges

In the sequel, we consider the cases of built-in and simply supported edges, which are closer to experimental conditions. Here, the linear isotropic hardening modulus is given by the new value $h = 5.10^{10}$ Pa in order to obtain non-axisymmetric modes with dimples spreading over the entire surface of the cylinder (“diamond” modes) or limited in regions near the ends.

Under prescribed stresses, with built-in edges, two “diamond” modes are obtained with the thicknesses $t = 5$ cm and $t = 3$ cm, as shown in Fig. 10. These “diamond” modes display both axial and circumferential waves, whose numbers increase with the ratio $R/t$, as in the case of free edges.

Again, the load-displacement curves in Fig. 11 show that the bifurcation point occurs just before the limit load. It is noteworthy that the ratio between the two stresses at the bifurcation points is almost equal to the ratio of thicknesses, which reminds us the proportionality between the analytical stress in Eq. (58) and the thickness of the cylinder.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Primary modes (“diamond” modes) of an elastoplastic cylinder under axial compression (built-in edges, prescribed stresses).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Elastoplastic cylinder under axial compression (built-in edges, prescribed stresses): load-displacement curves for primary modes.</figcaption></figure>
**Remark.** Prescribed displacements give rise to similar results.

The secondary modes are obtained with lower hardening moduli, and only in the case of prescribed dis-
placements, because the localization of the axisymmetric mode seems to appear more easily in such conditions.
As in the case of the primary modes, the elastoplastic response is not homogeneous in the cylinder since the
edges are clamped. The cylinder is partly plastified before the bifurcation occurs, so that the bifurcation does
occur in plasticity. But unlike the primary mode case, here the plastification onset and the bifurcation are very
close. It is more difficult to computationally detect the bifurcation points for secondary modes. The step size
should be smaller than usual and some numerical devices have to be added in order to correctly deal with this
bifurcation type. Fig. 12 shows two “elephant foot” transitional secondary modes with $h = 5.10^8$ Pa and med-
ium thicknesses. As in the primary mode case, the circumferential wave number mainly depends on the radius-
to-thickness ratio.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Secondary modes (“elephant foot” modes) for an elastoplastic cylinder under axial compression (built-in edges, prescribed displacements).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Elastoplastic cylinder under axial compression (built-in edges, prescribed displacements): load-displacement curves for secondary modes.</figcaption></figure>
The load-displacement curves for secondary modes are shown in Fig. 13. It is worth mentioning the presence of a so-called snap-back (a limit point in displacement) in the case of the thinner shell.

### 5.3.3. Simply supported edges

Consider now a cylinder with simply supported edges and thickness $t = 2$ cm. In the case of prescribed stresses on the ends, a non-axisymmetric mode is obtained with dimples localized at the ends of the cylinder, as shown in Fig. 14a. As shown in Eq. (48a,b), the axial stress is prescribed on the edges through the control parameter $\lambda$. Here, since the numerical computations show that there can be some limit load, one has to use the arc-length method (see Section 4.2) to control the prescribed stress.

Fig. 14b shows that in the case of prescribed displacements on the ends, the deformation is more significant and there are more dimples in the middle of the cylinder than in the case of prescribed stresses. This fact could be accounted for as follows: since the lower and upper circumferences are constrained to remain plane under prescribed displacements, the global stiffness increases as compared with the case of prescribed stresses, thus increasing the dimples number too.

Fig. 15 shows that the load-displacement curves are very close before the critical point. Contrary to the case of prescribed stresses, the post-critical branch displays several successive limit points in the case of prescribed displacements, which must appear in the same time as additional dimples.

## 5.4. Summary

Table 8 summarizes the encountered bifurcation modes and the post-bifurcation behaviors versus the cylinder thickness and the boundary conditions. The major influence of the shell thickness on the axial and circumferential wave numbers has been highlighted for all encountered mode types.

(a) In the case of sufficiently thick shells and for some combinations of $R/t$ and stress-strain properties, the axisymmetric mode degenerates into a non-axisymmetric one after the limit load of the axisymmetric response and the second bifurcation corresponds to the so-called “elephant foot” mode. Note that in the cases of simply supported or built-in edges, the pre-critical deformations are found to be negligible and an axisymmetric mode is still possible.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Primary modes (“diamond” modes) of an elastoplastic cylinder under axial compression with simply supported edges.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Elastoplastic cylinder under axial compression (simply supported edges, $t = 2$ cm): compressive stress versus axial shortening.</figcaption></figure>

<table><caption>Table 8
Bifurcation modes and post-bifurcation behaviors of a cylinder under axial compression versus the thickness and the boundary conditions</caption><thead><tr><th>Cylinders</th><th>Simply supported edges</th><th>Built-in edges</th><th>Free edges</th></tr></thead><tbody><tr><td>Thick</td><td>Sinusoidal axisymmetric mode<br>Post-bif: localization at the ends<br>(Fig. 4)</td><td>Sinusoidal axisymmetric mode<br>Post-bif: localization somewhat far from the ends (Fig. 4)</td><td>Damped sinusoidal axisymmetric mode (Fig. 1)</td></tr><tr><td>Moderately thick</td><td>(Results between those for thick and thin cylinders)</td><td>Axisymmetric mode (+localization)<br>"Elephant foot" secondary mode (Fig. 12)</td><td>(Results between those for thick and thin cylinders)</td></tr><tr><td>Thin</td><td>"Diamond" mode (Fig. 14)</td><td>"Diamond" mode (Fig. 10)</td><td>Damped sinusoidal non-axisymmetric mode (Figs. 6 and 8)</td></tr></tbody></table>

(b) For thin shells, it has been found that there are no more axisymmetric modes (as if the two bifurcation points coalesce into one) and one gets straight to the so-called "diamond" mode, before the limit load of the "diamond" response.

The hardening modulus value has a great influence on the solution, as it may give rise to secondary bifurcation modes or not. Only quite high hardening modulus values (and the thicknesses considered) lead to interesting post-critical deformed shapes as presented in this work.

## 6. Conclusions

This work deals with the numerical solution of the elastoplastic thin shell buckling and post-buckling problem.

The proposed shell finite element has been validated by showing that the numerical results on the plastic buckling of axially compressed cylinders compare well with the analytical solutions derived from the plastic bifurcation theory.

The numerical computations have been carried out considering the influence of the geometry, the material properties and the boundary conditions. For the cylinder under axial compression, it has been shown that a large range of primary and secondary buckling modes can be computed in a unified way without using imperfections, and that the post-critical behavior can be described in an advanced deformation state, exhibiting the localization phenomenon as observed in experiments.
## Acknowledgements

The authors thank the reviewers for their instructive comments which have helped to clarify many statements in the manuscript.

## References

Bardi, F.C., Kyriakides, S., 2006. Plastic buckling of circular tubes under axial compression—part I: Experiments. International Journal of Mechanical Sciences 48, 830–841.

Bardi, F.C., Kyriakides, S., Yun, H.D., 2006. Plastic buckling of circular tubes under axial compression—part II: Analysis. International Journal of Mechanical Sciences 48, 842–854.

Batdorf, S.B., 1949. Theories of plastic buckling. Journal of Aeronautical Sciences 16, 405–408.

Batterman, S.C., 1965. Plastic buckling of axially compressed cylindrical shells. AIAA Journal 3 (2), 316–325.

Batterman, S.C., 1968. Free-edge plastic buckling of axially compressed cylindrical shells. Journal of Applied Mechanics 35, 73–79.

Bergan, P.G., Horrigmoe, G., Krakeland, B., Soreide, B., 1978. Solution techniques for non-linear finite element problem. International Journal for Numerical Methods in Engineering 12, 1677–1696.

Bijlaard, P.P., 1949. Theory and tests on the plastic stability of plates and shells. Journal of Aeronautical Sciences 16, 529–541.

Brank, B., Peric, D., Damjanic, F.B., 1997. On large deformations of thin elastoplastic shells: implementation of a finite rotation model for quadrilateral shell element. International Journal for Numerical Methods in Engineering 40, 689–726.

Bushnell, D. 1982. Plastic buckling, Pressure vessels and piping: Design technology. In: Zamrik, S.Y., Dietrich, D. (Eds.), A decade of Progress, Book No. G00213, 47–117.

Bushnell, D., 1985. Computerized Buckling Analysis of Shells. M. Nijhoff.

Christoffersen, J., Hutchinson, J.W., 1979. A class of phenomenological corner theories of plasticity. Journal of the Mechanics and Physics of Solids 27 (5–6), 465–487.

Crisfield, M.A., 1991. Non-linear finite element analysis of solids and structures. In: Essentials, 1. John Wiley & Sons.

Gellin, S., 1979. Effect of an axisymmetric imperfection on the plastic buckling of an axially compressed cylindrical shell. Journal of Applied Mechanics 46, 125–131.

Goto, Y., Zhang, C., 1999. Plastic buckling transition modes in moderately thick cylindrical shells. Journal of Engineering Mechanics 125 (4), 426–434.

Green, A.E., Naghdi, P.M., 1965. A general theory of an elastic–plastic continuum. Archive for Rational Mechanics and Analysis 18, 251–281.

Green, A.E., Naghdi, P.M., 1971. Some remarks on elastic–plastic deformation at finite strain. International Journal of Engineering Sciences 9, 1219–1229.

Halphen, B., Nguyen, Q.S., 1975. Sur les matériaux standard généralisés (On the generalized standard materials). Journal de Mécanique 14 (1), 39–63.

Hill, R., 1958. A general theory of uniqueness and stability in elastic–plastic solids. Journal of Mechanics and Physics of Solids 6, 236–249.

Hutchinson, J.W., 1974. Plastic buckling. Advances in Applied Mechanics 14, 67–144.

Kyriakides, S., Bardi, F.C., Paquette, J.A., 2005. Wrinkling of circular tubes under axial compression: effect of anisotropy. Journal of Applied Mechanics 72, 301–305.

Lam, W.F., Morley, C.T., 1992. Arc-length method for passing limit points in structural calculation. Journal of Structural Engineering 118 (1), 169–185.

Peek, R., 2000. Axisymmetric wrinkling of cylinders with finite strain. Journal of Engineering Mechanics 126 (5), 455–461.

Riks, E., 1979. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures 15, 524–551.

Riks, E., 1991. On formulations on path-following techniques for structural stability analysis. In: New Advances in Computational Structural Mechanics. Elsevier, 65–80.

Seydel, R., 1994. Practical Bifurcation and Stability Analysis. From Equilibrium to Chaos. Springer-Verlag.

Shanley, F.R., 1947. Inelastic column theory. International Journal of Aeronautical Sciences 14, 261–267.

Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Springer.

Tvergaard, V., 1983a. On the transition from a diamond mode to an axisymmetric mode of collapse in cylindrical shells. International Journal of Solids and Structures 19 (10), 845–856.

Tvergaard, V., 1983b. Plastic buckling of axially compressed circular cylindrical shells. Thin-Walled Structures 1 (2), 139–163.

Wagner, W., Wriggers, P., 1988. A simple method for the calculation of post-critical branches. Engineering and Computations 5, 103–109.

Wempner, G.A., 1971. Discrete approximations related to non-linear theories of solids. International Journal of Solids and Structures 7, 1581–1599.