The effect of residual stress on the stability of a circular cylindrical tube

Residual stresses in an unloaded configuration of an elastic material have a significant influence on the response of the material from that configuration, but the effect of residual stress on the stability of the material, whether loaded or unloaded, has only been addressed to a limited extent. In this paper we consider the level of residual stress that can be supported in a thick-walled circular cylindrical tube of non-linearly elastic material without loss of stability when subjected to fixed axial stretch and either internal or external pressure. In particular, we consider the tube to have radial and circumferential residual stresses, with a simple form of elastic constitutive law that accommodates the residual stress, and incremental deformations restricted to the cross section of the tube. Results are described for a tube subject to a level of (internal or external) pressure characterized by the internal azimuthal stretch. Subject to restrictions imposed by the strong ellipticity condition, the emergence of bifurcated solutions is detailed for their dependence on the level of residual stress and mode number.

A theoretical basis for non-linearly elastic solids with residual stress was developed in a series of papers by Hoger, in particular in [3,4], while for a homogeneously initially stressed material (as distinct from a residually stressed material) further contributions detailing the effect of the initial stress on the propagation of waves were provided by Shams et al. [5], and Ogden and Singh [6], for example.
Although the basic analysis of the effect of residual stress on material response is now well developed, very little has been done on the effect that residual stress has on stability, an exception being the paper by Ciarletta et al. [7], which was concerned with the stability (in two dimensions) of an unloaded circular annulus carrying both radial and circumferential residual stresses.
In the present paper we examine the effect of axial load, internal and external pressure and residual stress on the stability of a circular cylindrical tube of incompressible elastic material based on the linear theory of elastic increments superimposed on the (finitely deformed) circular cylindrical configuration, with attention restricted to increments confined to the radial-circumferential plane.
We begin in Sect. 2 by providing the basic equations of elasticity with residual stress, while in Sect. 2.1 the constitutive equation is specialized in terms of invariants. The theory is then applied in Sect. 3 to the basic configuration of a residually stressed circular cylindrical tube subject to axial load and internal or external pressure. Section 4 summarizes the required incremental equations, which are specialized in detail in Sect. 5 for the case of planar increments (in the radial-circumferential plane). Specific forms of residual stress and constitutive equation to be used for illustration are also provided along with consideration of the restrictions imposed by the strong ellipticity condition. The numerical procedure used for solving the incremental equations and boundary conditions is described in Sect. 6.
Results of the calculations are illustrated in Sect. 6.3 to show the relative influences of the residual stress, the tube thickness, the deformed inner radius of the tube (as a measure of internal or external pressure) and the mode number on the onset of bifurcation, with particular reference to the strong ellipticity condition.
Some brief concluding remarks are provided in Sect. 6.4.

Residual stress and elasticity
We consider a solid continuum which, when unloaded, is designated its reference configuration, denoted B r . Within this configuration it possesses a residual stress distribution, with the residual (Cauchy) stress tensor denoted τ . By definition [3,4] τ satisfies the equilibrium equation and zero traction boundary conditions Div τ = 0 in B r , τ N = 0 on ∂B r , where ∂B r is the boundary of B r , which has unit outward normal N, Div being the divergence operator in B r , i.e. with respect to material position X ∈ B r . It is assumed that there are no intrinsic couple stresses so that τ is symmetric. Note that residual stresses are necessarily non-uniform, i.e. depend on X, and the material that supports them is inhomogeneous. When loads are applied, the resulting deformation is measured from B r with deformation gradient F, associated with which are the right and left Cauchy-Green deformation tensors, denoted C and B, respectively, and defined by The deformed configuration is denoted B and its boundary by ∂B, within which the material point X is located at x, which is given by the deformation function χ , such that x = χ (X) and F = Grad x, Grad being the gradient operator with respect to X.
In this paper we shall assume that the material response relative to B r is incompressible, so the constraint is satisfied for each X ∈ B r . We suppose that the response to applied loads is purely elastic and characterized in terms of a strain-energy function W defined per unit volume in B r and measured therefrom. It depends not only on the deformation gradient F but also the residual stress τ . Objectivity requires that W depends on F through C. The dependence on τ is included explicitly in the arguments of W and we write which is automatically objective since τ is unaffected by rotations in the deformed configuration B. Note that the elastic properties of the material relative to B r are anisotropic, i.e. τ has an effect on the constitutive law analogous to that of a structure tensor associated with preferred directions in fibre-reinforced materials, as for example, in [8].
In the special case that τ specializes to a rank-one tensor, say M ⊗ M, the invariants associated with transverse isotropy are recovered. In general, we may refer to τ as a generalized structure tensor. The appropriate invariants for a general τ will be listed in the following section.
Since W is measured from B r , we impose the condition The formulas for the nominal and Cauchy stresses are the same as in standard non-linear elasticity, except by the dependence of W on τ , which has a passive role since it is independent of the deformation. Thus, for the considered incompressible elastic material the nominal stress tensor S and Cauchy stress tensor σ are given by where p is a Lagrange multiplier necessitated by the constraint (3) and I is now written for the identity tensor in B (which we do not distinguish from that in B r ). For equilibrium in B the equations Div S = 0, div σ = 0 have to be satisfied in the absence of body forces. When F = I, each of the expressions in (6) reduces to where p (r) is the value of p in B r . The latter condition imposes restrictions on W in B r , which will be made explicit in the following.

