Design considerations for composite cylindrical shells on elastic foundations subject to compression buckling

Elastic boundary conditions play an important role in the buckling analysis of cylinders under compressive loading. These structures are used widely in aerospace applications and are highly sensitive to geometrical, material, loading, and boundary imperfections. In fact, the presence of these imperfections can lead to catastrophic failure. In 1968, NASA reported relations for obtaining the Knockdown Factor (KDF) based on an empirical method that is valid for isotropic and orthotropic materials; however, these relations do not consider the effect of elastic boundaries that can lead to highly conservative values of KDF. In design practice, a universal KDF of 0.65 has been used for recent designs by NASA, which may not be applicable to new types of structural configuration with different loading and boundary conditions. Therefore, there is a need for robust design factors for future designs which reduce the dependency on testing during preliminary design phases and speeds up the product development process. The availability of up‐to‐date and different KDF expressions for different structural configurations would help engineers to design lighter structures with improved load carrying capacity and reliability. The main objective of this work is to identify the buckling load sensitivity of cylindrical shells due to their boundary conditions and develop KDF relations considering elastic boundaries. To achieve this goal, the effect of axial, radial and tangential support stiffness on a quasi‐isotropic cylinder under axial compression is investigated. A data‐driven design approach is used to develop new KDF empirical relations for a quasi‐isotropic cylinder on different elastic foundations. The accuracy of these relations is within 5% for any elastic foundation considered.


