Wrinkling of Elastic Cylinders With Material Properties Varying in Radial Direction

Although the instability of graded elastic cylinders has been analyzed by many researchers, most of them focused on the core-shell cylinders and film-substrate structures with inhomogeneous Young’s modulus. For a radially graded elastic cylinder subjected to the axial compression, the variation of Poisson’s ratio may result in the radial and circumferential stresses and thereby affects the critical condition of instability. By assuming linear elasticity with nonlinear kinematics, the governing equation for the incremental stress field is developed for instability analysis of the cylinder with radially graded material properties (Young’s modulus and Poisson’s ratio). Considering the arbitrariness of material properties, the state space technique is implemented and a semi-analytical solution is acquired. The obtained solution is validated by the finite element results. Numerical examples show that the critical condition of instability for graded elastic cylinders is related to whether Poisson’s ratio is assumed to be constant.


INTRODUCTION
The morphology characteristics caused by surface instability are often observed in nature and biological tissues, such as fruits, vegetables, animal brains, human guts (Dai and Liu, 2014;Budday et al., 2014;Balbi et al., 2015), and have been exploited for a broad range of applications, e.g., testing of mechanical properties for ultrathin films , active surfaces (Tokarev and Minko, 2009), micro-optics , microfluidic devices (Sugiura et al., 2007), and soft electronics and actuators (Rogers, et al., 2010;Yang et al., 2010). Due to curiosity about the natural phenomena and motivation of the development and utilization based on these configurations, the intrinsic mechanism of surface instability has attracted much attention of a large number of researchers over the last decades.
For the surface wrinkling instability, there are many previous works on core-shell cylinders and film-substrate structures. The pioneering work on film-substrate structures originated from the delamination analysis of such systems with a large unbounded flaw under compressive stress (Evans and Hutchinson, 1984;Hutchinson and Suo, 1992). The instability of an elastic film on a viscous substrate was then studied profoundly (Sridhar, et al., 2001;Huang and Suo, 2002a;Huang and Suo, 2002b). The relevant analyses of stiff films on compliant substrates were later presented ceaselessly in a large numbers. Huang et al. (2005) provided a nonlinear analysis of a film bonded to a substrate with finite thickness. Using the finite deformation theory, Song et al. (2008) studied the mechanics of thin buckled films on compliant substrates and obtained an analytical solution. Jia et al. (2012) investigated the wrinkling problem of a bilayer film on a soft substrate via theoretical analysis and finite element simulations. Xu et al. (2015) studied the occurrence and evolution of 3D instability patterns in thin films on hyperelastic substrates based on a fully nonlinear 3D finite element model. In comparison with the filmsubstrate structure, the stability analysis of cylindrical shells with soft cores was considered earlier by Seide (1962) using Donnell's equations. About the same time a theoretical solution for the determination of buckling characteristics was presented for a long cylinder shell with a solid or elastic soft core under axial compression with the sinusoidal buckling assumption (Yao, 1962). Inspired by optimization of natural structures, Karam and Gibson (1995) theoretically analyzed the buckling of coreshell structures and found that the buckling resistance of a hollow cylindrical shell could be improved significantly by infilling a compliant elastic core. Hutchinson and He (2000) addressed the buckling of cylindrical sandwich shells with foamed metal cores for the optimal design for structural weight and carrying capacity, and a similar comprehensive analysis was provided by Dawson and Gibson (2007). By energy method, Arani et al. (2007) studied the elastic axisymmetric buckling of a thin simply supported cylindrical shell with an elastic core under axial compression. To further investigate the effect of the filled core rigidity and thickness on the buckling behavior, Ye et al. (2011) developed a simple formula to predict the critical buckling stress by using the Rayleigh-Ritz approximation. Based on the nonlinear Donnell-Mushtari-Vlassov shell theory, Zhao et al. (2014) investigated surface wrinkling of a core-shell cylinder along the axial and circumferential directions. Recently, the axisymmetric and non-axisymmetric mode transitions in the post-buckling regime in axially compressed core-shell cylinders were investigated based on a nonlinear 3D finite element model (Xu and Potier-Ferry, 2016) and nonlinear shell finite element (Lavrencic et al., 2020).
The aforementioned literature only involves the film-substrate and core-shell systems with uniform substrates and cores. In nature and some engineering problems, the material properties of substrates and cores are commonly inhomogeneous, and thus the surface wrinkling of such systems is necessary to be studied. However, the inhomogeneity of materials makes it difficult to obtain a completely analytical solution. In the past few years, some theoretical researches on surface instability have been devoted to graded elastic materials. Cao et al. (2012a) derived the critical condition of wrinkling for a stiff thin layer on a semiinfinite substrate with the power-law and exponential grading moduli, respectively. For elastic layers with material properties varying from the surface to the interior, a semi-analytical solution was obtained by the theoretical analysis combined with the finite element method (Lee et al., 2008) and the state space method (Wu et al., 2014); similar semi-analytical solutions for onset of surface instability of elastic cylinders with radially graded Young's modulus were performed by Jia et al. (2014) and Han et al. (2017). Based on Floquet's principle, most recently an explicit semi-analytical isogeometric analysis method was suggested for predicting wrinkling instability of a stiff film on a graded compliant cylinder (Li et al., 2020). We note in the existing studies that Poisson's ratio in graded elastic layers and cylinders was commonly assumed to be a constant. In practice, however, for a radially graded elastic cylinder subjected to the axial compression, the variation of Poisson's ratio may result in radial and circumferential stresses, which must affect the critical condition of instability.
In addition to be directly subjected to compression, surface instability of film-substrate or core-shell systems may be induced by swelling or growing, and surface stability of soft elastic layers may be also influenced by surface tension (Kang and Huang, 2010a;Wang, 2020). Upon swelling in a solvent, the surface of thin hydrogel layers on rigid substrates may become unstable and further evolves into various morphologies (Guvendiren et al., 2009;Guvendiren et al., 2010;Kang and Huang, 2010b;Wu et al., 2017). The tubular organ of animals, such as esophagus, can be seen as cylindrical bilayers with the stiffer outer layer (muscle) and softer inner layer (mucosa). When the inner mucosa grows, the outer muscle restricts its expansion, eventually leading to wrinkling, folding or creasing on the inner surface of mucosa (Jin et al., 2011;Li et al., 2011;Razavi et al., 2016). Similar instability was also found on the outer surface of core-shell cylinders as the outer shell grows or swells on a stiff core Cao et al. (2012b). The surface instability related to growth and swelling has been adequately studied in the past 10 years.
In this paper, we focus on an inhomogeneous elastic cylinder and assume that both Young's modulus and Poisson's ratio vary arbitrarily in the radial direction. A state space method is employed to investigate the onset of wrinkling for elastic cylinders under axial compression. The effects of the variation of Poisson's ratio on stress distribution and critical condition of instability are discussed. The rest of this paper is arranged as follows. In Section 2, the instability theory of elastic cylinders is introduced. The solution of stress field in fundamental state is obtained by using the state space method in Section 3. In Section 4, the instability analysis for elastic cylinders with material properties varying arbitrarily in the radial direction is carried out and the critical condition of instability is presented.
Two instability examples, a cylinder covered by a bilayer and a homogeneous cylinder with a linearly graded stiff layer, are discussed in Section 5. Concluding remarks are presented in Section 6.

