# Nonlinear Euler buckling

Alain Goriely, Rebecca Vandiver, Michel Destrade

## Abstract

The buckling of hyperelastic incompressible cylindrical tubes of arbitrary length and thickness under compressive axial load is considered within the framework of nonlinear elasticity. Analytical and numerical methods for bifurcation are developed using the exact solution of Wilkes for the linearized problem within the Stroh formalism. Using these methods, the range of validity of the Euler buckling formula and its first nonlinear corrections are obtained for third-order elasticity. The values of the geometric parameters (tube thickness and slenderness) where a transition between buckling and barrelling is observed are also identified.

Keywords:

## 1. Introduction

Under a large enough compressive axial load an elastic beam will buckle. This phenomenon known as elastic buckling or Euler buckling is one of the most celebrated instabilities of classical elasticity. The critical load for buckling was first derived by Euler in 1744 (Euler 1744, 1759; Oldfather et al. 1933) and further refined for higher modes by Lagrange in 1770 (Lagrange 1770; Timoshenko 1983). Both authors reached their conclusion on the basis of simple beam equations first derived by Bernoulli (Todhunter 1893; figure 1). Since then, Euler buckling has played a central role in the stability and mechanical properties of slender structures from nano- to macrostructures in physics, engineering, biochemistry and biology (Timoshenko & Gere 1961; Niklas 1992). Explicitly, the critical compressive axial load N that will lead to a buckling instability of a hinged–hinged isotropic homogeneous beam of length L is(1.1)where ‘π is the circumference of a circle whose diameter is one’ (Euler 1759); E is Young's modulus; and I is the second moment of area, which, in the case of a cylindrical shell of inner radius A and outer radius B, is .

Figure 1

Euler problem: (a) illustrations from Euler (1744). (b) Lagrange solutions (1770), modes 1, 2 and 3.

There are many different ways to obtain this critical value and infinite variations on the theme. If the beam is seen as a long slender structure, the one-dimensional theory of beams, elastica, or Kirchhoff rods, can be used successfully to capture the instability, either by bifurcation analysis, energy argument (Timoshenko & Gere 1961) or directly from the exact solution, which in the case of rods can be written in terms of elliptic integrals (Nizette & Goriely 1999). The one-dimensional theory can be used with a variety of boundary conditions, it is particularly easy to explain and generalize and it can be used for large geometric deflections of the axis (Antman 1995). However, since material cross sections initially perpendicular to the axis remain undeformed and perpendicular to the tangent vector, no information on the elastic deformation around the central curve can be obtained. In particular, other modes of instability such as barrelling cannot be obtained. Here, by barrelling, we refer to axisymmetric deformation modes of a cylinder or a cylindrical shell. These modes will typically occur for sufficiently stubby structures.

The two-dimensional theory of shells can be used when the thickness of the cylindrical shell is small enough. Then, the stability analysis of shell equations such as the Donnell–von Kármán equations leads to detailed information on symmetric instability modes, their localization and selection (Hunt et al. 2003). However, the theory cannot be directly applied to obtain information on the buckling instability (asymmetric buckling mode).

The three-dimensional theory of nonlinear elasticity provides, in principle, a complete and exact description of the motion of each material point of a body under loads. However, due to the mathematical complexity of the governing equations, most problems cannot be explicitly solved. In the case of long slender structures under loads, the buckling instability can be captured by assuming that the object is either a rectangular beam (Biot 1962; Levinson 1968; Nowinski 1969) or a cylindrical shell under axial load. Using the theory of incremental deformations around a large deformation-stressed state, the buckling instability can be recovered by a bifurcation argument, usually referred to, in the nonlinear elasticity theory, as small-on-large, or incremental, theory. By taking the proper asymptotic limit for long slender structures, the Euler criterion can then be recovered. In comparison to the one- and two-dimensional theories, this computation is rather cumbersome as it is based on non-trivial tensorial calculations, but it contains much information about the instability and the unstable modes selected in the bifurcation process.