Introduction
Thin-walled cylinders are widely used in the aerospace industry. These structures are well-known to exhibit high load carrying capacity, but also to be prone to buckling and imperfection sensitivity. In fact, deviations from the perfect geometry, material properties and loading and boundary conditions can lead to large reductions of the buckling load in these structures. Therefore the design of thin-walled cylinders under compression loading is usually driven by buckling considerations. The first studies on the buckling strength of isotropic thinwalled cylinders under axial loading were presented by Lorenz [1] (1908), Timoshenko [2] (1910) and Southwell [3] (1914). These solutions were restricted to typical perfect cylinders with simple boundary conditions and assumed a uniform membrane stress state on elastic materials prior to buckling. From the early 1920s to the 1970s several scientists carried out experiments for understanding the buckling behaviour of cylinders under compressive loading. Notable work in this direction was given by Lundquist [4] in 1933, followed by Donnell [5] in 1934. Subsequently, von Kármán and Tsien [6,7] investigated the lack of consistency in the predicted buckling load. In the 1950s and 1960s, Guist [8] collected the data obtained from an extensive campaign of experimental test conducted on many cylinders. In this study the theoretical buckling loads of the cylinders under consideration were shown to be significantly lower than corresponding experimental data. However, it was Koiter [9,10] who quantified the discrepancy between the theoretical buckling load and its corresponding experimental value. Moreover, he established how the actual behaviour of these structures can be significantly affected by the presence of initial imperfections. In 1991, Geier et al. [11] delivered a major contribution to the literature publishing results of experimental buckling tests for many cylinders with different layups confirming the discrepancy between experimental and corresponding classical buckling loads. In conclusion, thin-walled cylinders under compression buckle at loads lower than the classical linear buckling load and this discrepancy is due to the presence of geometrical imperfections (thickness variations, uneven outer surfaces, ply overlaps, dents, cut-outs and voids) [21,26,33]; material imperfections (elastic properties variations); loading imperfections (uneven loading, load eccentricity and loading angle) [27,28] and boundary imperfection (clamped, simply supported, elastic boundaries and variable boundaries) [32,35,14]. Currently, the premise for the pre-design of thin and moderately thick cylinders under compression loading is given by the NASA SP-8007 guidelines [17]. This set of guidelines was published in 1968 and is still used for evaluating the KDF through empirical relations. KDF given in NASA SP-8007 guidelines is defined as the ratio of critical buckling load from the test to the classical buckling load. Put simply, KDF is the ratio of the buckling load of the imperfect to perfect cylinder [23,22]. However, structures designed using this guideline are heavy and conservative leading to high costs and low service-life of the structure itself [18]. Moreover, this guideline does not take into account the full potential of new materials, fabrication process, structural concepts and effect of boundary conditions. In 2007, NASA established the Shell Buckling Knockdown Factor (SBKF) Project. The SBKF Project had the goal of developing and experimentally validating improved (i.e. less-conservative, more robust) analysis-driven shell buckling design factors (KDFs). Moreover, developing design recommendations for launch vehicle structures and determine whether the conservatisms applied to KDFs during the Apollo era were still warranted today. In this context, Hilberger [20] highlighted the need of analysis studying the structural response of launch structures under different design scenarios (e.g. stiffener variation, cylinder geometry, weld lands and boundary conditions).
Recent research has shown high sensitivity of the buckling load to boundary conditions [12][13][14][15][16] (i.e. 50-60% reduction of the buckling load when compared to that predicted with idealised boundary conditions). Launch structures comprise several cylinders welded and/or bolted together (see Fig. 1). Therefore, to predict the buckling load of the single cylindrical component belonging to these structures should consider realistic boundary conditions. This means that a proper evaluation of the experimental buckling load of a single cylinder needs to have a deep control of the stiffness of the surrounding structure used in the experiment [16]. In fact, this elastic foundation can reduce the buckling load significantly from fully clamped conditions to the simply supported in an axial direction and this variation indicates a strong dependency of buckling load on the surrounding structure. Simply put, to consider the real stiffness of the surrounding structure leads to a more reliable evaluation of the buckling behaviour of cylinders under compressive loading. Therefore, the need of rela-tions which take into account the effect of boundary conditions and are used in pre-design phases is of vital importance [20].
The objective of this work is to develop new expressions for knockdown factor (KDF) that include the contribution of elastic boundaries that mimic realistic boundary conditions. These expressions can help engineers and designers during the preliminary design-phase of cylinders used for aerospace applications. More precisely, the focus of the current work is to develop robust analysis-driven KDF for calculating the buckling load of cylinders supported on elastic foundations. These new relations also provide an option to switch from elastic boundaries to different classical boundary conditions. Users can change the radial, axial or tangential stiffness of support to perform hand calculations which would be useful in the preliminary design phase and speed up the design process by reducing the number of finite element (FE) simulations and tests. To achieve this goal, firstly, the results obtained from FE data for classical boundary conditions are correlated with analytical calculations and test data reported in the literature [25,11]. Then, the FE model is modified to include the elastic foundations. The elastic foundations are modelled with spring elements in a cylindrical coordinate system. To explore the large design space Python code was used to generate multiple FE models with spring stiffness varied from 1 N/m (free) to 10 GN/m (fixed). The buckling load of the cylinder was found to change sharply when the spring stiffness is comparable to that of the cylinder. The results also show a significant effect of surrounding structure which may reduce the buckling load by 59% compared to fixed conditions. The paper is organised as follows: Section 2 describes the structural system under consideration and the importance of considering elastic foundations for studying its buckling behaviour. Then, the procedure for obtaining empirical relations of KDF for three different configurations of elastic boundary is given. Section 3 shows numerical results and KDF curves for all cases under consideration. Finally, conclusions are drawn in Section 4.