INSTABILITY THEORY OF ELASTIC CYLINDERS
An elastic hollow inhomogeneous cylinder with traction free surfaces is considered as shown in Figure 1A, where a Cartesian coordinate system is set up with the reference coordinates X 1 and X 2 at one cross-section and X 3 along the axis of the cylinder. The inner and outer radii of the cylinder in the stress-free state are A and B, respectively. When subjected to an axial compression, the compressive stress inside the cylinder may cause surface or internal instability of the cylinder. Prior to surface or internal instability, the compressed cylinder is seen as in the fundamental state with the inner and outer radii a and b ( Figure 1B), and the corresponding current coordinates are denoted as (x 1 , x 2 , x 3 ).
In the present study, we consider an elastic cylinder with material properties varying in the radial direction, especially emphasizing the influence of Poisson's ratio on stress distribution and structural stability. The material is assumed to be linear elastic and isotropic with a quadratic strain energy density function in terms of Green-Lagrange strain: where E IJ = (F kI F kJ − δ IJ )/2, F kJ = zx k /zX J , δ IJ is the Kronecker delta, and the elastic modulus C IJKL is a function of X 1 and X 2 and possesses the isotropic symmetry. The first and second Piola-Kirchhoff stresses P iJ and S IJ are accordingly in the Cartesian coordinate system as (Wu et al., 2014) In the fundamental state, the equilibrium equation is and the boundary conditions at the inner and outer surfaces are P iJ N J 0, at X 2 1 + X 2 2 A and B (4) where the notation ( ) ,J stands for differentiation with respect to X J in the reference state, N J represents the direction cosine of the outer normal relative to the coordinate X J ; the Einstein summation convention is implied over repeated indices unless noted otherwise. Under the axial compression, the imposed axial nominal strain ε 0 = F 33 -1 is identical everywhere in the fundamental state. In the process of instability analysis, it is assumed that an incremental displacement Δu i (i = 1, 2, 3) occurs from the fundamental state, then the corresponding incremental stress field must satisfy the following equilibrium equation and boundary conditions: Assuming that the strain in the fundamental state is small, we have F iK ≈ δ iK and the equilibrium Eq. 5 and boundary condition (6) may be translated as that ΔP ij,j C ijkl Δu k,l + P kj Δu i,k ,j 0 ( 7 ) where ( ) ,j denotes differentiation with respect to x j in the current state, and n j represents the direction cosine of the outer normal relative to the coordinate x j . Considering axisymmetry of the problem, the cylindrical coordinates in the reference state (R, Θ, Z) and in the current state (r, θ, z) are applied as shown in Figure 1. With the coordinate transformation from the Cartesian coordinates to the cylindrical coordinates, the equilibrium Eq. 3 in the fundamental state for radially graded elastic cylinders can be converted into the following form: with u R representing the radial displacement component. The axial stress in the fundamental state is then determined by The boundary condition (4) becomes In general, Poisson's ratio changes with the change of Young's modulus of graded elastic materials. Poisson's ratio commonly changes very little and thus it is often approximated as a constant for convenience. When the radially graded elastic cylinder is subjected to axial compression, the axial compressive stress is generated definitely. However, the radial and circumferential stresses do not necessarily occur, which depend on Poisson's ratio. If Poisson's ratio remains constant in the radial direction, transversely homogeneous strain generates and no radial and circumferential stresses occur. In contrast, if Poisson's ratio varies in the radial direction, the radial and circumferential strains are inhomogeneous and corresponding stresses must be produced, which will inevitably affect the axial and circumferential stabilities of cylindrical structures. The instability analysis for Poisson's ratio fixed as a constant has been carried out by many researchers (Jia et al., 2014;Han et al., 2017). In the present study, we assume that Poisson's ratio varies in the radial direction as that of Young's modulus. By the coordinate transformation, the equilibrium equation of the incremental stress field in Eq. 7 can be converted into the following form in the cylindrical coordinate system: The boundary condition of the incremental stress field in Eq. 8 becomes ΔP rr ΔP θr ΔP zr 0, at r a and b (26) Equations 9-12 and Eqs 14-25 constitute the governing equations of instability analysis for radially graded elastic cylinders, which together with the boundary conditions (13) and (26) lead to an eigenvalue problem. If there is a nontrivial solution for the incremental displacement field, the critical condition for onset of instability is obtained successfully.
It should be noted that if the cylindrical structure is a solid cylinder, the boundary condition (13) at R = A needs to be replaced with and the boundary condition (26) at r = a replaced with ΔP θr ΔP zr 0 and Δu r 0, at r 0 (28)