Here, we are concerned with the case of a cylindrical shell under axial load. This problem was first addressed in the framework of nonlinear elasticity in a remarkable 1955 article by Wilkes who showed that the linearized system around a finite axial strain can be solved exactly in terms of Bessel functions. While Wilkes only analysed the first axisymmetric mode (n=0, see below), he noted in his conclusion that the asymmetric mode (n=1) corresponds to the Euler strut and doing so, opened the door to further investigation by Fosdick & Shield (1963), who recovered Euler's criterion asymptotically from the solution of Wilkes. These initial results constitute the basis for much of the modern theory of elastic stability of cylinders within the framework of three-dimensional nonlinear elasticity (Haughton & Ogden 1979a,,b; Simpson & Spector 1984; Duka et al. 1993; Pan & Beatty 1997; Bigoni & Gei 2001; Dorfmann & Haughton 2006). The experimental verification of Euler's criterion was considered by Southwell (1932) and by Beatty & Hook (1968).

The purpose of this article is threefold. First, we revisit the problem of the stability of an incompressible cylindrical shell under axial load using the Stroh formalism (Stroh 1962) and, based on the solution of Wilkes, we derive a new and compact formulation of the bifurcation criterion that can be used efficiently for numerical approximation of the bifurcation curves for all modes. Second, we use this formulation to obtain nonlinear corrections of Euler's criterion for arbitrary shell thickness and third-order elasticity. Third, we consider the problem of determining the critical aspect ratio where there is a transition between buckling and barrelling.

## 2. Large deformation

We consider a hyperelastic homogeneous incompressible cylindrical tube with isotropic cross sections of initial inner radius A, outer radius B and length L, subjected to a uniaxial constant strain λ3 and deformed into a shorter tube with current dimensions a, b and l. The deformation bringing a point at (R, Θ, Z), in cylindrical coordinates in the initial configuration, to (r, θ, z) in the current configuration is(2.1)where and λ3=l/L. Since the material is isotropic in the cross sections, the physical components of the corresponding deformation gradient F are(2.2)showing that the principal stretches are the constants λ1, λ2=λ1, λ3; and that the pre-strain is homogeneous. Owing to incompressibility, det F=1, so that(2.3)The principal Cauchy stresses required to maintain the pre-strain are (Ogden 1984)(2.4)(no sum), where p is a Lagrange multiplier introduced by the internal constraint of incompressibility and W is the strain energy density (a symmetric function of the principal stretches). In our case, σ2=σ1 because λ2=λ1. Also, σ1=0 because the inner and outer faces of the tube are free of traction. It follows that:(2.5)where Wi≡∂W/∂λi, and we conclude that the principal Cauchy stresses are constant.

## 3. Instability

To perform a bifurcation analysis, we take the view that the existence of small deformation solutions in the neighbourhood of the large pre-strain signals the onset of instability (Biot 1965).

### (a) Governing equations

The incremental equations of equilibrium and incompressibility can be written as (Ogden 1984)(3.1)where is the incremental nominal stress tensor and u is the infinitesimal mechanical displacement. They are linked by(3.2)where is the increment in the Lagrange multiplier p and 0 is the fourth-order tensor of instantaneous elastic moduli. This tensor is similar to the stiffness tensor of linear anisotropic elasticity, with the differences that it possesses only the major symmetries, not the minor ones, and that it reflects strain-induced anisotropy instead of intrinsic anisotropy. Its explicit non-zero components in a coordinate system aligned with the principal axes are (Ogden 1984)(3.3)(no sums), where . Note that some of these components are not independent because here λ1=λ2. In particular, we have(3.4)

### (b) Solutions

We look for solutions that are periodic along the circumferential and axial directions, and have yet unknown variations through the thickness of the tube, so that our ansatz is(3.5)where n=0, 1, 2, … is the circumferential number; k is the axial wavenumber; the subscripts (r, θ, z) for u and refer to components within the cylindrical coordinates (r, θ, z); and all upper-case functions are functions of r alone.