Cylinders on elastic foundations
Elastic foundations play an important role in the evaluation of the buckling behaviour of cylinders under compression loading. In reality, these structures are supported by surrounding structures, which have finite stiffness. However, these stiffnesses are often not considered in FE simulations of the shell component and boundary conditions are modelled assuming their behaviour as being infinitely rigid (or free to deform). These choices can lead to inaccurate evaluations of the buckling behaviour of cylinders under compression loading. To have a better idea of the role that elastic foundations play, let us consider a cylinder as shown in Fig. 2. The cylinder shown in Fig. 2 is supported on axial springs, which have finite elastic stiffness. Then, by assuming the cylinder comprises a series of vertical strips, which are connected   in series and by neglecting the interaction between them, then the equivalent axial stiffness of the entire system K Eq can be obtained as where K Cyl A and K Spr u represent the axial stiffness of the cylinder and the axial stiffness of springs, respectively. In particular, K Cyl A is the sum of individual axial stiffness of each strip, while K Spr u is the axial stiffnesses for all springs. From Eq. (1), it is noted that the support stiffness plays a crucial role and can alter the value of the stiffness of the entire system. In fact, if the support stiffness has infinitesimal value compared to the cylinder's stiffness, then the overall stiffness of the system significantly reduces and the system buckles at low loads. In contrast, if the support is significantly stiffer than the cylinder, then the overall stiffness of the system tends to the stiffness of the cylinder itself. This means that the cylinder exhibits buckling behaviour similar to classical boundary conditions (i.e. u ¼ 0) and the buckling load increases. Therefore, the buckling behaviour is significantly affected by the value of the stiffness of the support structure. To develop best practice in this area, further study is needed. In this paper, the effect of axial, radial and tangential elastic foundations on the buckling behaviour of a quasi-isotropic (QI) (with inherent yet small degrees of flexural anisotropy) cylinder under compression loading is investigated. In this study, the cylinder is considered to be supported on elastic foundations at the bottom, while is loaded at the top. To introduce the boundary conditions to the top circle a multi-point constraint (MPC) approach is used. All nodes of the top circle (slave nodes) are connected to the centre node (reference point) through pin joints (u ¼ v ¼ w ¼ 0). The reference point is allowed to move in the axial direction (u=free), while the remaining five degrees of freedom are constrained at this reference point leading to a state of pure compressive loading. KDF expressions are deduced for three test cases. Case-1 considers the cylinder to be supported by axial springs and is free to deform in both radial and tangential directions. Case-2 considers the axial displacement to be fully restrained with conventional boundary conditions, while the radial stiffness is elastic and is modelled with springs and the cylinder is free to deform tangentially. In Case-3, the axial and tangential displacements are fully restrained, while the tangential stiffness is elastic and is modelled with springs. By doing so, physical insights regarding the contribution of each elastic component of the boundary stiffness is gained. Subsequently, for each case under consideration an empirical relation of the critical buckling load varying with respect to the stiffness of the elastic support is sought. To obtain these relations, a Python script is used for varying the stiffness of the springs and obtain the corresponding values of buckling loads. Then, these values are used as data for the Curve-Fitting function implemented in Matlab. For all cases under consideration, the shape function best fitting those values of buckling loads was found to be exponential. In order to simplify these functions the stiffness of the cylinder is described using the Equivalent Fully Isotropic Laminate (EFIL) [24] and material invariants. The Young's modulusÊ, shear modulusĜ and Poisson's ratioν of the EFIL can be respectively obtained aŝ where U 1 , U 4 and U 5 are material invariants, detailed in [19].

Validation
At first, using FE analysis, the buckling behaviour of a thin-walled cylinder under axial compression is validated against results given in the literature. In particular, the buckling behaviour of a cylindrical shell under uniform axial compression is studied for both SS1 Uniform compression is applied at the top of the cylinder using a MPC approach. Therefore, all nodes of the top circle (slave nodes) are connected to a control node (see Fig. 5) and a compression force is applied to the control node point. More precisely, connections between the control point and slaves nodes are obtained with pinned constraints. For verification purpose, the geometry and material properties of the cylinder are extracted from Geier et al. [26,11]. The radius and height are R ¼ 250 mm and H ¼ 510 mm, respectively. The ply thickness is 0:125 mm. The material used for the simulation is carbon/epoxy for which GPa and G 23 ¼ 3:40 GPa, respectively. Seven different layups have been considered (see Table 1). The cylinder is meshed with S4R elements using commercial ABAQUS software. For all layups, the mesh convergence study indicates that 320 elements along the circumference and 100 elements along the height are needed for accurate evaluations of mode shapes and buckling loads. For each layup the buckling load is obtained with a linear eigenvalue analysis using the subspace algorithm. Subsequently, these results are compared with experimental results given by Geier et al. [26,11] and/or analytical solutions based on classical laminated-plate theory (CLPT) [26] and first-order shear deformation theory (FSDT) [11] (see Table 1). Table 1 shows that buckling loads under CLPT and FSDT for SS2 boundary conditions match closely with those obtained with Abaqus FE. Under SS1 boundary conditions the experimental results agree well with FE results except for those for the cylinder Z11/Z23. This indicates that the test boundary condition with potting/cascade is close to idealised SS2 or SS1 boundary. In actual application scenarios, cascade dimensions or boundary imperfections could lead to variations in buckling load and further investigation is required for other classical and elastic boundaries, as studied in Section 3.3. It is also interesting to note that the buckling load under SS2 and SS1 conditions are similar. The maximum percentage difference between analytical CLPT and results obtained with Abaqus is 5.50%. Fig. 4 shows a comparison of buckling loads for thin-walled cylinders with layups given in Table 1. Three different boundary conditions are considered. Firstly, the cylinder is subjected to u ¼ 0 (axially constrained), then Fig. 3. Fig. 4 shows that the buckling loads with SS1 and CC1 conditions are nearly equal for all layups. For a given cylinder, the buckling load reaches its maximum value when all translational degrees of freedom (DoF) are constrained and the addition of restrained tangential rotation does not increase the buckling load. This can be explained by considering that the tangential rotation is approximately given by the derivative of radial displacement to axial axis. It is also interesting to note that under only constrained axial translation the buckling load drops by up to 55.32% (for Z17/Z25 layups) compared with CC1 condition; meanwhile for other layups, e.g. Z12/Z24 and Z22, a reduction of 22% of the buckling load is observed.