Invariant formulation
Based on the general theory of Spencer [9] it can be seen that, for an incompressible material, in three dimensions W (C, τ ) depends on 9 invariants generated by the two tensors C and τ . With the invariant I 3 = det C omitted because of the constraint (3), these are typically taken to be, for C, for the combination of C and τ , while for τ we denote the three invariants collectively as The invariants of τ are not affected by the deformation, while in the configuration B r the other invariants reduce to We emphasize that the above set of 9 invariants, or an equivalent set of alternative invariants, forms a complete set of invariants of C and τ in three dimensions for an incompressible material. When the dimension of the considered problem is reduced from three to two, such as for plane strain, the number of independent invariants is reduced.
In terms of the invariants the expressions in (6) for the nominal and Cauchy stresses can be expanded in the form where I is the index set {1, 2, 5, 6, 7, 8}, the notationW i = ∂W /∂ I i , i ∈ I has been adopted, and W is written W (I 1 , I 2 , I 4 , I 5 , I 6 , I 7 , I 8 ) to reflect the dependence on the invariants. Note that the derivative of I 4 with respect to F vanishes and so is not included in the above expression, but I 4 is included in the functional dependence ofW . The expressions for ∂ I i /∂F, i ∈ I, required in (13), are given by the symmetry of τ having been used. The Cauchy stress in (13) then expands in full as wherein we have introduced the notation Σ = Fτ F T for the Eulerian tensor which is the push-forward of τ . We also recall that B = FF T is the left Cauchy-Green tensor. When evaluated in B r the latter reduces to the specialization of (8) as where eachW i , i ∈ I, is evaluated for the invariants given by (12). Thus, following [5] with a slightly different notation, we obtain the restrictions on the strain-energy function in B r .
Further details of the residual stress formulation are given in [5] and, for a material also containing a preferred direction, in [6], while an application to a prototype problem involving an inhomogeneous deformation, that of (plane strain) azimuthal shear of a circular cylindrical tube, for a residually stressed material, is provided in [1]. In the following section we consider the specialization of the above equations to a circular cylindrical configuration with the residual stress confined to the cross-sectional plane of the cylinder.