The specialization of the governing equations (3.1) to this type of solution has already been conducted in several articles (see Wilkes 1955; Fosdick & Shield 1963; Mack 1989; Pan & Beatty 1997; Negron-Marrero 1999; and Dorfmann & Haughton 2006 for the compressible counterpart). Here we adapt the work of Shuvalov (2003a) on waves in anisotropic cylinders to develop a Stroh-like formulation of the problem (Stroh 1962). The central idea is to introduce a displacement–traction vector,(3.6)so that the incremental equations can be written in the form(3.7)where G is a 6×6 matrix, with the block structure(3.8)Here the superscript ‘+’ denotes the Hermitian adjoint (transpose of the complex conjugate) and G1, G2 and G3 are the 3×3 matrices(3.9)respectively, with(3.10)

As it happens, there exists a set of six explicit Bessel-type solutions to these equations when n≠0. This situation is in marked contrast with the corresponding set-up in linear anisotropic elastodynamics, where explicit Bessel-type solutions exist only for transversely isotropic cylinders with a set of four linearly independent modes and do not exist for cylinders of lesser symmetry (Martin & Berger 2001; Shuvalov 2003b). As mentioned in §1, the six Bessel solutions are presented in the article by Wilkes (1955; for a derivation see Bigoni & Gei 2001).

First, denote by , the roots of the following quadratic in q2:(3.11)Then the roots of this quartic in q are ±q1 and ±q2, and it can be checked that the following two vectors are solutions to (3.7):(3.12)where q=q1, q2 in turn and In is the modified Bessel function of order n. Similarly, we checked that the following vector η(3):(3.13)is also a solution when q3 is the positive root of the quadratic equation(3.14)Finally, we also checked that the vectors η(4), η(5) and η(6), obtained by replacing In with the modified Bessel function Kn in the expressions above, are solutions too.

Next, we follow Shuvalov (2003a) and introduce (r) as a fundamental matrix solution to (3.7):(3.15)It clearly satisfies(3.16)Let M(r,a) be the matricant solution to (3.7), i.e. the matrix such that(3.17)It is obtained from (r) (or from any other fundamental matrix made of linearly independent combinations of the η(i)) by(3.18)and it has the following block structure:(3.19)say.

### (c) Boundary conditions

Some boundary conditions must be enforced on the top and bottom faces of the tubes. Considering that they remain plane (Uz=0 on z=0, l) and free of incremental shear tractions ( on z=0, l) leads to(3.20)where m=1, 2, 3, … but, since the equations depend only on k, we can take m=1 without loss of generality.

The other boundary conditions are that the inner and outer faces of the tube remain free of incremental tractions. We call the traction vector and the displacement vector. We substitute the condition S(a)=0 into (3.17) and (3.19) to find the following connection:(3.21)is the (Hermitian) 3×3 impedance (Shuvalov 2003a). Since S(b)=0, a non-trivial solution exists only if the matrix z(b, a) is singular, which implies the bifurcation condition(3.22)This is a real equation since z=z+ (Shuvalov 2003a) that applies independently of the nature (i.e. real or complex (Pan & Beatty 1997), simple or double (Dorfmann & Haughton 2006)) of the roots q1, q2 and q3.

We are now in a position to use the bifurcation condition (3.22) to compute explicitly bifurcation curves for each mode n. We note that the components of 0 depend on the strain energy density W and on the pre-strain, which by (2.3) depends only on λ3; so do q1, q2, q3, by (3.11) and (3.14). According to (3.12) and (3.13), the entries of M(b, a) thus depend (for a given W) on λ3, n, ka and kb only. For a given material (W specified) with a given thickness (b/a=B/A specified), the bifurcation equation (3.22) gives a relationship between a measure of the critical pre-stretch: , and a measure of the tube slenderness: , for a given bifurcation mode (n specified). That is, for a given tube slenderness, what is the axial strain necessary to excite a given mode?

While this bifurcation condition is formally clear, it has not been successfully implemented to compute all bifurcation curves. Indeed, for mode n>1, the root finding of det(z) becomes numerically unstable and numerical methods become unreliable (as observed in Dorfmann & Haughton (2006) for a similar problem) and, in explicit computations, most authors do not use the exact solution by Wilkes but use a variety of numerical techniques to solve the linear boundary-value problems directly (such as the compound matrix method (Haughton & Orr 1997), the determinantal method (Ben Amar & Goriely 2005) or the Adams–Moulton method (Zhu et al. 2008)). Note that from a computational perspective, the Stroh formalism is particularly well suited and well behaved (Biryukov 1985; Fu 2005) and if numerical integration was required it would provide an ideal representation of the governing equation.