Effect of elastic foundation and empirical formulation
In this section, the elastic foundation of the cylinder is modelled in FE analysis as a linear extensional spring (SPRING1 element in Aba-qus), which is located at each node of the bottom circumference of cylinder using an in-house built Python script. Subsequently, the spring stiffness is varied from an infinitesimal to a high value and buckling load values obtained in doing so are compared with those obtained using idealised boundaries. Then, the spring stiffness is varied to investigate the effect of axial, radial and tangential elastic foundations on the buckling behaviours of the cylinder. Fig. 5 shows some details of the FE model.
In particular, a QI cylinder with stacking sequence ½0=90=45= À 45=0=90=45= À 45 s made of IM7/8552 carbon-epoxy pre-preg [28,29] is used in this study. The relevant material properties are  Table 2 and Fig. 6 show the buckling loads and mode shapes of the considered cylinder under various boundary conditions, respectively. Fig. 6 shows that the buckling mode shape changes from local to global as more degrees of freedom are constrained on the base of the cylinder. In particular, the cylinder depicted in Fig. 6(a) is axially constrained but allowed to move in radial and tangential directions, which leads to a wavy local buckling mode shape. In Fig. 6(b), the nodes at the bottom of the cylinder are radially constrained in a cylindrical coordinate system, which forces the peak deformation of mode shape upwards. At the same time, the buckling load is higher than that for the axially constrained cylinder, rising from 449 kN to 546 kN (21.6% increase). In Fig. 6(c), the cylinder is also constrained tangentially, which makes the base fully constrained translationally. The addition of tangential support causes the waves to spiral, which shows the coupling effects from shearing and torsion. In the final configuration, the waves are inclined at −45 degrees due to the fact that the −45 degree plies are nearest to the symmetric plane of laminate. When the tangential springs are added, the buckling load changes from 546 kN to 809 kN showing an increase of 35.5% compared to the case of axial and radial constrained displacements. In the following sections, more details of mode shape changes across different classical boundary conditions are investigated by considering the vertical, radial and tangential stiffness of the base as elastic foundations.

