Buckling of Thin-Walled Cylinders from Three Dimensional Nonlinear Elasticity

The famous bifurcation analysis performed by Flügge on compressed thin-walled cylinders is based on a series of simplifying assumptions, which allow to obtain the bifurcation landscape, together with explicit expressions for limit behaviours: surface instability, wrinkling, and Euler rod buckling. The most severe assumption introduced by Flügge is the use of an incremental constitutive equation, which does not follow from any nonlinear hyperelastic constitutive law. This is a strong limitation for the applicability of the theory, which becomes questionable when is utilized for a material characterized by a different constitutive equation, such as for instance a Mooney-Rivlin material. We re-derive the entire Flügge’s formulation, thus obtaining a framework where any constitutive equation fits. The use of two different nonlinear hyperelastic constitutive equations, referred to compressible materials, leads to incremental equations, which reduce to those derived by Flügge under suitable simplifications. His results are confirmed, together with all the limit equations, now rigorously obtained, and his theory is extended. This extension of the theory of buckling of thin shells allows for computationally efficient determination of bifurcation landscapes for nonlinear constitutive laws, which may for instance be used to model biomechanics of arteries, or soft pneumatic robot arms.


Introduction
Buckling of thin-walled cylinders subject to axial thrust represents one of the most famous problems in mechanics and a fascinating question in bifurcation theory.In fact it is well-known 1 that the critical load for buckling (calculated in a linearized context by Lorenz (1908), Timoshenko (1910), Southwell (1914), von Mises (1914), Flügge (1932), and Donnell (1933) and elegantly reported by Flügge (1962) and Yamaki (1984)) provides only an overestimation of the carrying capacity which can experimentally be measured on real cylinders.This overestimation was explained in terms of post-critical behavior in a number of celebrated works (among which, von Kármán and Tsien (1941), Koiter (1945), Hutchinson (1965), Hutchinson andKoiter (1970), andTsien (2012)).Fifty years later, the mechanics of thin shells remains a prosperous research topic (Lee et al. (2016), Jiménez et al. (2017), Elishakoff (2014), and Ning and Pellegrino (2016)), also embracing recent applications to nanotubes (Waters et al. (2005)) and soft materials, the latter developed as a key to understand biological systems, (Lin et al. (2018) and Steele (1999)) or towards mechanical applications, (Leclerc and Pellegrino (2020) and Shim et al. (2012)).
The bifurcation analysis performed by Flügge is based on a series of approximations, among which, the incremental constitutive equations do not follow from a finite strain formulation of any hyperelastic material.In particular, it is shown that the equations relate, through a fourth-order isotropic elastic tensor, the Oldroyd increment of the Kirchhoff stress to the incremental Eulerian strain.These equations, involving Lamé moduli λ and µ are certainly valuable in an approximate sense, but how this approximation may be tied to a rigorous theory of nonlinear elasticity remains unknown.
The focus of the present article is the incremental2 bifurcation analysis of an axially-loaded thinwalled cylinder, characterized by rigorously-determined, nonlinear hyperelastic constitutive equations.Our analysis generalizes and rationalizes the famous derivation performed by Flügge not only from the point of view of the constitutive equations, but also because it allows to either rigorously prove, or clearly elucidate other assumptions.In particular, the Flügge derivation is based on the smallness assumption for the thickness of the cylinder wall.This represents an approximation on the one hand and a simplification on the other.There are only three alternatives to circumvent this approximation, namely: (i.) a numerical approach (for instance through a finite element code), but numerical solutions are approximated and far from providing the deep insight and the generality intrinsic to a theoretical determination; (ii.) a direct approach from three-dimensional incremental elasticity (as for instance pursued by Wilkes (1955), Haughton and R. Ogden (1979), Bigoni and Gei (2001) and Chau (1995)), but the numerical solution of the bifurcation condition involved in this technique becomes awkward in the thin-walled limit; (iii.) a reduction (if possible) of the nonlinear elastic constitutive laws to a small-strain version, based on Lamé constants to be used in the Flügge equations, but in doing so, an unknown approximation is introduced.
The three above-mentioned alternatives are abandoned in this article (except for the 'direct approach' that will be used to validate the obtained results), in favor of a re-derivation of the buckling of a thin-walled cylinder, pursued from a different perspective.First, the incremental equilibrium equations are rigorously derived in terms of mean quantities, represented by generalized stresses (holding true regardless of the thickness of the cylinder), through a generalization of the approach introduced by Biot (1965) for rectangular plates.The incremental kinematics is postulated as a deduction from the deformation of a two-dimensional surface, again in analogy with the incremental kinematics of a plate.Our treatment of thin-walled cylinders allows the use of every nonlinear constitutive law.In particular, two different nonlinear elastic constitutive equations describing compressible neo-Hookean materials (Pence and Gou, 2015) are rigorously used.While the linearized kinematics adopted coincides with that used by Flügge, the incremental equilibrium equations derived in this article are different from Flügge's corresponding equations, but are shown to reduce to the latter by invoking smallness of the cylinder wall thickness.The equations obtained for the incremental deformation of prestressed thin cylindrical shells (Sections 2-4) are general and can be used for different purposes, so that the ensuing bifurcation analysis represents only an example of application, while other problems can be pursued, such as for instance, the torsional buckling.When compared (Section 5), the bifurcation landscapes obtained from our formulation and that given by Flügge are shown to be almost coincident and perfectly consistent with results obtained through the 'direct approach', where the fully three-dimensional problem is solved (which is also a new result presented here in Section 7).Finally, the following formulas are rigorously obtained as limits of our approach: (i.) the surface instability, in the short longitudinal wavelength limit; (ii.) the wrinkling, occurring as axial buckling of a mildly long cylindrical shell, characterized by the well-known formula obtained by Flügge, (iii.) the Euler rod buckling for a long cylindrical shell (Section 6).
The re-derivation of the Flügge formulation within a three-dimensional finite elasticity context, including the calculations of the bifurcation loads and the determination of the famous formula for buckling of a mildly long cylindrical shell, is important from two different perspectives.First, the validity of the Flügge theory, considered a reference in the field, is now confirmed.Second, the new derivation is applicable to soft materials, characterized in the framework of nonlinear elasticity by constitutive equations different from those used by Flügge.Therefore, the determination of the buckling stress is now possible for a cylindrical shell made up of an Ogden or a neo-Hookean compressible elastic material (Levinson and Burgess 1971;R. W. Ogden 1972), or for an artery obeying the Holzapfel et al. (2000) constitutive law.