SOLUTION OF STRESS FIELD IN FUNDAMENTAL STATE
Considering arbitrariness of the material properties varying in the radial direction, it is impossible to obtain the analytical solution of stress field in the cylinder, and the state space technique is thereby employed in this paper. From Eq. 10, we have Substituting Eq. 29 into Eq. 11 and then into Eq. 9, we get zP rR zR and the circumferential and axial nominal stresses are, respectively, where λ and μ are Lame's constants of isotropic elasticity, relating to Young's modulus E and Poisson's ratio ] as follows: Combining Eqs 29, 30, a set of inhomogeneous differential equations are obtained in a matrix form as where the coefficient matrix D and the inhomogeneous item K are, respectively, Since the coefficient matrix and inhomogeneous item in Eq. 35 contain the radial coordinate R and the radially varying material parameters, the analytical solution to Eq. 34 cannot be obtained directly. For this reason, we first divide the hollow cylinder into n subshells and number them from the inner subshell to the outer subshell. The thickness of the ith subshell is denoted by H i (i 1, 2, ···, n) as illustrated in Figure 1D, and its inner and outer radii are R i−1 and R i , respectively. Thus there are R i A + H 1 + H 2 + / + H i , R 0 A and R n B. The material parameters are accordingly dispersed into each subshell as piecewise constant functions E i , ] i . In each subshell, the material parameters are assumed as their corresponding values at its mid-surface, i.e., When n → ∞ and H i → 0, the virtual laminated cylinder with material properties varying as piecewise constant functions infinitely approaches the originally graded cylinder. Taking R (R i−1 + R i )/2 in Eq. 35, a set of the first-order inhomogeneous differential equations with constant coefficients for the ith subshell, i.e. the state equation, are obtained as z zR The relation between the state vectors at outer surface and the inner surface of the ith subshell is determined as Thanks to the continuity of the state vector across all interfaces between the adjacent subshells, the relationship between the state vectors at the outer and inner surfaces of the cylinder can be acquired as By using the boundary condition (13), i.e.
u R1 (0) 0 and P rRn (R n ) 0 the state variables at the inner and outer surfaces are subsequently obtained as functions of the axial strain ε 0 , and u Ri and P rRi (i 1, 2, /, n) at any interface may be calculated successively by Eq. 37. The other two stress components P θΘi and P zZi are in turn computed from Eqs 31, 32. As a result, all displacement and stress components in fundamental state are determined as functions of ε 0 .