Effect of axial elastic foundation
In this case, the stiffnesses of axial springs embedded under each node are varied from 1:0 N/m (u is free or K Su ≈0:0) to 10 GN/m (u ¼ 0 or K Su ¼ 1). The buckling load increases from 333 kN to 449 kN. It is noted that the buckling load does not change further when the spring stiffness is smaller than 1 kN/m or greater than 0:1 GN/m. Hence, it is reasonable to conclude that 1 kN/m represents the free condition while 0:1 GN/m represents simply supported condition in the axial direction.
It is also worth noting that the number of half-waves in Fig. 7 changes from 16 in the case of a soft axial foundation to 18 in the case of a hard axial foundation.
A static analysis for the cylinder under axial load reveals the total axial stiffness of the cylinder to be 0.69 GN/m. In fact, assuming the cylinder comprises 320 parallel strips and noting the axial stiffness of each individual strip is 2.16 MN/m, then the total stiffness is obtained as the sum of the stiffness of each strip. Subsequently, a parametric buckling analysis shows the influence of spring stiffness on buckling load. Results of this parametric analysis are shown in Table 3, which reveals an interesting insight. Table 3 shows a significant variation in buckling load with changes in spring stiffness. In particular, when the axial support stiffness is similarly valued to the axial stiffness Table 1 FE model verification and validation of buckling loads. Nomenclature of the specimens used in the first column is the same as that used in [11,26,27]. Loads are expressed in kN. One test reports two values [26] because it is obtained from two different tests on two cylinders with same layup.      of the cylinder, the buckling load increases significantly. This is due to increasing interaction of the elastic foundation with the cylinder which increases the buckling load. Table 3 also shows that when ratio K Su =K CA varies from 1 to 50, then the buckling load varies most and reaches a maximum value of 449 kN. Finally, buckling loads shown in Table 3 are used to formulate an empirical expression for KDF as follows. The Curve Fitting function of Matlab is used for obtaining an empirical relation for KDF. To obtain this, firstly an exponential function is used because it provides the best fitting function for data shown in Table 3 as follows Fig. 9.
where the constant terms a, b, c and d are defined as c ¼ P umax À P umin ; ð8Þ The axial stiffness of the single strip of the cylinder can be calculated as where t is the thickness of the laminate, while N is the number of strips. Therefore, the axial stiffness of the cylinder can be evaluated as the sum for K CA of all strips or with static analysis. Finally, the empirical expression for KDF can be written as where λ u ¼ d (see Eq. (9)). Finally, the upper and lower bounds of the buckling load (P umax and P umin ) can be evaluated analytically or by FE analysis.

Effect of radial elastic foundation
In this case, the axial direction is simply supported (u ¼ 0 or K Su ¼ 1) and then radial stiffness is varied from 1.0 N/m (w is free or K Sw ≈0:0) to 10 GN/m (w ¼ 0:0 or K Sw ¼ 1). As explained at the beginning of Section 3.3, a localised wavy mode shape is observed where the radial and tangential stiffness is relatively small (Fig. 8a). When the radial spring stiffness is increased to high values, the wavy shape at the bottom disappears and a full circle shape is observed at this position (Fig. 8b). Moreover, when radial stiffness is added to the axial stiffness, the buckling load increases from 449 kN to 546 kN.
From Fig. 10, it is observed that the buckling load starts increasing when the radial stiffness is 0.1 kN/m and attains the maximum value at 1 MN/m. Similar to the axial elastic foundation case, a static analy-    sis was performed that showed the bending stiffness of the cylinder is 3.93 kN/m about the longitudinal axis of cylinder. It is interesting to note that when the radial spring stiffness is comparable to the bending stiffness of the cylinder, the cylinder interacts with the spring and the buckling load changes. More specifically, when the ratio K Sw =K CB reaches 1.0 the KDF w increases and when the ratio is 24.0 no further increases in KDF w are observed. The disappearance of the wavy mode shape confirms this result when radial stiffness is high. The KDF w expressions for the elastic spring in the radial direction can be expressed as where The vertical deflection of the cylinder can be calculated approximately using It is also possible to calculate this stiffness in Abaqus noting K CB ¼ F CB =δ v . Maddux et al. [30] and Gangamwar et al. [34] studied the defection of closed rings and curved beams, respectively. The bending stiffness of composite cylindrical shells can be calculated either from Abaqus static analysis or from Eq. (15) as where I is the second moment of area given by Fig . 11-13 show comparisons between the bending stiffness evaluated with FE analysis and theory for different lengths, thicknesses and radii showing good correlation. When the length of cylinder increases beyond a certain length, such as 1000 mm as shown, in Fig. 11, bending behaviour changes from global to local and the contribution from    the longitudinal direction starts dominating, which introduces some deviation between theory and FE analysis. In the FE model, both ends of the cylinder were constrained in the axial direction (longitudinal to cylinder) and a unit concentrated load was applied at one end and deflection measured to calculate the bending stiffness of the cylinder. The second moment of area was taken about the longitudinal axis of the cylinder and calculated using Eq. (17). Fig. 11 shows the bending stiffness where thickness and radius of the cylinder are 2.096 mm and 250 mm respectively while the height of cylinder varies. Similarly, Fig. 12 shows the bending stiffness where the thickness of cylinder varies and height and radius of cylinder are each 250 mm for all cases. Finally, Fig. 13 shows the bending stiffness for where the radius of cylinder varies and thickness and length of cylinder are 2.096 mm and 250 mm, respectively. It is observed that the thickness of the cylinder greatly influences stiffness. Moreover, the thickness is also important for KDF and critical buckling load calculations. For fully homogenised QI cylinders then the buckling load is directly proportional to the square of the thickness. Note, the maximum buckling load for the current QI layup is 809 kN and buckling load changes from 449 kN to 546 kN when radial stiffness is null and infinite, respectively. Thus, the total buckling load increases due to the radial elastic foundation is approximately 12%.