Incremental field equations in terms of generalized stresses
The undeformed stress-free configuration is described by means of cylindrical coordinates (r 0 , θ 0 , z 0 ), being the z 0 -axis aligned parallel to the axis of revolution of the shell.Along its fundamental path before bifurcation, the shell is assumed to undergo a homogeneous, axisymmetric compression in its longitudinal direction z 0 , preserving the circular cylindrical geometry.A uniaxial stress is generated in the form where K = JT represents the Kirchhoff stress tensor, with J = det F, being F the deformation gradient and T the Cauchy stress, while G = e z ⊗ e z (e z is the unit vector singling out the z-axis).The current configuration is described through coordinates (r, θ, z) by means of the principal stretches {λ r , λ θ , λ z } as r = λ r r 0 , θ = θ 0 , z = λ z z 0 , with λ r = λ θ following from axial symmetry.Therefore, the deformation gradient and the left Cauchy-Green deformation tensor read as F = diag{λ r ; λ θ ; λ z } and B = FF ⊺ = diag{λ 2 r ; λ 2 θ ; λ 2 z }, respectively.The incremental equilibrium equations are derived, governing the bifurcations of a cylindrical shell of current length l, external radius r e and internal radius r i .The cylinder, whose thickness is denoted by t = r e − r i , is not assumed to be thin for the moment, therefore all results presented in this Section are rigorous in terms of mean values of the incremental field quantities.The geometrical descriptors adopted here are the mid-radius a = (r e + r i ) 2, defining the 'mid-surface' of the shell, and the so-called reduced radius r = r − a.A standard notation is used, where bold capital and lower case letters denote tensors and vectors, respectively.
Adopting the relative Lagrangean description, with the current configuration assumed as reference, such that F = I, and neglecting the body forces, the incremental equilibrium of a pre-stressed solid is expressed through Ṡ, the increment of the first Piola-Kirchhoff stress tensor S, as div Ṡ = 0. (2.2) The cylindrical shell is subject to traction-free surface boundary conditions on its lateral surface, so that Ṡir = 0 as r = ±t 2 (i = r, θ, z). (2. 3) The increment of the Kirchhoff stress K can be related to Ṡ through equation S = KF − ⊺ , namely (2.4) where L = grad v is the gradient of the incremental displacement field v.In a relative Lagrangean description, equation (2.4) becomes Ṡ = (tr L) T + Ṫ − TL ⊺ . (2.5) Introducing the uniaxial pre-stress, Eq. (2.1), into Eq.(2.5), the following relations between the components of the incremental first Piola-Kirchhoff stress tensor Ṡ are derived: (2.6)

Generalized stresses
In the shell theory, it is common to introduce the so-called 'generalized stresses', namely, stress resultants per unit length referred to the mid-surface of the shell.For a cylinder of current uniform wall thickness t = λ r t 0 , the following definitions are adopted for the increments of forces and moments: where the subscript ⋅ stands for r, θ, or z in turn, while 'stress ⋅z ' and 'stress ⋅θ ' represent the ⋅z and the ⋅θ component of a generic Eulerian stress measure.The superimposed ⋆ identifies a suitable increment, in particular here symbols ⋅ and ○ are used to denote material time derivative and Oldroyd derivative, respectively.The factor 1 + r a in ⋆ n ⋅z and ⋆ m ⋅z is the consequence of the integration over a circular sector.
The following generalized stresses play a role hereafter: ( Ṡrθ + Ṡθr ) r dr = 0. (2.9) The derivatives of the generalized moments ṁθθ and ṁθz according to Eqs. (2.7) can easily be recognized in the above equation, while an integration by parts allows to transform the first term as so that, exploiting Eq. (2.6) 1 , Eq. (2.9) becomes ṁθθ,θ + a ṁθz,z (2.10) The same procedure can be applied to Eq. (2.8) 3 after multiplication by r and subsequent integration to generate the next rotational equilibrium equation (2.11) where P = K zz t represents the pre-stress load per unit length along the mid-circular surface, multiplied by J. From a mechanical point of view, Eqs.(2.10) and (2.11) enforce the equilibrium of moments about the z-and θ-axes, respectively.The three translational equilibrium equations for the generalized stresses are obtained in a similar vein, through a direct through-thickness integration of Eqs.(2.8) with an integration by parts (2.12) Enforcing the boundary conditions, Eq. (2.3), on Eqs.(2.10)-(2.12), the full system of equilibrium equations is finally obtained (2.13) A substitution of Eq. (2.13) 4 and Eq.(2.13) 5 into Eq.(2.13) 1 and Eq.(2.13) 2 allows to remove the shear forces, thus leading to the following equations: (2.14)