Application to a thick-walled tube with residual stress
We now consider a thick-walled circular cylindrical tube which is subject to a uniform axial stretch and radial deformation. Let the cross-sectional geometry of the tube in the configuration B r be defined by 0 < A ≤ R ≤ B, 0 ≤ Θ ≤ 2π , with axial coordinate Z . After deformation the internal and external radii A and B become a and b and since the material is incompressible the deformation is given in terms of cylindrical polar coordinates r, θ, z by where the constant λ z is the axial stretch, while b is given by The circumferential stretch, denoted λ, is given by λ = r/R, and by incompressibility the radial stretch is λ −1 λ −1 z . The residual stress in B r is taken to consist of radial and circumferential components τ 33 and τ 11 , respectively, with axial component τ 22 = 0. For the most part the components of tensors are signified with indices p, q, α, β, ..., etc., which take values 1, 2, 3 in the general case. The components τ 33 and τ 11 satisfy the equilibrium equation with the boundary condition (1) 2 specializing to Note that τ 11 and τ 33 are the only components of the residual stress since the considered geometry can support neither a shear nor an axial component of residual stress without Z dependence. The first two invariants in I 4 in (11) reduce to just τ 11 + τ 33 and τ 11 τ 33 and the third vanishes, while I 1 and I 2 in (9) become and I 5 , . . . , I 8 can be expressed in terms of τ 11 , τ 33 and the stretches as Since the invariants depend on only two independent stretches λ and λ z together with the residual stress components τ 11 and τ 33 , we introduce the reduced form of the strain-energy function, denotedW (λ, λ z , τ 11 , τ 33 ), and using (17) we then obtain the Cauchy stress differences the indices 1, 2, 3 corresponding to the coordinates θ, z, r , respectively. In (25) the subscripts λ and λ z signify partial derivatives.
For the considered geometry the equilibrium equation (7) 2 reduces to For the moment we assume that the radial deformation is accompanied by both internal and external pressures, denoted P a and P b on r = a and r = b, respectively. Thus, σ 33 satisfies the boundary conditions Integration of (26) then gives Note that positive (negative) P corresponds to an effective internal (external) pressure.
Since τ 22 = 0, when the expressions (25) are evaluated in B r , where λ = λ z = 1, they reduce to which provide restrictions onW that specialize those in (19). The corresponding expression for the reduced axial load on a cross section of the tube is The latter provides an expression for the reduced axial load required to maintain a fixed value of the axial stretch, while (28) provides an expression for the pressure difference for a prescribed value of the inner radius a and the axial stretch λ z , on recalling that the outer radius is given by

Incremental equations
In order to examine the linear stability of the basic configuration we need to consider incremental deformations relative to that configuration. Incremental quantities are signified by a superposed dot, so that the incremental displacement is denotedẋ =χ (X) and the corresponding increment in the deformation gradient byḞ, witḣ F = Gradẋ = LF, where L = grad u is the displacement gradient, u being the Eulerian counterpart ofẋ defined by u(x) =χ (χ −1 (x)). Then, on taking the increment of (6) 1 we obtain the increment in the nominal stresṡ whereṗ is the incremental Lagrange multiplier and the fourth-order tensor of elastic moduli A is defined in symbolic and index notation by The incremental equilibrium equation is DivṠ = 0.
On updating the reference configuration from B r to B the incremental constitutive equation (30) becomeṡ whereṠ 0 is the push-forward version ofṠ defined byṠ 0 = FṠ and satisfying the Eulerian form divṠ 0 = 0 of the incremental equilibrium equation, while A 0 is the corresponding push-forward of A defined in index notation by The incremental traction per unit area of the boundary ∂B isṠ T 0 n and will be made explicit later for the considered problem, n being the outward unit normal to ∂B.
Let e 1 , e 2 , e 3 be unit basis vectors in an orthogonal curvilinear coordinate system. Then, in component form, the equilibrium equation divṠ 0 = 0 yields the three scalar equationṡ as in Haughton and Ogden [10], in which summation over repeated indices j and k from 1 to 3 is implied and the subscript notation, j represents the derivative associated with the jth curvilinear coordinate. These components will be made explicit in the following section for the considered cylindrical polar coordinates. Equation (34) is valid for any elastic material whether or not it carries residual stress.
In respect of the invariant expansion, the components A 0 piq j have been given in [5] in full generality for three dimensions. These are very lengthy and not needed here, but we do note the general connection coming from incremental rotational balance equation (see [5]), which will be used here together with the incremental incompressibility condition in the form div u = 0. (36)