Rather than integrating the original linear problem numerically, we now show how to use an alternative form of (3.22) to compute all possible bifurcation curves. This method bypasses the need for numerical integration and reduces the problem to a form that is manageable both numerically and symbolically, to study analytically particular asymptotic limits. The main idea is to transform condition (3.22) by factoring non-vanishing factors. We start by realizing that since the fundamental solutions {η(i), i=1, …, 6} are linearly independent, the matrix (r) is invertible for all r∈[a,b], which implies that the elements of M(r, a) are bounded for r∈[a,b]. Therefore, det(M1(r, a)) is uniformly bounded away from zero and det z=0 implies det(M3(b, a))=0. Instead of expressing det(M3(b, a)) as the determinant of a 3×3 submatrix of a matrix obtained as the product of two 6×6 matrices, we first decompose (r) as(4.1)say, where each block is a 3×3 matrix. We also rewrite equation (3.18) as(4.2)and write explicitly the two entries 3(r) and 4(r), which are(4.3)(4.4)which implies(4.5)Using again the fact that the entries of are bounded, we have that the bifurcation condition det(M3(b, a))=0 implies that(4.6)where(4.7)and adj(A) is the adjugate matrix of A, i.e. the transpose of the cofactor matrix (which in the case of an invertible matrix is simply adj(A)=det(A)A−1). This new bifurcation condition is equivalent to the previous one but has many advantages. The matrix Q involves only products of 3×3 matrices and is polynomial in the entries of , i.e. det Q(b, a) is a polynomial of degree 18 in Bessel functions and has no denominator (hence no small denominator). Both numerically and symbolically, this determinant is well behaved, even in the limits a→0, which corresponds to a solid cylinder, and n=0, which corresponds to the first barrelling mode (and usually requires a special treatment). We will refer to the use of this form of the bifurcation condition as the adjugate method.

### (a) Numerical results

As a first test of the stability of the numerical procedure, we consider a neo-Hookean potential , where we set C1=1 without loss of generality and consider the typical value B/A=2. We compute the critical value of λλ3 as a function of the current stubbiness kb=πb/l (the initial stubbiness is ) for the first 9 modes (n=0–8) as shown in figure 2. The known classical features of the stability problem for the cylindrical shell are recovered, namely for slender tubes, the Euler buckling (n=1) is dominant and becomes unavoidable as the slenderness increases; there is a critical slenderness value at which the first barrelling mode n=0 is the first unstable mode (in a thought experiment where the axial strain would be incrementally increased until the tube becomes unstable); and for very large kb, the critical compression ratio tends asymptotically to the value λ=0.444, which corresponds to surface instability of a compressed half-space (Biot 1962).

Figure 2

Bifurcation curves (stretch as a function of stubbiness) of a homogeneous neo-Hookean cylindrical tube for modes n=0–8 with b/a=B/A=2 and C1=1.

For a second test, we consider very thin neo-Hookean tubes with B/A=1.01. Here we are interested in the mode selection process. As the stubbiness increases, the buckling mode rapidly ceases to be the first excited mode and is replaced by different barrelling modes. From figure 3, it appears clearly that as kb increases, modes n=1–9 are selected (modes n=0 and 10 remain unobservable). There is one particularly interesting feature in these two sets of bifurcation curves. Depending on both the tube thickness and the stubbiness, the instability mode of a tube transition occurs from buckling to barrelling, the material transition from either the one-dimensional behaviour of slender column to the two-dimensional behaviour of a thin short tube, or the three-dimensional behaviour of a thick short tube. Accordingly, we will refer to these particular geometric values where transition occurs as dimensional transitions and obtain analytical estimates for them in §5.

Figure 3

Bifurcation curves (stretch as a function of stubbiness) of a homogeneous neo-Hookean cylindrical tube for modes n=0–10 with b/a=B/A=1.01 and C1=1.

## 5. Asymptotic Euler buckling