INSTABILITY ANALYSIS
To acquire the critical condition of instability for cylinders with material properties varying arbitrarily in the radial direction, the virtual dividing model in the prior section is still demanded in the following solving process. From Eqs 14-25, we can get the differentiations of incremental displacements and stresses with respect to the radial coordinate r as follows: zΔP rr zr λ + 2μ + P θθ − λ 2 λ + 2μ + P rr Δu r r 2 − μ + P θθ − μ 2 μ + P rr z 2 Δu r r 2 zθ 2 + λ − λ 2 λ + 2μ + P rr zΔu z rzz With the assumption of small strain in the fundamental state, we have P rr ≈ P rR , P θθ ≈ P θΘ , and P zz ≈ P zZ , they are functions of the axial strain ε 0 and determined by the state space method in Section 3. According to Eqs 41-46, the incremental displacements and stresses may be assumed to predict the axial and circumferential instabilities as Δu r U r (r) cos ω 1 θ cos ω 2 , zΔu θ U θ (r) sin ω 1 θ cos ω 2 z, Δu z U z (r) cos ω 1 θ sin ω 2 z (47) ΔP rr T r (r) cos ω 1 θ cos ω 2 z, ΔP θr T θ (r) sin ω 1 θ cos ω 2 z, ΔP zr T z (r) cos ω 1 θ sin ω 2 z where ω 1 and ω 2 refer to the wave number along the circumferential and axial directions, respectively. Substituting Eqs 47, 48 into Eqs 41-46, we obtain a set of homogeneous differential equations in the matrix form as dG(r) dr MG(r) where the state vector G(r) = [U r (r) U θ (r) U z (r) T r (r) T θ (r) T z (r)] T , and M is a 6×6 coefficient matrix with the following nonzero elements: Taking E E i , ] ] i , and r (r i + r i−1 )/2 in the coefficient matrix M, the state equation for the i th subshell can be obtained as Integrating two sides of Eq. 50, the state vector of the outer surface relating to that of the inner surface is derived as By means of the continuity of the state vector across all interfaces, the relationship between the state vectors for the outer and inner surfaces of the whole cylinder is determined as where N 1 i n N i (h i ). Corresponding to Eq. 52, the boundary condition Eqs 26, 28 can be rewritten as that Since the coefficient matrix M i includes the wave numbers ω 1 and ω 2 , and the stress components P rRi , P θΘi , and P zZi (i 1, 2, /, n) are functions of the axial strain ε 0 , the relationship between the state vector G n (r n ) and G 1 (r 0 ) by Eq. 52 therefore depends on ω 1 , ω 2 and ε 0 . Using the boundary conditions Eqs 53, 54, the critical strain for onset of instability of radially graded elastic cylinders can be found and eventually written in the following form:

RESULTS AND DISCUSSION
For a validation of the present study, we first analyze the stress distribution for a graded elastic cylinder under axial compression, and then discuss the critical condition of instability for a soft homogeneous cylinder covered by a bilayer or including a linearly graded stiff layer, respectively. The numerical simulation is computed by the finite element analysis code ABAQUS. The stress field and the critical condition of instability for graded elastic cylinders with a varied Poisson's ratio are addressed and compared to that with a fixed Poisson's ratio.

Validation of Stress Field in Fundamental State
For a radially graded elastic cylinder subjected to axial compression, from the theoretical analysis we know that the radial and circumferential stresses will be induced due to the variation of Poisson's ratio. Consider a solid cylinder with the outer radius B = 100 mm, both Young's modulus and Poisson's ratio are assumed to vary linearly in the radial direction, i.e.
where y stands for Young's modulus E or Poisson's ratio ], y 0 and y B represent E or ] at R = 0 and B, respectively. Figure 2 depicts the radial displacement u R , radial stress P rR , circumferential stress P θΘ , and axial stress P zZ as functions of the radial coordinate R for a solid cylinder with E 0 = 1 MPa and E B = 100 MPa, comparing with the numerical results by ABAQUS. From Figure 2, it is observed that the present results are in excellent agreement with the results by ABAQUS. Figures 2A,D reveal that both the varied Poisson's ratio (] 0 = 0.45 and ] B = 0.35) and the fixed one (] 0 = 0.4 and ] B = 0.4) have little effectiveness on the radial displacement and the axial stress. However, the variation of Poisson's ratios has significant influence on the radial and circumferential stresses (Figures 2B,C). When Poisson's ratio is varied and ranges from 0.45 to 0.35, the radial and circumferential compressive stresses increase firstly and then decrease in the radial direction. The radial stress is consistently compressive within the cylinder and equal to zero at the surface, satisfying the traction free boundary condition, whereas the circumferential stress turns to be tensile as r > 60 mm, and the maximum tensile stress reaches 0.035 MPa. Comparatively, when Poisson's ratio is fixed as 0.4, the radial and circumferential stresses are equal or very close to zero, consisting with the theoretical analysis results. During the following instability analyses, Poisson's ratio is assumed to be varied, and thus the radial and circumferential stresses will affect the critical condition of instability.