Planar bifurcation
We now order the coordinates as θ, z, r , associated with the stretches λ, λ z , λ r , respectively, and for the associated components of vectors and tensors we adopt the indices 1, 2, 3 in the same order, as already noted. The unit basis vectors associated with the cylindrical polar coordinates θ , z, r are denoted e 1 , e 2 , e 3 , and the derivatives (·) ,k in (34) denoted by subscripts with preceding commas become ∂/r ∂θ, ∂/∂z, ∂/∂r for k = 1, 2, 3, respectively. Here, however, we restrict attention to planar increments in the (r, θ) plane with all incremental quantities independent of z. Then, for the cylindrical polar coordinates the only non-zero scalar products e i · e j,k in (34) are In the absence of the incremental axial displacement the incremental displacement u may be written with components u and v corresponding to the radial and circumferential directions, respectively, and are in general dependent on both r and θ . The matrix of components of L = grad u with respect to the basis vectors e 1 , e 2 , e 3 is where the subscripts θ , r following a comma indicate the corresponding partial derivatives. The incremental incompressibility condition (36) specializes to which is satisfied by introducing a function φ(θ, r ) such that The incremental equilibrium equations are now used to obtain the equation for φ. First, we note that the equilibrium equation for i = 2 in (34) is satisfied identically, while the equations for i = 1, 3 yielḋ For the considered underlying cylindrical configuration the components ofṠ 0 in the above two equations are given bẏ where the components of the elastic modulus tensor A 0 are specialized below. On substitution of these expressions into (42) and (43) and use of the incompressibility condition (40) with (39) we obtaiṅ respectively, where a prime denotes differentiation with respect to r . It is now convenient to define the reduced notations Then, by eliminatingṗ r andṗ θ from (48) and (49) by cross differentiation and then substituting for the expressions (41) in the result, followed by some rearrangements, we obtain the required equation for φ, namely γ r 4 φ ,rrrr + 2βr 2 φ ,rrθθ + αφ ,θθθθ + 2(r γ + γ )r 3 φ ,rrr + 2(rβ − β)r φ ,r θθ The transition to this equation makes use of the equilibrium equation (26), the connections from (35), and We now turn to the boundary conditions associated with Eq. (51). For the considered pressure boundary conditions the incremental traction is given bẏ no incremental pressure being admitted. Hence, in terms of components this giveṡ By considering the boundary conditions involvingṠ 031 with the expression (46) and the connection from (52) 1 with the notation (50) 3 and using the boundary conditions (27), we obtain In general, apart from isolated configurations, we may assume that γ = 0. Then, in terms of φ, (56) becomes Next, on use of (47) and the incremental incompressibility condition in the form L 11 + L 33 = 0 the boundary conditions involvingṠ 033 can be rewritten as The term inṗ is then eliminated by differentiating (58) with respect to θ and using (48). Then, by using the expressions for the components L i j in terms of the function φ, the boundary conditions (27), and an application of (57) this yields the boundary condition We now have a partial differential equation (51) for φ as a function of r and θ and two boundary conditions in each of (57) and (59). To make further progress the system of equations and boundary conditions needs to be solved numerically. Towards this we consider solutions for φ of the form n being a non-negative integer. Fig. 1 Plots of the scaled circumferential and radial residual stress components, τ 11 /ν (continuous curve) and τ 33 /ν (dashed curve), respectively, for (51) is now expressed in terms of the function f n (r ) and its derivatives as The boundary conditions (57) and (59) become The foregoing equations and boundary conditions apply for any radial deformation with λ > 0 and axial stretch λ z > 0. However, to avoid buckling-type instabilities (z-dependent) associated with axial compression we restrict attention to λ z ≥ 1 since we are only considering prismatic bifurcations here.

Specialization of the residual stress and constitutive equation
At this point we specialize the form of the residual stress in order to provide explicit illustrations of the effect of residual stress on the material response. First, we choose a simple form of τ 33 satisfying the boundary conditions (22) and then use (21) to determine τ 11 . Specifically, we take and hence where ν is a constant, which may be positive or negative. Plots of τ 11 /ν and τ 33 /ν are shown in Fig. 1. For ν > 0, their behaviours are very similar to those arising from the residual stresses calculated for a single-layer artery wall from the so-called 'opening angle' method [11] or the assumption that the circumferential stress at a typical physiological pressure is uniform [12], and are therefore realistic in this context. However, it is appropriate to allow for the possibility that ν be negative in other contexts.
In order to obtain explicit results we now consider a simple form of energy function, namelȳ where μ > 0 is a constant and κ is a non-negative constant, as introduced in [5] in a slightly different notation, and also adopted in [1] and [2], for example. Thus,W = 0 and (19) 1 is satisfied with p (r) = μ in B r . From Sect. 3, for the considered deformation we haveW (λ, λ z , τ 11 , τ 33 ), which becomes and the stress differences (25) become Noting that the energy function (67) is inhomogeneous, the total energy per unit length of the tube, denoted W, is and this can be obtained explicitly on substitution from (67), use of (64), (65), (20) 1 and λ = r/R for fixed λ z and given material constants. The resulting expressions are quite lengthy and not therefore included here. For the case κ = 0, Fig. 2a shows example plots of the total energy W divided by 2πμ, denotedW, for four values of the dimensionless residual stress parameterν = ν A 2 /μ with λ z = 1 and B/A = 1.25. Thus, the total energy is positive, although for some R and λ the local energy (67) can be slightly negative. For κ > 0 the additional term is positive and the energy is greater than for the κ = 0 case. In dimensionless formκ = κμ. This is illustrated in Fig. 2b forκ = 0.5 for four values ofν. In this case the local energy (67) is positive for all λ except for a very small negative value near λ = 1. Example plots of the pressure from Eq. (28) for the energy function (67) with λ z = 1 using (64) and (65) are shown in Fig. 3.