Effect of tangential elastic foundation
In this case, the cylinder is supported in axial and radial directions as simply supported (u ¼ w ¼ 0 or K Su ¼ K Sw ¼ 1) and the tangential stiffness varies from 1 N/m (v is free or K Sv ≈0:0) to 10 GN/m (v ¼ 0 or K Sv ¼ 1). From the mode-shape it is clear that when tangential displacement is free then the half-waves Fig. 14(a) around the cylinder circumference aligns vertically. When the tangential spring stiffness increases to a high value, these half waves incline at −45 deg as shown in Fig. 14(b) and resemble a twisted cylinder [23]. The mode shape also indicates twisting behaviour and indicates that the shear stiffness is an important consideration when the cylinder is tangentially constrained. Fig. 15 shows that the buckling load increases when the radial stiffness is 1 kN/m and attains the maximum value when it is 0.1 GN/m. When the tangential stiffness of the spring approaches the shear stiffness of the cylinder then the buckling load starts to increase and attains a near maximum value when the tangential spring stiffness is 120.0 times the cylinder's shear stiffness. When the ratio of K Sv =K CG is 1 the KDF v starts to increase until the ratio reaches 120.0 and then no further increases in KDF v are observed. The KDF v expression for tangentially translational elastic foundation in empirical form is where The shear stiffness of the single strip of the cylinder is given by Detailed studies regarding shear stiffness calculations of tubular structure are given by Hoogenboom et al. [31], while shear deflection of beams is documented in [35]. The shear stiffness of the single strip of the cylinder is given by K CG ¼ V=δ ¼ĜA=H, where V is the shear force, G the shear modulus, A is the cross-sectional area and H is the height of the cylinder. The shear stiffness of the cylinder can be calculated either with a static analysis or as sum of the shear stiffness K CG of all strips.
To calculate shear stiffness from the FE model, load was applied tangentially at one end and deformation recorded. The other end of the cylinder was constrained in radial, axial and tangential directions. Comparison of shear stiffnesses calculated from FE and analytical solu-  tions for different cylindrical configurations are shown in Figs. 16-18, showing a good match. In particular, Fig. 16 shows the effect of the variation when height varies, whilst the thickness and radius are 2.096 mm and 250 mm, respectively. Whereas, Fig. 17 shows effects when height and radius of the cylinder are both 250 mm and the thickness of cylinder varies. Similarly, Fig. 18 shows KDF results when thickness and length of cylinder are 2.096 mm and 250 mm respectively and the radius of cylinder varies. Classical buckling load is directly proportional to the modulus of elasticity and the square of the thickness and therefore thickness strongly affects buckling load calculations. However, the maximum buckling load for the QI layup is 809 kN, and buckling load increases from 546 kN to 809 kN when the tangential stiffness is null and infinite, respectively. Thus, the total buckling load increase due to tangential elastic foundation is approximately 32.5%. Initially, when axial support is infinitesimal (K Su ¼ 1 N/m) the buckling load is 333 kN and increases with an increase in axial stiffness of spring and attains a maximum value of 449E kN when support is large (K Su ¼ 10 GN/m). When the radial stiffness is added then the buckling load changes from 449 kN (K Sw ¼ 1 N/m) to 546 kN (K Sw ¼ 10 GN/m). Finally, tangential stiffness is added to current configurations and found further increases in buckling load from 546E kN (K Sv ¼ 1 N/m) to 809 kN (K Sv ¼ 10 GN/m). Fig. 19 shows that the buckling load increase due to the axial elastic foundation is approximately 11-12% compared to the highest possible critical buckling load of 861 kN (homogeneous quasi-isotropic lay-up). Similarly, the buckling load increase due to the radial elastic foundation is also about 11-12%. However, the buckling load due to tangential elastic foundations is 32-35%. Thus, the tangential elastic foundation contributes most to the critical buckling load when the axial and radial displacements are eliminated. It is also     19. Comparison of FE and empirical critical buckling load for axial, radial and tangential elastic foundation. For the red curve, the stiffnesses of axial springs embedded under each node vary from 1:0 N/m (u is free or K Su ≈0:0) to 10 GN/m (u ¼ 0 or K Su ¼ 1). For the purple curve, the axial direction is simply supported (u ¼ 0 or K Su ¼ 1) and then radial stiffness varies from 1.0 N/m (w is free or K Sw ≈0:0) to 10 GN/m (w ¼ 0:0 or K Sw ¼ 1). For the blue curve, the cylinder is supported in axial and radial directions as simply supported (u ¼ w ¼ 0 or K Su ¼ K Sw ¼ 1) and the tangential stiffness varies from 1 N/m (v is free or K Sv ≈0:0) to 10 GN/m (v ¼ 0 or K Sv ¼ 1).