Incremental equilibrium equations: spatial formulation
In a relative Lagrangean description the incremental equilibrium equations (2.14) can equivalently be expressed by means of a new set of generalized stresses, based on the Oldroyd increment (Oldroyd 1950) of the Kirchhoff stress K, namely K = Ṡ − LK . (2.15) The traction-free incremental boundary conditions (2.3) can be re-expressed through K as and a new set of generalized stresses is obtained from the initial definition, Eqs.(2.7).In fact, by introducing the components of K given by Eqs.(2.15), the first three Eqs.(2.14) are given the following 'spatial' format: (2.17)

Rotational equilibrium about axis r
A sixth incremental equilibrium equation expressing the rotational equilibrium about axis r can be obtained from a through-thickness integration of Eq. (2.6) 3 after multiplication by (1 + r a) The introduction of the generalized stresses in Eq. (2.7) leads to the expression in material formulation while the spatial version in terms of Oldroyd increments reads Note that all equations obtained until now, in particular Eqs.(2.14), (2.17), (2.19), and (2.20) do not involve any approximation and thus are rigorous.

The Flügge approximation
As already mentioned, all equations derived so far, to be used in the following elaboration, are exact.Interestingly, the corresponding equations provided by Flügge (1962) can be recovered as an approximation of Eqs.(2.17), when the assumption is introduced that the cylinder wall thickness t is small.In fact, a Taylor series expansion allows to show that and therefore the equations introduced by Flügge (1962) are recovered: (2.21) In addition to the above equations, Flügge used also Eq. (2.20), albeit he did never explicitly mention the use of either the Oldroyd increment or the Kirchhoff stress measure.

Incremental deformations of a prestressed shell
As a premise, the Euler-Bernoulli beam theory is briefly discussed on the basis of the standard assumptions (Love 1906).The kinematics of a beam in a plane is described through the displacement ū(x 01 ) of a generic point lying on its centroidal axis, singled out by the material coordinate x 01 along the straight reference configuration, x 0 = x 01 e 1 .Assuming that the centroidal axis behaves as the Euler's elastica corresponding to the evolution of an extensible line, the unit normal vector n (counterclockwise rotated π 2 with respect to the tangent Bigoni (2012) and Bigoni (2019)) For any point of the beam in its spatial configuration, x = x 1 e 1 + x 2 e 2 , having x 0 = x 01 e 1 + x 02 e 2 as material counterpart with x 02 ∈ [−t 2, +t 2], the following displacement field is postulated: If the derivatives of the displacement components (3.2) are negligible compared to unity (i.e.u 1,1 and u 2,1 ≪1), the linearized kinematics of the Euler-Bernoulli beam is recovered, namely, The kinematics of the incremental deformations in a prestressed cylindrical shell is illustrated as an extension of the development outlined above for the beam, following the standard assumptions discussed, among others, by Love (1906), Flügge (1932), Podio-Guidugli (1989), Steigmann and R. W. Ogden (2014), andGeymonat et al. (2007).In a cylindrical coordinate system, the prestressed shell configuration and its evolution after superposition of an incremental deformation are described through the geometry of the midsurface (Malvern (1969), R. W. Ogden (1984), and Chapelle and Bathe ( 2011)), respectively as where a is the radius of the prestressed cylindrical midsurface, while v r (θ, z), v θ (θ, z) and v z (θ, z) represent its incremental displacement components.The unit normal to the deformed surface is defined as where ⋅ represents the norm of its vector argument, while the derivatives read as It is useful to consider two new vectors parallel to x ′ ,θ and x ′ ,z , respectively, (3.7) Up to the leading-order, assuming the incremental displacement components v r and v θ to be small (negligible if compared to radius a), and the incremental displacement gradient to be negligible with respect to unity, the following approximations can be introduced (3.8) while the unit normal to the cylindrical surface n follows as (3.9) Paralleling the beam theory assumption, Eqn.(3.2), the incremental kinematics of a cylindrical shell is represented in the form (Chapelle and Bathe 2011) (3.10) On the basis of the above-described linearized kinematics, the gradient of the incremental displacement becomes so that the components of the Eulerian strain increment tensor D = (L + L ⊺ ) 2 are (3.12)

Two constitutive equations for compressible hyperelasticity
Two hyperelastic material models, both isotropic in their undeformed state, are considered, for which the strain energy functions are provided by Pence and Gou (2015, their Eqs. (2.11) and (2.12)).
Adopting the same notation proposed by those authors, the strain energy functions W a and W b are adopted, namely, and where I 1 (B) = tr B and I 3 (B) = det B, while µ and κ represent the shear and bulk moduli of the material in its unstressed state, related to the Young modulus E and Poisson's ratio ν through the usual formulae, namely, The strain energy function (4.1) is a special form of the general Blatz-Ko material model, in contrast with the strain energy function (4.2), which allows instead a separation between the pure volumetric effects and other contributions from the deformation.Both the models describe compressible neo-Hookean materials and satisfy, in the undeformed state, the stress-free condition, as well as the consistency with the classical linearized elasticity theory.Therefore, with reference to the generic strain energy function W , the following conditions hold true where the derivatives W ,i = ∂W (I 1 , I 3 ) ∂I i are to be evaluated for I 1 = I 2 = 3 and I 3 = 1 (Horgan and Saccomandi 2004).
The Cauchy stress, in general defined according to assumes for the strain energy (4.1) the expression while, for the strain energy (4.2), it reads as Through the relative Lagrangean description, in which the current configuration is assumed as reference, the Oldroyd increment of the Kirchhoff stress turns out to be related to the strain energy density of a hyperelastic material W as (Bigoni 2012) where E (2) denotes the Green-Langrange strain tensor, while the tensor product ⊠ is defined as Inserting the form (4.1) for the strain energy function W a into Eq.(4.7), the following expression for the elastic fourth-order tensor H is derived where S is the fourth-order symmetrizer tensor, leaving D unchanged because of symmetry, namely, Therefore the Oldroyd increment of the Kirchhoff stress (4.5) for the model with strain energy (4.1) becomes while, paralleling the procedure for the model with strain energy W b , Eq. (4.2), the following expression is obtained: It is noteworthy to point out that the constitutive equation used by Flügge can be recovered from both Eqs.(4.8) and (4.9), assuming the pre-stressed and unstressed configurations to be coincident, so that F = I, a condition leading to K = (κ − 2 3 µ) (tr D) I + 2µD, (4.10) which represents the incremental law used by Flügge.