Surface Instability for a Cylinder Covered by a Bilayer
As an example for surface instability analysis for a graded elastic cylinder, we consider a soft cylinder covered by a bilayer as shown in Figure 3. Similar to the example in Jia et al. (2014) and Han et al. (2017), the radius of the inner soft cylinder is h s = 90 mm, the thicknesses of the surface shell and the intermediate shell are h f = 1 mm and h i = 10 mm, respectively. Young's moduli for the surface shell and the inner cylinder are E f = 1 GPa and E s = 1 MPa. However, in this paper Poisson's ratios for three different layers are not fixed as 0.4. We here assume that Poisson's ratios for the surface shell and the inner cylinder to be ] f = 0.4 and ] s = 0.48, and that for the intermediate shell varies with its Young's modulus and is given by where E i is Young's module for the intermediate shell.
The critical strain for onset of surface instability ε c and the corresponding wavelength λ c 2π/ω c , normalized by h f , are plotted in Figure 4 for such a cylinder with the modulus ratio of the intermediate shell to the soft cylinder E i /E s varying from 1 to 1,000, and the varied Poisson's ratio ] i from 0.48 to 0.4 according to Eq. 57. In comparison with the results by Han et al. (2017) in which Poisson's ratio is fixed as ν i = ] f = ] s = 0.4, it is found that the critical strain for the varied Poisson's ratio is very close to that for the fixed one, and their corresponding critical wavelengths are also close to each other in most part of the interval of modulus ratio. Both of the critical wavelengths for the two cases have a sudden rise. The sudden rise point for the varied Poisson's ratio happens at E i /E s = 11.502, the dimensionless wavelength increases from 24.7 to 37.1, different from E i /E s = 17.025 for the fixed Poisson's ratio with the dimensionless wavelength from 18.0 to 57.5. Obviously, the increase amplitude of the wavelength for the varied Poisson's ratio is greatly reduced in comparison with that for the fixed one. The results show that the radial and circumferential stresses induced by the variation of Poisson's ratio inside the cylinder have little effect on the critical strain, but have some effect on the critical wavelength, especially on the point of critical wavelength switching for such a system.
It should be noted that during searching for the critical stain, there are two local minimums for the candidate of the critical strain ε c , corresponding to a short wavelength and a long wavelength (see Han et al., 2017). The smaller local minimum strain and the corresponding wavelength are the true critical strain and critical wavelength. The dotted lines in Figures 4A,B illustrate the bigger local minimum strain and the corresponding wavelength for the system with varied Poisson's ratio at a possible metastable state of surface instability.
To verify the correctness of the theoretical solution, we select two computation models with two different lengths for numerical simulation: one has the length 162.8 mm with a varied Poisson's ratio, equivalent to four times of the wavelength at E i /E s = 12.0; the other has 175.28 mm with a fixed Poisson's ratio, equivalent to three times of the wavelength at E i /E s = 17.5. Figures 5A,B present the surface wrinkling patterns for the cylinder with a varied Poisson's ratio, and their modulus ratios are E i /E s = 11.0 and 12.0, respectively. Figures 5C,D provide the surface wrinkling patterns for the cylinder with a fixed Poisson's ratio and their E i /E s = 16.5 and 17.5. It can be found that the wave numbers are 6.5 in Figures 5A, 4 in Figure 5B, respectively, indicating that the modulus ratio for critical wavelength switching must be within 11.0 and 12.0 for the varied Poisson's ratio. In Figures 5C,D, the wave numbers are 8.5 and 3, respectively, illustrating that the modulus ratio for critical wavelength switching is within 16.5 and 17.5 for the fixed Poisson's ratio. The numerical simulation manifests that the modulus ratios 11.502 and 17.025, predicted for critical wavelength switching by the theoretical analysis, are correct and enough accurate.