The elasticity tensor and strong ellipticity
For the strain-energy function (66), we obtain, by specializing the general expressions given in [5],  (50) and (71) that These expressions are required in Eq. (61) and the boundary conditions (63). Additionally, in (61), recalling that a prime indicates differentiation with respect to r , we require the expressions Explicit forms of some of these expressions are fairly lengthy and not listed here but are easily obtained from (72)-(74). The strong ellipticity condition imposes restrictions on the material response (see, for example, [5]). In general the strong ellipticity condition can be written A 0 piq j n p n q m i m j > 0, where m i and n i , i = 1, 2, 3, are the components of non-zero vectors m, n satisfying the incompressibility condition m · n = 0. For m and n restricted to the 1, 3 (θ, r ) plane we take n = (n 1 , 0, n 3 ), m = (m 1 , 0, m 3 ) = (n 3 , 0, −n 1 ). This yields for all non-zero n 1 , n 3 . This requires α > 0 and γ > 0, and for the expressions given in (72) it follows that β > 0 (note that more generally the inequality β > − √ αγ would be required [13], but is automatically satisfied here).

Numerical solution
In order to solve Eq. (61) we form the system of first-order equations with dependent variables y = (y 1 , y 2 , y 3 , y 4 ), where y 1 = f n , y 2 = y 1 , y 3 = y 2 , y 4 = y 3 , the matrix M being given by The corresponding boundary conditions (62) and (63) become where the 2 × 4 matrix B is
Note thatν andκ were introduced earlier but are included here for completeness. The dimensionless forms of the required derivatives arē where the prime attached to a quantity with an overbar denotes a derivative with respect tor . Equation (79) is then written in the dimensionless form where the components ofȳ are given byȳ 1 = y 1 ,ȳ 2 = Ay 2 ,ȳ 3 = A 2 y 3 ,ȳ 4 = A 3 y 4 , and wherē B = n 2 −rr 2 0 (2β +γ )n 2 −(2β +γ )n 2r 0γr 3 . (94) The coefficients in (92) have the forms 6.2 Restrictions from strong ellipticity From Sect. 5.2 and (72), the requirement α > 0 in dimensionless form yields and γ > 0 entails In each case the inequalities should be satisfied for allR ∈ [1,B] and govern the allowed ranges of values of ν, λ, λ z . For κ = 0 these are independent of the deformation and reduce to which restricts significantly the allowable values ofν. ForB = 1.2, for example, this gives −25/6 <ν < 5, and the restrictions become tighter asB increases, as Fig. 4 shows. Note that the upper condition on the right-hand side of (101) was omitted from equation (62) in [2]. When κ = 0 the restrictions are different and depend on the deformation. Figure 5 provides example plots ofν versusā (the circumferential stretch on the inner boundary) showing the (shaded) regions where the strong ellipticity condition is satisfied (i.e. α > 0, γ > 0). In particular, the representative plots are for λ z = 1 andB = 1.2 for several values ofκ. Each value ofR ∈ [1, 1.2] yields a different region, and in each case the region shown is that for which strong ellipticity holds for allR ∈ [1, 1.2]. The character of the results is similar for other values ofB. Thus, in some cases, the quadratic term in the energy function increases the ranges of values ofν andā for which strong ellipticity holds compared with the case κ = 0, while in other cases the ranges of allowable values ofν ofā are reduced. Asκ increases the area of the region where strong ellipticity holds becomes smaller and smaller.
The analogues of the caseB = 1.2 in Fig. 5b forB = 1.1, 1.5 are shown in Fig. 6. As is the case with κ = 0 (Fig. 4), the region where strong ellipticity holds becomes smaller asB increases.