Axisymmetric pre-stress
The axisymmetric ground-state assumption prescribes the coincidence of the radial and circumferential stretches, λ r = λ θ , as well as the vanishing of the radial and circumferential stress components, T rr = T θθ = 0. Therefore, with regard to strain energy function W a , the following condition is obtained from Eq. (4.5) for the components of the Cauchy stress in the trivial configuration: Solving Eq. (4.11) for the radial stretch λ r yields Note that both δ and λ r are real for ν ∈ [0, 0.5].The stress tensor in Eq. (4.5) can be simplified by means of Eq. (4.12), so that its only nonzero component turns out to be . (4.13) A substitution of Eq. (4.12) into Eq.(4.8) yields the following expressions for the diagonal (denoted assuming the repeated indices i not to be summed over) and the out-of-diagonal components of the Oldroyd increment of the Kirchhoff stress, tensor K: Similarly to Eq. (4.11), the condition of axisymmetric pre-stress for a material admitting the strain energy function W b , Eq. (4.2), can be written as which can be solved numerically to compute λ r as a function of the pre-stretch λ z for ν ∈ [0, 0.5).The corresponding axial stress T bzz , as well as the components of the Oldroyd increment of the Kirchhoff stress Kb ij are finally evaluated by means of Eq. (4.6) and Eq.(4.9), respectively.Noteworthy, in the incompressible limit, ν → 0.5, the radial stretch tends to the incompressibility constraint λ r = λ −1 2 z with dimensionless axial stress T zz µ = (λ 3 z − 1) λ z for both materials with the strain energy functions W a and W b , Eq. (4.1) and Eq.(4.2).
The radial stretch λ r and the dimensionless axial stress T zz µ are reported in Fig. 1, for the two models characterized by the strain energies W a (blue lines) and W b (red lines), as functions of the axial stretch λ z .For the strain energy function W a , Eqs. (4.12) and (4.13) have been used, while for the strain energy function W b , Eq. (4.15) has to be numerically solved.Three values of ν are reported The curves demonstrate the high non-linearity of the models and the differences in the mechanical response to stretch.Note that when ν = 0, the radial stretch is constant and equal to unity, λ r = 1, for the strain energy W a , while for W b , λ r remains close to 1 for values of λ z > 0.7.

Incremental plane stress assumption
For cylinders having 'sufficiently' thin walls, the assumption of plane stress becomes reasonable and is hereafter extended to the bifurcation state as well, namely, Ṡrr = 0 , ∀r ∈ [−t 2; t 2] . (4.16) Recognizing that Krr = Ṡrr as a result of the assumed structure of the pre-stress in Eq. (2.1), together with Eq. (2.15), the enforcement of Eq. (4.16) for the material with the strain energy function W a in Eq. (4.1), yields to be further simplified through the introduction of Eq. (4.12) as Under the constraint represented by Eq. (4.18), the incremental constitutive equations (4.14) assume the following expression: For the material with strain energy function W b in Eq. (4.2), the fulfillment of the plane stress requirement, Eq. (4.16), instead of Eq. (4.18), leads to where Finally, the substitution of Eq. (4.20) into Eq.(4.9), after introducing the implicit relation λ r (λ z ) represented in Fig. 1 (a) that aims to satisfy Eq. (4.15), allows to determine the components of tensor Kb , whose expression remains in implicit form for the model with strain energy W b .