We are now in a position to look at the asymmetric buckling mode (n=1) corresponding to the Euler buckling in the limit λ→1. The asymptotic form of the Euler criterion cannot be obtained for a general strain-energy density. This is why we choose the Mooney–Rivlin potential, which, for λ close to 1, corresponds to the most general form of third-order incompressible elasticity (see §6 and Rivlin & Saunders (1951)),(5.1)where C1≥0 and C2>0 are material constants; ; and . Close to λ=1, we introduce a small parameter related to the stubbiness ratio(5.2)and look for the critical buckling stretch λ as a function of ϵ of order M,(5.3)Similarly, we expand in powers of ϵ,(5.4)and solve each order dm=0 for the coefficients λm. This is a rather cumbersome computation. It can be checked that λm vanishes identically for all odd values of m and that the first non-identically vanishing coefficient appears at order 24. A computation to order 28 is necessary to compute the correct expression for a, which is found to be to order 6 in ϵ(5.5)with(5.6)(5.7)(5.8)where ρB/A=b/a. It is of interest to compare the different approximations. We recover the Euler formula by keeping only the term up to ϵ2, which we denote by Euler2. We define similarly Euler4 and Euler6 by keeping terms up to orders 4 and 6 in ϵ. We show the different approximations as a function of ϵ2 for ρ=1.01 (figure 4a) and ρ=10 (figure 4b). The classical Euler formula is well recovered in the limit ϵ→0, but the Euler4 and Euler6 approximations clearly improve the classical formula for larger values of ϵ. It also appears from the analysis of Euler4 that for C2≥0 the classical Euler formula always underestimates the critical stretch for instability.

Figure 4

Comparison of the different Euler formulae obtained by expanding the exact solution to order 2 (the classical Euler buckling formula), Euler4 and Euler6 for a neo-Hookean potential C1=1, C2=0. For comparison purpose, we show the critical stretch for mode n=1 versus ϵ2 in which case the graph becomes linear in the limit ϵ→0. (a) ρ=b/a=B/A=1.01, (b) ρ=b/a=B/A=10.

## 6. Nonlinear Euler buckling for third-order elasticity

The analytical result presented in §5 was formulated in terms of parameters and quantities natural for the computation and the theory of nonlinear elasticity. In order to relate this result to the classical form of Euler buckling, we need to express equation (5.5) in terms of the initial geometric values A, B, L, the axial load acting on the cylinder and the elastic parameters entering in the theory of linear elasticity.

We first consider the geometric parameters. We wish to express the critical load as a function of the initial stubbiness ν=B/L and tube relative thickness ρ=B/A. Recalling that ϵ=πb/l and λ=l/L, b=λ−1/2B, we have(6.1)To express ϵ as a function of ν, we expand ϵ in powers of ν to order 6, and solve (6.1) to obtain(6.2)where λ(2) and λ(4) are defined in (5.6) and (5.7) and come from the expansion of λ in powers of ϵ.

Second, we want to relate the axial compression to the actual axial load N. To do so, we integrate the axial stress over the faces of the tubes, i.e.(6.3)Since σ3 is constant and given by (2.5), we have(6.4)

Third, we relate the elastic Mooney parameters C1 and C2 to the classical elastic parameters. Here, we follow Hamilton et al. (2004; Destrade & Saccomandi 2005) and write the strain-energy density to third order for an incompressible elastic material as(6.5)where μ is the usual shear modulus, or second Lamé parameter, and n3 is a third-order elasticity constant; μ is related to Young's modulus by E=3μ; also, in Murnaghan's notation, n3=n and in Landau's notation, n3=A (see Norris (1998) for other notations). In (6.5), i1, i2, i3 are the first three principal invariants of the Green–Lagrange strain tensor, related to the first three principal invariants I1, I2, I3 of the Cauchy–Green strain tensor by(6.6)Since I3=1, we can solve this linear system for i2 and i3 and write the strain-energy density (6.5) as a function of I1 and I2, i.e.(6.7)which by comparison with (5.1) leads to(6.8)

To write the nonlinear buckling formula, we consider (6.4) and first expand λ in ϵ using (5.5), then expand ϵ in ν using (6.2) and, finally, substitute the values of the moduli in terms of the elastic parameters, which yield(6.9)While it is not surprising, it is comforting to recover to order ν2 the classical Euler buckling formula (1.1) (using ρ=B/A, ν=B/L and μ=E/3).