Results
The aim now is to determine the numerical solution , outside these regions either α < 0 or γ < 0, or both ofB at the point of onset of bifurcation is determined for different mode numbers, applied (internal or external) pressure (as measured by the dimensionless internal radiusā) and fixed axial stretch (exemplified by λ z = 1). The system is solved for the dimensionless forms (90) with (93) using the routine NDSolve in Mathematica [14].
The results are illustrated for the model (66), first of all with κ = 0, in Fig. 7, for λ z = 1, where the curves of ν versusā for which bifurcation is possible are plotted for a series of values of the mode number n and for each ofB = 1.1, 1.2, 1.5 separately. Results for other values of λ z (> 1) are very similar and not shown separately. For example, with λ z = 1.5 the curves are essentially just shifted to the left with minor numerical differences. With reference to the inequalities (101) the limits imposed by strong ellipticity are shown as horizontal red lines in Fig. 7. The valid region in (ā,ν) space is between the lines in each case.
As can be seen from Fig. 7 most of the bifurcation curves in the strongly elliptic region are whereā < 1, i.e. when there is no internal pressure, which is consistent with the situation where there is no residual stress discussed in [10]. However, it can be seen from Fig. 7a, b that there is a small range of negative values ofν where bifurcation can occur forā > 1, at least for mode number n = 2. This is not the case in Fig. 7c. In the domainā < 1, corresponding to external pressure, and subject to strong ellipticity, asā is reduced from 1 the n = 2 mode occurs first, followed by n = 3, 4, ... in sequence until there is some interchange of priorities for higher modes. This is of minor theoretical interest since the mode n = 2 dominates. Note, in particular, with reference to Fig. 7a, b, that within the strongly elliptic region bifurcation may arise for an unloaded cylinder, i.e. whenā = 1, for limited mode numbers. In the regionā > 1 and outside the strongly elliptic domain, bifurcation is possible for higher-order mode numbers, depending on the value of the residual stress parameterν. To avoid overcrowding, each panel of Fig. 7 shows a representative selection of curves for a limited set of mode numbers n. Figure 8 illustrates the corresponding results for the model (66) for a non-zero value ofκ, exemplified byκ = 0.5, again forB = 1.1, 1.2, 1.5. There are significant similarities with Fig. 7 and this is also the case for larger values of κ (not shown). One difference is that the strong ellipticity domain depends onā and spreads into a larger area wherē a > 1.2, as shown in more detail in Figs. 5b and 6 (forκ = 0.5), the pressure then being internal. Whenā = 1 (the unloaded case), the situation is different from that for κ = 0 as no bifurcation curves arise. Higher modes can occur within the strongly elliptic region forā greater than approximately 1.2, but these are not shown. The area of the strongly elliptic domain decreases asB increases, as with the case κ = 0.

Concluding remarks
The influence of residual stress on the incremental response and bifurcation of a circular cylindrical tube subject to axial extension and internal or external pressure has been investigated. Incremental deformations superimposed on the cylindrical configuration have been restricted to the planar cross section of the tube, i.e. independent of the axial coordinate z. Results have been presented for an unstretched tube (λ z = 1) only as results for λ z > 1 are qualitatively similar, while those for λ z < 1 require consideration of possible buckling modes (axisymmetric or asymmetric) that depend on z. The results are described for a prototype model of elasticity that incorporates residual stress, and they quantify the dependence of the dimensionless residual stress parameterν on the dimensionless internal radius a of the tube, which provides a measure of the internal or external pressure, for different tube thickness ratiosB and selected bifurcation mode numbers n. Restrictions imposed by the strong ellipticity condition form a framework for the interpretation of the results. In particular, limits on the range of the allowable residual stress magnitude are determined. It has been found that, with limited exceptions, it is the case that bifurcation can only occur when the pressure is external, and lower mode numbers have priority, as is the case for the situation without the presence of residual stress [10], where it was shown that bifurcation under internal pressure was excluded. In particular, mode n = 2 occurs first as the deflation of the tube proceeds from the point where the internal radius is at its initial value (ā = 1). For the considered model, higher modes, i.e. for n greater that about 48, can in principle arise within the strongly elliptic domain, but these are of minor interest and can occur only for small values ofν. It is worth noting in passing that for very large values of n → ∞ the only solution of the governing equations is the trivial solution f n (r ) → 0. For the special case in which the tube is unloaded (λ z = 1 andā = 1) bifurcation (wrinkling) solutions were obtained for a different model from that adopted here and illustrated in the paper by Ciarletta et al. [7]. We have also obtained solutions for their model and found results which are very similar to those in Figs. 7 and 8.
In a subsequent paper attention will be devoted to axisymmetric and asymmetric modes of bifurcation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.