Bifurcation of an axially-compressed thin-walled cylinder
The bifurcation problem for an axially-compressed thin-walled cylinder is set up on the basis of the kinematical conditions (3.12), the equilibrium equations (2.17), expressed in terms of generalized incremental stresses, and the constitutive relations • Eq. (4.19), for the material obeying the strain energy function W a , • Eq. (4.9) together with Eq. (4.20), and the implicit relation λ r (λ z ) satisfying Eq. (4.15), for the material obeying the strain energy function W b .
The pre-stress load per unit length P in equations (2.17) can be evaluated for the two materials by means of Eq. (4.5) and Eq.(4.6), respectively.
In the following of this article, explicit calculations will be presented with reference to the constitutive law following from the strain energy function W a , Eq. (4.1), with the index a omitted (for the sake of conciseness).The analogous calculations we have performed for the function W b in Eq. (4.2) are not reported here.Final computations of the bifurcation solution and asymptotic derivations of limit loads will be presented for both models.
As standard in the incremental bifurcation analysis of elastic solids (Hill and Hutchinson 1975), the following ansatz is introduced for the incremental displacements at bifurcation, corresponding to a free sliding condition along perfectly smooth rigid constraints on the upper (z = l) and lower (z = 0) faces: (5.1) where n = 0, 1, 2, ... and η = mπa l (m = 1, 2, ...) represent the circumferential and the longitudinal wave-numbers, respectively, singling out the bifurcation mode, while the amplitudes are collected in the vector c = {c 1 , c 2 , c 3 } T .The incremental displacement field, Eq. (5.1), constant throughout the thickness of the shell, enforces the conditions of null incremental force nθz and moment mθz at the ends z = 0 and z = l.In Flügge (1973) the boundary conditions for the lower and upper ends were modelled as simple supports, therefore preventing radial and circumferential incremental displacements, while no restrictions were imposed on the axial incremental displacement.However, both the boundary conditions assumed by us and by Flügge lead to the same bifurcation conditions.Through the introduction of Eqs.(5.1) into the kinematical conditions (3.12) and the substitution into the constitutive relations (4.19), the final form of the three incremental equilibrium equations (2.17) is obtained, with the generalized stresses defined according to Eqs. (2.7).The bifurcation condition is eventually expressed in the standard form as M c = 0, where matrix M is a function of the axial stretch λ z (while λ r is replaced through Eq. (4.12)), the dimensionless thickness of the shell τ = t a, the material parameter ν and the wave-numbers n and η.Bifurcation occurs when the coefficient matrix is singular, det M = 0, (5.2) a condition that allows to define the critical stretch λ z for bifurcation (and therefore the corresponding dimensionless axial compressive load p z = −P D, with D = Et (1 − ν 2 ) representing the extensional stiffness of the shell), as a function of the geometrical variable τ , the material parameter ν and the wave-numbers n and η.Fig. 2 shows the buckling diagram obtained for ν = 0.3 and r e r i = 1.05, so that τ = 0.0488 (note that both the radii ratio r e r i and the dimensionless thickness τ remain constant during the prebifurcation deformation, while the cylinder deforms maintaining its shape).The critical axial stretch is plotted as a function of the longitudinal wave-number η for different values of the circumferential wave-number n.The critical modes are illustrated as continuous lines, while the dashed lines represent the modes corresponding to high axial stresses that cannot be reached when the load is continuously increased from zero.As expected, for small values of η, corresponding to very slender cylinders, the mode representative of the Euler buckling, characterized by n = 1, becomes dominant.
A selection of the bifurcation eigenmodes for a thin shell, corresponding to different values of circumferential and longitudinal wave-numbers n and m, is displayed in Fig. 3, where the colours highlight the peculiar bulges of the buckled shell geometry.
Critical envelopes of the intersecting buckling curves are shown in Fig. 4, for different ratios r e r i and various circumferential wave numbers n.Depending on the ratio l (ma), the bi-logarithmic plot shown in Fig. 4 highlights the sequence of three different behaviour ranges, found by Flügge and now recovered for two exact models of compressible elasticity: • Cylinders with very small curvature (region on the left), tending to behave as plates, therefore the bifurcation condition approaches the plate buckling.Fig. 4 highlights how the bifurcation solution pertaining to a thin plate (denoted by the letter S in the figure), tends to progressively dissociate from the bifurcation solution for a thin-walled cylinder at increasing cylinder wall thicknesses.This analysis will be addressed in Sect.6.1; • Moderately long cylinders (intermediate region) present an almost constant buckling load, independent of both the circumferential and longitudinal wave-numbers.This load is denoted in Fig. 4 by the letter W and analyzed in Sect.6.2; • Cylinders with high slenderness (on the right) approach the Euler buckling solution, denoted in Fig. 4 by the letter E. A detailed investigation of this case is presented in Sect.6.3.
The results presented above are based on a large strain approach with a constitutive equation assuming the strain energy W a , Eq. (4.1).We have obtained similar results with the strain energy W b , Eq. (4.2), not reported here for conciseness.Both cases are different from the small strain analysis performed by Flügge, which is based on a constitutive equation not following from a potential.Nevertheless, results in terms of critical loads for bifurcation turn out to be only marginally dependent on the constitutive equations, because bifurcation occurs at low stretch.Therefore, a comparison between the approach pursued in this paper and the solution obtained by Flügge shows almost coincident results; the comparison is not reported here as the curves are scarcely distinguishable.
The accuracy of the current 2D approach (developed on the basis of two models provided by Pence-Gou and presented in §4) will definitely be assessed though a comparison with the 3D full-field solution for bifurcation on the basis of the constitutive model with strain energy W a , Eq. (4.1) (Fig. 5 in Sect.7).

Limiting cases via asymptotic analysis
Three crucial limiting cases are analyzed in this Section.The well-known solutions for cylinders with a very small curvature and for moderately long cylinders are rigorously derived from the finite elasticity approach developed in this article on the basis of both the constitutive models, Eqs (4.1) and (4.2).The problem of an Euler rod consisting in a hollow cylindrical shaft is finally addressed.In all cases the asymptotic solutions are obtained for both the elasticity models considered here.Again, for the sake of conciseness, all the results will be presented only for the material whose strain energy function is W a , Eq. (4.1).At this point it is convenient to introduce the relationships expressing the push-forward operation