## 7. Dimensional transition

Finally, we use the buckling formula to compute the transition between modes as parameters are varied. That is, to identify both the geometric values and the axial strain for which there is a transition between buckling and barrelling modes. Here we restrict again our attention to the neo-Hookean case (with C1=1). From figures 2 and 3, it appears clearly that for ϵ small enough there is a transition (depending on the value of ρ) from either mode n=1 to mode n=0 (large ρ), or from mode n=1 to mode n=2 (ρ close to 1) as ϵ increases. We refer to this transition as a dimensional transition, in the sense that the material mostly behaves as a slender one-dimensional structure when it buckles according to mode n=0 and mostly as a two-dimensional structure when it barrels with mode n=2. Indeed both modes of instability can be captured by, respectively, a one- or a two-dimensional theory. For ρ close to unity, the transition n=0→n=1 occurs for small values of ϵ. Therefore, in this regime, we can use the approximation (5.5) for the barrelling curve and substitute it in the bifurcation condition of mode n=2. Expanding again this bifurcation condition in ϵ as well as ρ, one identifies the values ρt of ρ and λt of λ at which the transition occurs, as(7.1)In terms of the initial stubbiness ν=B/L, we have(7.2)This relationship also provides a domain of validity for the Euler buckling formula. For sufficiently slender tube (ν small), the buckling mode disappears when ρ>ρt at the expense of the n=2 barrelling mode. For stubbier and fuller tubes, this approximation cannot be used. To understand the dimensional transition, we solve numerically the bifurcation condition, using the adjugate method, for the intersection of two different modes. That is, for a given value of ρ*, we find the value of ϵ* such that both the bifurcations for either modes n=1 and 2, or modes n=1 and 0 are satisfied. If the corresponding value λ* is the largest value for which a bifurcation takes place, the pair (ϵ*, ρ*) is a transition point. The corresponding transition point in terms of the initial parameters is . In figure 5, we show a diagram of all such pairs for both transitions.

Figure 5

Dimensional transition for a neo-Hookean cylindrical tube of initial length L and initial radii A and B. All tubes in the n=1 regions will become unstable by buckling. As the tubes get stubbier or thinner (arrows), it will not buckle but instead will be subjected to a barrelling instability. Note that only the transition curves from mode n=0 are shown. Tubes in the barrelling regions may be subjected to other unstable modes.

## 8. Conclusion

This article establishes a reliable and effective method to study the stability of tubes based on the exact solution of the incremental equations proposed by Wilkes (1955) within the Stroh formalism. It then puts the method to use, to obtain the first geometric and material corrections to the Euler buckling. The method can be also used to obtain the transition between buckling and barrelling modes when a tube becomes unstable.

The method presented here can be easily generalized to different materials and different boundary conditions. For instance, using the exact solution of the incremental equations proposed in Dorfmann & Haughton (2006) for compressible materials and the adjugate method, an explicit form of the bifurcation condition in terms of Bessel functions can be obtained by following the steps presented here and various asymptotic behaviours can be obtained. Similarly, a variety of boundary-value problems can be analysed by the adjugate method, such as the stability problem of a tube under pressure and tension (Han 2007; Zhu et al. 2008), the problem of a tube embedded in an infinite domain (Bigoni & Gei 2001) and the problem of a tube with coating (Ogden et al. 1997). In all these cases, useful asymptotic formulae for the buckling behaviour could be obtained by perturbation expansions.

It is also enticing to consider the possibility of performing an analytical post-buckling analysis of the solutions. Since the solutions of the linearized problem can be solved exactly, a weakly nonlinear analysis of the solution should be possible to third order. This would yield, in principle, an equation for the amplitude of the unstable modes containing much information not only about the actual amplitude of the unstable modes but also on the localization of unstable modes after bifurcation. We leave this daunting task for another day.

## Acknowledgments

This study is based in part upon work supported by the National Science Foundation under grant no. DMS-0604704 (A.G.) and made possible by a CNRS/USA Collaborative Grant from the French Centre National de la Recherche Scientifique (M.D.).