Comparison of axial, radial and tangential elastic foundation
shown that the buckling load at the right end of the axial (red) curve is the same as the left end of the radial (pink) curve. This means that P umax ¼ P wmin . Similarly, the critical buckling load at the right end of the radial (pink) curve is the same as that at the left in of the tangential (blue) curve. Therefore, P wmax ¼ P vmin . Finally, the right end of the tangential curve indicates the highest buckling load and this load approaches the classical buckling load and becomes equal to P cl ¼ 1:2πEt 2 [8,23], when there is no anisotropy of the homogenised layup.
Since the classical buckling load of QI laminates is a quadratic function of thickness, the proposed empirical relations are validated for different laminate thicknesses in the following example. Fig. 20 shows the comparison of buckling load from FE and empirical relations for 3 mm laminate thickness. The trends indicate good correlations between FE and empirical buckling loads.
Similarly, these KDF relations are validated for 4 mm and 5 mm thickness and are in good agreement as shown in Figs. 21 and 22, where the change in buckling load between corresponding curves across figures confirms the quadratic dependency on thickness. It is observed that even as thickness increases, the major contributor in Fig. 20. Comparison of axial, radial and tangential elastic foundations (FE and empirical) on 3 mm thick cylindrical shell. For the red curve, the stiffnesses of axial springs embedded under each node vary from 1:0 N/m (u is free or K Su ≈0:0) to 10 GN/m (u ¼ 0 or K Su ¼ 1). For the purple curve, the axial direction is simply supported (u ¼ 0 or K Su ¼ 1) and then radial stiffness varies from 1.0 N/m (w is free or K Sw ≈0:0) to 10 GN/m (w ¼ 0:0 or K Sw ¼ 1). For the blue curve, the cylinder is supported in axial and radial directions as simply supported (u ¼ w ¼ 0 or K Su ¼ K Sw ¼ 1) and the tangential stiffness varies from 1 N/m (v is free or K Sv ≈0:0) to 10 GN/m (v ¼ 0 or K Sv ¼ 1). Fig. 21. Comparison of axial, radial and tangential elastic foundations (FE and empirical) on 4 mm thick cylindrical shell. For the red curve, the stiffnesses of axial springs embedded under each node vary from 1:0 N/m (u is free or K Su ≈0:0) to 10 GN/m (u ¼ 0 or K Su ¼ 1). For the purple curve, the axial direction is simply supported (u ¼ 0 or K Su ¼ 1) and then radial stiffness varies from 1.0 N/m (w is free or K Sw ≈0:0) to 10 GN/m (w ¼ 0:0 or K Sw ¼ 1). For the blue curve, the cylinder is supported in axial and radial directions as simply supported (u ¼ w ¼ 0 or K Su ¼ K Sw ¼ 1) and the tangential stiffness varies from 1 N/m (v is free or K Sv ≈0:0) to 10 GN/m (v ¼ 0 or K Sv ¼ 1). Fig. 22. Comparison of axial, radial and tangential elastic foundations (FE and empirical) on 5 mm thick cylindrical shell. For the red curve, the stiffnesses of axial springs embedded under each node vary from 1:0 N/m (u is free or K Su ≈0:0) to 10 GN/m (u ¼ 0 or K Su ¼ 1). For the purple curve, the axial direction is simply supported (u ¼ 0 or K Su ¼ 1) and then radial stiffness varies from 1.0 N/m (w is free or K Sw ≈0:0) to 10 GN/m (w ¼ 0:0 or K Sw ¼ 1). For the blue curve, the cylinder is supported in axial and radial directions as simply supported (u ¼ w ¼ 0 or K Su ¼ K Sw ¼ 1) and the tangential stiffness varies from 1 N/m (v is free or K Sv ≈0:0) to 10 GN/m (v ¼ 0 or K Sv ¼ 1).  buckling load is tangential stiffness when axial and radial stiffnesses are infinite. Subsequently, a similar cylinder with a fully homogenised QI layup is studied with stacking sequence ½60= À 60 2 =0 3 =60 2 =0= À60=60 2 = À 60 3 =0 2 =60 s . The thickness of the laminate is t ¼ 2:096 mm, while material properties are the same as those used for the previous layup (½0=90=45= À 45=0=90=45= À 45 s ). This new configuration eliminates all anisotropies [23] and, as a result, has the highest possible buckling load. Fig. 23 shows a comparison between the buckling load of both layups for an axial elastic foundation. Similarly, Fig. 24 shows the comparison for radial elastic foundations for both laminates. Finally, Fig. 25 compares the tangential elastic foundation. In all cases, good agreement is observed between results obtained from FE analysis and those calculated with the empirical expressions. The ratio of highest buckling load for QI to homogenised QI with zero anisotropy also known as the ideal case is 0.94 and this slight reduction reflects the relatively small amount of anisotropy present in the original QI layup. The maximum buckling load for the homogeneous, with zero flexural/twist coupling, or isotropic layup can be calculated as P cl ¼ 1:2πEt 2 , which gives a buckling load of 861 kN and is in good agreement with results shown in Fig. 25 (see black curves).
To further validate our KDF expressions, three different composite layups with relatively small amounts of anisotropy is taken from Weaver [23] as shown in Table 4. Layup A, F, and B are fully homogenised QI layups that contains a small amount of flexural/twist anisotropy. The magnitudes of this anisotropy is quantified as which, for a homogeneous, QI, symmetric layup (D 11 ¼ D 22 ) becomes To further validate the KDF relations, a similar cylinder with different stacking sequences is validated from FE, and comparisons of buckling loads are shown in Fig. 26. It is interesting to note that, the presence of small amounts of flexural/twist anisotropy does not influence the axial and radial elastic foundations. In other words, a small amount of flexural/twist anisotropy has a negligible effect on the upper and lower bound of axial and radial foundations. However, there is a small influence on the upper bound of tangential elastic foundation on critical buckling load and presumably is due to the shear induced by flexural/twist coupling at the boundary which can interact with the tangential stiffness of the elastic foundation. For larger degrees of flexural/twist coupling or where layups exhibit other types of anisotropies further work is required to characterise responses.   Table 4 Quasi-Isotropic Laminates with relatively small amounts of flexural/twist anisotropy [23].