Cylinders of very small curvature: surface instability
If the reference geometry of the shell is altered, increasing the radius a 0 , while keeping constant both the length l 0 and the thickness t 0 , a hollow cylinder of very small curvature is generated.The latter exhibits the surface instability of a plane plate strip with two free and two constrained opposite edges (endowed with simple supports, or, equivalently, sliding clamps), subject to an in-plane dead load.The bifurcation solution for such plate strip is known, see Timoshenko and Gere (1961) and Flügge (1973), where ) the flexural and extensional stiffnesses of the shell in its reference configuration (note that in current configuration, k = K (D a 2 ) = k 0 ).From the analysis of the critical pairs {λ z , η} obtained for the constitutive law Eq.(4.1), it turns out that, as recognized by Flügge, at high values of the longitudinal wave-number η, the dimensionless critical load for the plate strip p z,S , Eq. (6.2), approximates the curves corresponding to n = 0 in the dimensionless bifurcation load envelopes shown in Fig. 4 for thin shells.In order to capture this limit behaviour, a Taylor-series expansion in λ z , truncated at the linear term about λ z = 1, is introduced into the bifurcation condition (5.2), with matrix M evaluated at n = 0.The following critical stretch is obtained and neglecting the terms becoming inessential at large values of η, it can be further simplified to where At large longitudinal wave-numbers η, Eq. ( 6.4), suitable for thin shells, which are characterized by small axial deformation before bifurcation, allows to compute the leading order approximation for the dimensionless pressure p z at bifurcation.An additional third-order series expansion around τ = 0 leads to The above detailed procedure, based on an approximation of the bifurcation condition truncated at linear order in λ z , was repeated assuming an expansion up to the third order, which led to a much more cumbersome equation with respect to Eq. (6.4), but yielded precisely the same result, Eq. (6.5).
To allow a comparison with the plate strip solution, Eq. (6.2), the asymptotic solution above, expressed in terms of current variables, as usual in bifurcation analysis, is to be restated in terms of reference variables, thus an approximated explicit version of Eq. (6.1) 2 is sought.This equation is conveniently restated as η 2 −η 2 0 λ 2 r λ 2 z = 0, which turns out to involve only η, η 0 , ν, τ 0 upon introducing Eqs.(4.12) and (6.4) for λ r and λ z .Finally, the development of the latter condition into a Taylor series around τ 0 = 0 up to the order 3, yields a bi-quadratic equation in η, whose solution gives the following approximated relationship, such that, η → η 0 as τ 0 → 0. Noteworthy, Eq. (6.6) turns out to be valid for both the material models characterized by the strain energies W a and W b .Considering Eq. (6.5), with the variables η and τ replaced by η 0 and τ 0 through Eqs.(6.6) and (6.1) 1 , respectively, a final Taylor series expansion around τ 0 = 0 up to the third order, leads to therefore, for large longitudinal wave-numbers η 0 , Eq. (6.2) is recovered asymptotically.The dimensionless critical pressure for the plate strip, Eq. (6.2), is superposed as a straight dashed line in the bi-logarithmic plot reported in Fig. 4 for different values of τ 0 .The conclusion is that at large values of η, the plate theory provides a good approximation to the critical load of thin-walled cylinders.

Medium length cylinders: wrinkling
As highlighted by Timoshenko and Gere (1961), experiments show that thin cylindrical shells under compression usually buckle into short longitudinal waves, at a large longitudinal wave-number η.The bifurcation diagrams reported in Fig. 4 display an intermediate region where the buckling loads are almost independent of the values of both wave-numbers n and η.This region, for mildly long shells, corresponds to the so-called 'wrinkling' (Zhao et al. (2014)), for which Flügge (1973) derived the critical load This classic solution can be rigorously recovered within the developed framework, by seeking the bifurcation condition as a minimum of the dimensionless axial pressure p z with respect to variable η.This corresponds to the stationarity of the bifurcation axial stretch evaluated for the mode n = 0, Eq. ( 6.3), leading to five solutions.Among these, one is trivial, two are purely imaginary conjugated roots and two are real with opposite signs.From the latter pair, the positive real root is selected, , (6.9) where, assuming ǫ = 17ν 2 − 20ν + 5, Equation (6.9) is now introduced into Eq.( 6.3) to evaluate the minimum axial stretch for the mode n = 0.The latter is finally used to compute the corresponding load, which is further expanded about τ = 0 to obtain, at first-order in τ , (recalling Eq. (6.1) 1 ) exactly the Flügge Eq. (6.8).