Instability of a Homogeneous Cylinder With a Linearly Graded Stiff Layer
Next we consider a homogeneous soft hollow cylinder containing a linearly graded stiff layer with the thickness h g = 1 mm as shown in Figure 6. The inner and outer radii of the cylinder are A = 40 mm and B = 100 mm, Young's modulus and Poisson's ratio are  In the numerical analysis, the stiff layer is considered to be the inner layer ( Figure 6A), the middle layer ( Figure 6B), and the outer layer ( Figure 6C) Figure 7 presents the critical strain and the corresponding wavelength versus the modulus ratio E m /E h ranging from 10 2 to 10 5 for soft cylinders with an inner, a middle, and an outer stiff layers. ] m = 0.4 and ] m = ] h = 0.48 are used for the stiff layer with the varied Poisson's ratio and fixed Poisson's ratio, respectively. It is seen that the critical strains for the stiff cylinder with varied and fixed Poisson's ratios are very close to each other in each case and even indistinguishable, while the critical wavelength for varied Poisson's ratio is obviously less than that for the fixed Poisson's ratio, showing that the change of Poisson's ratio mainly affects the critical wavelength. In addition, from the results for the three different cases we can find that the critical strain of the cylinder with an outer stiff layer is always less than that with an inner or a middle stiff layer, and the corresponding wavelength is greater than that of the other two cases overall. For the cylinders with an inner stiff layer and a middle one, there is an intersection point for both the critical strain and the corresponding wavelength.
Finally we provide three numerical simulations for the homogeneous hollow cylinders with an inner, a middle, and an outer stiff layers, respectively, as shown in Figure 8. As predicted in Figure 7B, there is an intersection point near the modulus ratio E m /E h = 10 3 for the critical wavelengths of the first and second cases. The exact dimensionless wavelengths at E m /E h = 10 3 for the three cases are 18.41, 18.68, and 24.91, respectively. Thus, there is an approximate common multiple for the three wavelengths, such as 74.73. In the numerical simulation, the lengths of cylinders are selected as 74.73 mm, equivalent to about three times the critical wavelength for the cylinder with an outer stiff layer, four times for the cylinder with an inner stiff layer or a middle one. From the profiles of the numerical simulations in    Figure 8, it can be clearly identified that there are four, four, and three wave numbers at the inner surface, middle face, and outer surface of the cylinder. It should be mentioned that in the above example analyses, the circumferential instability has not been found. In the practical numerical calculation, the critical strains for axial and circumferential instabilities are obtained by taking wave number ω 1 = 0 and ω 2 = 0 in Eq. 55, respectively. The critical strain for onset of circumferential instibility is found to be much greater than that of axial instability. Therefore, axial instability always occurs earlier than circumferential instability for the cylinder subjected to axial compression. We note that the circumferential instability has been investigated in Zhao et al. (2014) and Xu and Potier-Ferry (2016) for a stiff film on a soft cylinder. In their works, the postbuckling was analyzed and the circumferential instability was found to commonly occur after the formation of axial wrinkling. By further analysis, we can realize that when the cylinder wrinkles, the circumferential compressive stress at the trough of the axial wrinkling turns to be larger and larger with increasing of the axial deformation. Hence, it is not hard to understand that the circumferential instability always starts at the trough of the axial wrinkling. The relative postbuckling analysis for graded elastic cylnders is left for further studies.

CONCLUDING REMARKS
The critical condition for onset of instability of elastic cylinders was analyzed for their Young's modulus and Poisson's ratio varying in the radial direction. When the cylinder is subjected to an axial loading, the radial and circumferential stresses will be produced simultaneously due to the variation of Poisson's ratio, and thus both the axial and the circumferential instabilities may be caused individually. A semi-analytical solution was obtained for the prediction of instability in both the axial and the circumferential directions. The critical strain for onset of the circumferential instability is much greater than that of the axial instability, and consequently the axial wrinkling often occurs before the circumferential instability. For a soft cylinder covered by a bilayer with varied and fixed Poisson's ratios, the numerical results manifest that their critical strains are very close to each other, but their modulus ratios for critical wavelength switching are different obviously. For the soft cylinder with a linearly graded stiff layer, the critical wavelength for the varied Poisson's ratio is slightly smaller than that for the fixed one. Numerical examples show that the critical condition is different for the radially graded elastic cylinder as its Poisson's ratio is assumed to be varied or fixed as a constant. The results obtained may provide valuable reference for the design and evaluation of the relevant cylindrical structures.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
ZW contributed to the conception and design of the study. CZ performed the numerical analysis and prepared the first draft of the manuscript. All authors equally contributed to the manuscript revision, reading, and approved the submitted version.

FUNDING
This work is financed by the National Natural Science Foundation of China (Grant No. 11572108).