Conclusions
In this study, the effects of various types of elastic foundation on buckling behaviour of quasi-isotropic cylindrical shells without geometrical imperfections have been studied. Current focus has been on the effect of boundary support conditions on linear buckling loads. The deleterious effect of geometric imperfections on buckling is well understood but the effect of variation in elastic support is not so advanced. Therefore this study addresses this deficiency before future work will examine the combined effect of geometric imperfection and elastic support conditions on buckling response. Different spring types are used to obtain the buckling load and mode shapes of a quasiisotropic cylindrical shell under various boundary conditions. The buckling load of a thin-walled cylinder was found to reduce by up to 59% due to the presence of elastic boundary conditions. It is also interesting to note that the buckling load of a cylindrical shell changes as support stiffness becomes comparable in magnitude to the cylinder stiffness. Since the structure is mostly surrounded by finite stiffness of other structural components, there is a need for robust KDFs representing elastic foundations, which can help predict real buckling loads. Robust KDF expressions for axial, radial and tangential elastic foundation are proposed in this study to calculate buckling loads of cylinders supported on different types of elastic foundations. These relations can be useful in the preliminary design phase and speed up the design process by reducing the number of FE simulations and testing. It is also useful to investigate transitions of buckling mode shape and critical load among different classical boundaries.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.