Slender cylinders: Euler rod buckling
For a bar constrained with sliding clamps at both ends, assumed to be linearly elastic with Young modulus E, the Euler buckling load can be written in the form (6.10) where a 0 = r e0 − t 0 2, while α 0 = a 0 l 0 represents the stubbiness ratio, in other words the inverse of the slenderness ratio in the reference configuration (being α = a l its counterpart in the current configuration).Euler buckling, affecting slender shells, characterized by a small stubbiness ratio α 0 , corresponds to the anti-symmetric buckling mode with m = 1, n = 1.The Euler formula, Eq. (6.10), is recovered resorting to a perturbative technique (Golubitsky and Schaeffer 1985;Simmonds and Mann 1998;Kokotović et al. 1999;Holmes 2013) in the limit of vanishing longitudinal wave-number η.The approach followed by Goriely et al. 2008 for incompressible materials is generalized by expanding both the radial and longitudinal stretches λ r (α) and λ z (α) in power series about α = 0 up to the order M , (6.11) with coefficients λ ri and λ z i .The procedure is described here for the material with strain energy function W a in Eq. (4.1), but parallel computations have been performed for the strain energy W b , Eq. (4.2).The axisymmetry of pre-stress is enforced in Eq. (4.11) in an approximate form, developing for convenience λ 2 r λ z T arr in a Taylor-series about α = 0, with coefficients k j up to order M , ) .(6.12) The relation between the coefficients λ ri and λ z i is thus determined by requiring that the series in Eq. (6.12) vanishes at each order, so that the final approximation becomes To exemplify Eq. ( 6.13), when the order of approximation is M = 2, the following coefficients are computed: λ r1 = −νλ z1 and λ r2 = − 1 2 ν λ 2 z1 8ν 2 − 11ν + 2 + 2λ z2 .The approximations, Eqs.(6.13) and (6.11) 2 , are substituted in the bifurcation condition Eq. ( 5.2), with M computed for m = n = 1; note that an infinitely slender cylinder buckles at vanishing load, namely, {λ z , α} = {1, 0} represents a critical pair, therefore det M m=n=1 = 0 when α vanishes and thus λ r = λ z = 1.A further expansion into a Taylor series about α = 0 with coefficients d j up to order N , makes the buckling condition take the form In order to satisfy this condition at each order, all coefficients d j are enforced to vanish.This leads to a system of linear equations for the unknown parameters λ zi .As the coefficients d j (j = 1, 2) turn out to vanish, N = M + 2 is required to determine all the coefficients λ z i (i = 1, .., M ) in Eq. (6.11) 2 .It turns out that λ z i = 0 for all odd values of the index i = 1, .., M .Hence, the option N = 4 is sufficient to provide the asymptotic expansion up to the third-order of λ z (α), The critical axial stretch in Eq. (6.15) is defined with respect to the variables in the current configuration, and has to be related to the corresponding variables in the initial configuration to recover the critical load, Eq. (6.10).Therefore, the stubbiness ratio is expressed in terms of both current and initial variables as α = a l = α 0 λ r λ z , so that αλ z − α 0 λ r = 0. (6.16) The asymptotic expansions (6.13) and (6.11) 2 for λ r and λ z , respectively, plus a power series expansion about α 0 for the function α (with coefficients α k ), ) , (6.17) are introduced into Eq.(6.16).The obtained equation is solved at each order, thus obtaining the following expression (valid for P = 4) The longitudinal force resultant (positive when compressive) before bifurcation on the thin-walled tube can be finally computed as (6.19) so that inserting Eq. (4.5), and expanding the result in Taylor series about α 0 → 0 (slender columns), Eq. (6.19) becomes The buckling load asymptotically derived from finite elasticity under the assumption of plane stress, Eq. (6.20), can now be compared with the Euler buckling load, Eq. (6.10).It may be concluded that the two expressions for N z are identical at first-order in τ 0 , but differ at the third-order in τ 0 , because of a term depending on ν, so that the coincidence up to the fourth-order occurs only when ν = 0.This little discrepancy remains very small for ν ∈ [0, 0.5) when the dimensionless thickness τ 0 is small, i.e. for thin shells.In fact, the relative difference (N z − N z,Euler ) N z,Euler between the asymptotic approximation in Eq. (6.20) and the usual formula for Euler's critical load, Eq. (6.10), is an increasing function of ν and τ 0 , attaining a maximum of 0.42% as r e r i = 1.5 (τ 0 = 0.4); note that the latter is a value already far beyond the geometry of a thin shell.This is depicted for ν = 0.3 in Fig. 4, with the values of τ 0 spanning within the large interval [0, 0.4].
The asymptotic analysis has been repeated for the material with the strain energy function defined by Eq. (4.2), which allows for the separation of the volumetric effects.This analysis has yield the same asymptotic Euler buckling load up to order α 2 0 given by Eq. (6.20).It may be suggested that the 'discrepancy factor' multiplied by ν may be a consequence of both the incremental plane stress assumption and the simplified kinematics underlying the two-dimensional approach presented here.In fact, for the material with strain energy function W a in Eq. (4.1), the incremental plane stress assumption becomes exact for ν = 0 and the radial stretch becomes unity, λ r = 1, as depicted in Fig. 1 (a).For the material with strain energy W b , Eq. (4.2), the incremental plane stress assumption has an order of accuracy O(α 2 0 ), and the radial stretch at bifurcation is approximated by the unity, λ r = 1 + O(α 4 0 ).It should be noticed that in (Goriely et al. 2008), the classical Euler buckling formula is exactly recovered up to the order α 2 0 on the basis of a fully three-dimensional approach for an incompressible Mooney-Rivlin material.

3D bifurcation of a hollow (thick or not) cylinder
The fully three-dimensional solution for the bifurcation of a thick-walled cylinder made up of a hyperelastic material obeying the Pence-Gou model with strain energy W a , Eq. (4.1), is derived in this Section, following the procedure outlined by Chau (1995) (see also Chau (1993) and Chau and Choi (1998)) for a class of materials characterized by an incremental constitutive law in the form here the Zaremba-Jaumann (or corotational) rate of the Cauchy stress ▿ T = Ṡ − (tr D) T + T D − W T is adopted, as a function of D. The Pence-Gou model, Eq. (4.8), fits the incremental form (7.1), when the coefficients C ij are defined as Note that the incremental moduli (7.2) depend on the stretches λ r and λ z in the pre-bifurcation state and on the constitutive parameters κ and µ.
The three incremental equilibrium equations for the linearized bifurcation problem can be decoupled through the introduction of the two potentials Φ(r, θ, z) and Ψ(r, θ, z), such that 2) can thus be written as with ν 1 and ν 2 representing the roots of the characteristic equation The regimes can be classified, according to the nature of the roots ν 1 and ν 2 .The fulfillment of conditions B 2 − 4 A C > 0, A C > 0 and B > 0 define the elliptic-imaginary (EI) regime for the Pence-Gou material considered, where diffuse bifurcation modes are to be found (Chau 1995).
The following representation for diffuse eigenmodal bifurcations are introduced via the aboveintroduced potentials Φ(r, θ, z) = φ(r) cos (n θ) sin (η z) , Ψ(r, θ, z) = ψ(r) sin (n θ) cos (η z) , (7.6) where n and η maintain the same definitions as in Eqs.(5.1).This choice of the potential functions, automatically satisfying the boundary conditions of free sliding along the faces z = 0 and z = l, allows to write the equilibrium equations (7.4) in the form Figure 5: Comparison between the lower envelopes of dimensionless bifurcation loads p z for the buckling of cylindrical hollow cylinders, as a function of the longitudinal wave-number η in a bi-logarithmic plot.Results from the thin-shell approximation, developed for Pence-Gou compressible materials with strain energies W a and W b , are compared with the three-dimensional analysis performed for the Pence-Gou compressible materials with strain energy W a ; ν = 0.3 is adopted.
where H (1) n and H (2) n represent the Hankel functions of the first and second kind of order n, while I n and K n are the modified Bessel functions of the first and second kind of order n, with coefficients b i being complex numbers.
Enforcing the boundary conditions (2.3) of null tractions on both the inner and outer lateral surfaces of the pre-stressed cylinder, an eigenvalue problem in the form M b = 0 is obtained, with b = {b 1 , b 2 , b 3 , b 4 , b 5 , b 6 } T .Non-trivial solutions become possible when det M = 0.The latter condition only depends on the pre-bifurcation axial stretch λ z , the dimensionless thickness of the shell τ , the material parameter ν, as well as the circumferential and longitudinal wave-numbers n and η.For a given set of parameters ν, r e r i , n and η, the critical axial stretch can be found numerically.
A comparison is reported in Fig. 5 between the critical envelopes evaluated on the basis of the 3D approach and the thin-shell approximation.The 3D approach has been developed for a compressible Pence-Gou material with the strain energy W a , while results for the thin shell approximation are reported for both strain energies W a and W b .Geometry of the cylinder varies between very thin-, thin-and medium-walled.
It should be noted that the three-dimensional approach fully captures the nearly constant branch of the curve, corresponding to the asymptotic load derived by Flügge for medium length cylinders, Eq. (6.8).
The accuracy of the thin-shell approximation is evident from the comparison with the threedimensional solution described in the present Section.In particular, the critical modes characterized by small longitudinal wave-numbers η are neither altered by the chosen approach, nor by the consti-tutive model adopted, so that the curves reported in Fig. 5 are almost coincident within the most important part of the buckling landscape.On the contrary, for modes with small circumferential wave-numbers (in particular for n = 0, critical for large longitudinal wave-numbers η, Fig. 4) a noticeable difference between the curves becomes appreciable, becoming more evident when the thickness of the shell increases.This discrepancy is due to the fact that a surface bifurcation is approached and thus the thin-walled solution is no longer valid.
Of course the thin-shell approximation is much more efficient from the computational point of view (with CPU times for the single evaluation of a critical pair {λ r , η} according to the approximated approach getting as low as 1/300 of the times for the fully three-dimensional approach), however, the adopted hypothesis of incremental plane stress becomes unreliable as the shell thickness grows.

Conclusions
A complete re-derivation has been presented for the bifurcation of axially compressed thin-walled cylinders.The most important aspect of the new formulation is the independence from the constitutive equation used in the original formulation by Flügge, which does not stem from any strain potential and is now replaced by a generic nonlinear law of elasticity.Using two different hyperelastic constitutive laws, we have rigorously confirmed the results by Flügge, together with several limit formulae (for surface instability, wrinkling, and Euler rod buckling).The outlined approach allows now the precise and computationally efficient analysis of the bifurcation landscape for a thin-walled cylinder adopting all other related formulae for any constitutive law characterizing soft materials.

Figure 1 :
Figure 1: Uniaxial loading before bifurcation for a compressed cylinder, following from two the elastic models with the strain energies W a (blue lines) and W b (red lines).(a) The radial stretch λ r is determined as a function of the axial stretch λ z ; ν = {0, 0.3} are considered, together with the incompressibility limit, ν = 0.5.(b) The axial Cauchy stress T zz µ is determined as a function of the axial stretch λ z ; ν = 0.3 and ν = 0.5 are considered.

Figure 2 :
Figure 2: Critical stretch λ z of an axially-compressed thin-walled cylinder (r e r i = 1.05) made up of a Pence-Gou compressible material with strain energy function W a (ν = 0.3) as a function of the longitudinal wave-number η: the curves for different values of the circumferential wave-number n are denoted by n .Continuous lines represent the intersecting critical modes contributing to the buckling envelope, while the dashed lines correspond to modes arising at higher loads.The anti-symmetric mode labeled 1 represents Euler buckling (n = 1).

Figure 3 :
Figure 3: Different views for a selection of bifurcation eigenmodes for a shell with r e r i = 1.05.The material obeys the Pence-Gou model with strain energy W a and ν = 0.3.In particular, proceeding top-down, a surface instability mode, the Euler's buckling mode, and three different ovalization modes are shown.

Figure 4 :
Figure 4: Lower envelopes of the dimensionless load p z at bifurcation evaluated for a thin-walled cylinder made up of a Pence-Gou compressible material with strain energy W a (ν = 0.3) as a function of π η = l (ma) for different ratios r e r i between the external and internal radii of the cylinder (bi-logarithmic representation).The numbers adjacent to the curves indicate the critical circumferential modes of wave-number n, alternating along the envelopes for different values of π η.The dashed lines illustrate the asymptotic buckling loads for surface instability (S), wrinkling (W) and the Euler's column (E).