Document Type : Research Paper
Authors
^{1} Faculty of Materials and Manufacturing Technologies, Malek Ashtar University of Technology, Tehran, 177415875, Iran.
^{2} Department of Aerospace Engineering, Ferdowsi University of Mashhad, Mashhad, 9177948974, Iran
Abstract
Keywords
Buckling Analysis of Composite GridStiffened Cylindrical Shells Using a Generalized Equivalent Single Layer Theory
^{a }Faculty of Materials and Manufacturing Technologies, Malek Ashtar University of Technology, Tehran, 177415875, Iran
^{b} Department of Aerospace Engineering, Ferdowsi University of Mashhad, Mashhad, 9177948974, Iran
keywords 

ABSTRACT 
Buckling Cylindrical shells Composite Grid structure Reddy’s thirdorder shear theory

In this study, the buckling analysis of multilayer composite closed cylindrical shells, as well as composite grid cylindrical shells, was investigated using a highorder theory modified from Reddy’s thirdorder shear theory under simply supported conditions. The advantage of the present theory compared to other highorder theories is the calculation of the effect of the term of the shell section trapezoidal coefficient on relations related to displacement and strain fields, which improves the accuracy of the results. In grid shells, the discontinuous distribution of stiffness and mass of the shell between the reinforcing ribs and the distance between them is expressed by a suitable distribution function. In the case of integrated and grid cylindrical shells, the validation of the results was performed in comparison with other studies, as well as with the results of the numerical solution obtained using Abaqus software. It is shown that for the first buckling mode, the critical load first increases and then reaches a constant value and for the second buckling mode, the critical load first decreases and then reaches a constant value. Also, in the case of grid shells, by increasing the ratio of the cavity dimension to the dimensions of the whole shell, whether in singlecavity or multicavity mode, the present theory and finite element solution find more difference, indicating the higher accuracy of the present theory for integrated shells. 
Composite structures are used in various industries, including aviation and space industries, because of their high weight resistance and resistance to moisture, corrosion, and other unique properties. Composite shells are one of the most widely used structures and have been considered for many years. One way to improve the performance of composite shells is to use a variety of stiffeners. The nature of the lattice structure changes the direction of destructive loads around the damaged points and thus increases the damage tolerance in these structures. The buckling of lattice shells under various loads, including axial load, is one of the most important modes of their failure. Therefore, studying their buckling behavior is of particular importance. Jaunky et.al [1] proposed a method to calculate the buckling load of composite panels, called the smeared stiffener method. In this method, using mathematical methods, the stiffened plate panels are converted to an equivalent uniform plate panel, which has an anisotropic equivalent stiffness to the original panel. The interaction of the outer shell and the stiffeners was calculated by calculating the stiffness of the ribs and the shell in the joint area. The buckling load is calculated by placing the final equivalent stiffness in the RayleighRitz method. In 2003, Kidane et al. [2] used the same theory to provide an analytical method to examine a composite grid cylinder under buckling loads. The method provided a total critical buckling load. In this method, a single cell of the structure is considered, and the stress and strain matrix and its stiffness are calculated according to the strains of the middle layer of the shell and generalized to the whole structure. Then, the buckling critical load is calculated using the energy method, and finally, the effect of different parameters on the critical load is studied. In 2003, Wodesenbet et al. [3] examined the buckling of an isogrid shell using the same method. In this method, the equations are presented based on the middle plane of the shell, and the critical buckling load for the equivalent shell is calculated by minimizing the total potential energy. The effect of the stiffener stiffness on the stiffness of the whole structure is calculated using the effect of force and the moment of the stiffeners on a single element of the outer shell and integration on the whole structure. Also, using threedimensional modeling in ANSYS finite element software and experimental experiments, the results obtained from the analytical solution were validated. In 2009, Yazdani et al. [4] experimentally investigated the buckling of composite lattice cylinders. In this study, shells with different grids made by a fiber twisting machine were subjected to axial quasistatic loading. According to the test results, the critical load for cylinders with hexagonal and triangular grids was more compared to cylinders with rhombus grids and without grids. Rahimi et al. [5] in 2013 examined the effect of stiffener crosssection profile on the axial buckling load of the composite lattice cylinder shell. In this study, using ANSYS finite element software, peripheral and helix ribs were used to stiff the shell. In the same year, Ghasemi et al. [6], based on the smear method and using the firstorder shear theory, proposed a new method for buckling analysis of composite lattice shells, also taking into account the effect of transverse shear forces of the stiffener. Further, the effect of some geometric parameters on the buckling load was also investigated. In 2014, Kalantari and Fadaei [7] modified the smear method to examine the effect of some geometric parameters on the critical load of the stiffened metal shell. Kidane et.al [2, 3] in 2002 and 2003 presented an analytical method. This method is one of the methods that has been considered by many researchers, and despite the limitations and approximations in this method, it has been referred to many times in recent years. In 2017, Civalek [18] analyzed the buckling of composite panels and shells with different material properties by the discrete singular convolution (DSC) method. He investigated the effects of shell geometric quantities and material properties on buckling and results presented for isotropic, FGM, CNT reinforced FGM, and laminated composite conical and cylindrical shells and panels. In 2018, Shahgholian Ghahfarokhi and Rahimi [19] studied an analytical approach for the global buckling of composite sandwich cylindrical shells with lattice cores. The results show that the proposed approach has high prediction accuracy and low computational cost for the global buckling analysis of composite sandwich shells with lattice cores. Also, another advantage of the developed approach is the ability to predict the global buckling load of stiffened composite cylindrical shells with better accuracy. In 2020, Babaei Mahani et.al [20] investigated the thermal buckling of laminated NanoComposite conical shells reinforced with graphene platelets. They studied the effects of the volume fraction of GPLs, semi vertex angle, lengthtothickness ratio, and other parameters on the critical temperature.
In this study, a modified method for analyzing grid shells is proposed, and the stiffeners will be equivalent to a threelayer shell. These layers, along with the laminate layer of the main shell, will create an equivalent shell, the stiffness matrix of which can be easily calculated. After determining the equivalent stiffness for the shell and rib assembly, the critical load of the axial buckling was calculated using the Ritz method and based on the minimum potential energy principle. Abaqus finite element software was used to validate the results. Examination of the results shows that the proposed method is less different from the results of a finite element compared to previous methods, and it can be claimed that it has higher accuracy than existing methods.
In recent decades, many researchers have analyzed composite shells using the theory of threedimensional elasticity. Given that this theory has an extremely high computational volume when the shell of a material with a high degree of anisotropy is subjected to asymmetric and complex loading, so the analysis of the problem becomes very complex. Therefore, most researchers have analyzed the cylindrical shell using twodimensional and threedimensional highorder shell engineering theories.
To analyze cylindrical shells and derive equilibrium equations, the following hypotheses are first considered:
Figure 1 shows a composite cylindrical shell with medium radius R, thickness h, and length L with reference coordinates (conventional positive directions). The shell middle surface is considered as the reference surface, and the cylindrical coordinate system is installed on a shell head according to Fig. 1 [10, 11].
(1) 

In relation (1), terms u, v, and w are the displacements of the desired point with coordinates (x, j, z) in the multilayer space and the x, j, and z directions, respectively, and t is the time. Parameters u0 and v0 are inplane displacements, and w0 is the outofplane displacement of the desired point (x, j) on the multilayer middle surface. Functions and are the rotations of the line perpendicular to the middle surface of the shell element around the j and xaxes, respectively. Parameters , , , and are highorder terms in the Taylor expansion and represent the transverse deformation modes of the shell and are defined as follows:
(2) 

Fig.1. Composite cylindrical shell and reference coordinates [12]
As mentioned above, in the theory used in this study, the effect of transverse normal strains (along with the thickness) on the surface of the cylindrical shell is considered zero. This assumption makes the calculations of this theory easier, and the results of this theory can be compared with other theories. Similar to Reddy’s theory [9], traction forces on the internal and external surfaces are defined as follows:
(3) 

According to the definition of lateral strains and the placement of Eq. (1) in these definitions and finally zeroing them, a simpler definition of the displacement components is available [9]:
(4) 

(5) 


Finally, the displacement components are summarized as follows [9]:
(6) 

It should be noted that the presence of a trapezoidal coefficient in the main displacement equations, which is considered in the final programming and the present theory, is compared with the results of Reddy’s HOST theory [9], which is a continuation of analytical calculations without taking into account the trapezoidal coefficient.
By defining strains from the theory of linear elasticity for cylindrical shells, the general straindisplacement relations in the cylindrical coordinate system taking into account the effect of the trapezoidal term are as follows [13, 14]:
(7) 

In the present theory, the definition of the strain is more complete than in the present Reddy’s thirdorder shear deformation theory (RTSDT) [9]. This can be considered as one of the factors comparing the two theories [9]. The current RTSDT theory is defined as , which has fewer parameters and sections than the current Reddy’s higherorder shell theory (RHOST) theory. The reason is that in the present theory, is defined instead of .
Assuming that the main axes are defined in the material coordinate system (1,2,3), and the multilayer axes (x, j, z) are defined in the cylindrical coordinate system, the threedimensional stressstrain relationships for the kth orthotropic layer, according to the main axes for the theory that we want to develop, are defined based on different displacement models as follows:


(8) 

In relation (8), stiffness matrix elements are defined as follows [10]:
(9) 

In relation (9), is shear modulus, and and are Young’s modulus and Poisson’s coefficients, respectively, which are not independent of each other and are related to the following relations [10]:
(10) 

The main material axes of every single layer may not correspond to the main coordinate axes (x, j, z). Therefore, it is necessary to transfer the fundamental relations from the singlelayer axes (1,2,3) to the multilayer reference axes. The final relationships are as follows [10]:


(11) 

or in brief:
(12)
The coefficients of the matrix are the elastic constants of the orthotropic material related to the kth layer, according to the arrangement of the layers from the inside (K =1) to the outside (K = NL) [10]:
(13) 

In Eq. (13), s = sinq, c = cosq, and q are the angles of fiber orientation in each layer. The conventional positive direction of angle q is clockwise rotation relative to the reference direction of the fiber angle, i.e., the positive direction of the xaxis. By placing the strain relation (7) in relation (11) and integrating with respect to the thickness, Eq. (12) is reduced to the following equation:
(14)
In Eq. (14):
(15) 

In Eq. (14), the corresponding components of the stress resultant vector and the components of the middle surface strain vector are defined as follows:
(16) 

(17)


As can be seen in Eqs. (16) and (17), stress resultant vector and middle surface strain vector have 30 components, which, given that , , , , , , , , , and are zero, their resultants , , , , , , , , , and will be zero and removed from the equations. The use of these 19 components can be effective in making the results more accurate [15].
In relation (14), the stress resultant vectors for a multilayer with NL number of layers are calculated as follows:


(18) 

In relation (18), z_{i} and z_{i + 1 }are the distances between the inner and outer surfaces of each layer of the middle surface, respectively.
The extraction of governing equations was done using Hamilton’s principle and the principle of virtual work.
The analytical form of Hamilton’s principle can be expressed as follows [14]:
(19) 

In Eq. (19), U is the total energy caused by deformation and is defined as follows:
(20) 

In relation (20), the definition of the shell surface element is as follows [14]:
(21) 

In relation (19), W is the potential energy caused by external forces and is defined as follows:

(22) 
In relation (22), Wx is the work of edge forces at the boundaries, and W0 is the work of surface forces. Also, , and are constant static stresses of the edge, and q is wide pressures on the surface of the shell. In relation (19), K is the kinetic energy and is defined as follows:
(23) 

To solve the problem, Eq. (19) can be written as follows:
(24) 

In the following, each of the terms , and are calculated.
In this section, a brief explanation of the simulation provided for the integrated composite cylinder with Abaqus software for numerical analysis of finite element is provided. For finite element analysis using Abaqus software, a cylindrical shell model was modeled, and, for dynamic analysis, a buckling simply supportedsimply supported boundary conditions were provided. For buckling analysis, proper loading for this axial analysis is also required. Naturally, different dimensions and properties are selected for the model to provide different analyzes. Using a composite laminated layer, the thickness of the composite structure is applied to the structure. Then, using the meshing section, the structure meshed with the S8R element with an average of 17,535 elements, and finally, the present dynamic and buckling analyzes were performed for the structure.
Considering that one of the aims of the present study is the static study of cylindrical shells, in this section, the results of the buckling analysis of composite cylindrical shells under both external pressure loading and compressive axial force are compared with other theories. The result of this analysis is the recognition of the critical limit for the initial stresses applied to the shell. To compare the buckling analysis answers, according to Fig. 2, a composite cylindrical shell with inner radius a, critical radius b, thickness h, and length L is subjected to compressive axial force P and external pressure p.
The simplysimply boundary conditions of the shell are considered. Load interaction parameter S is defined as follows [16]:
(25) 

The concepts of P and p loads are shown in Fig. 2.
Fig. 2. Composite cylindrical shell under the combined loads of external pressure and compressive axial force [16]
Also, for comparison purposes, two parameters and are defined as follows [16]:
(26) 

In Tables 1 and 2, the values of critical buckling loads of P and p, along with the number of corresponding buckling modes obtained from the present study (MRTSDT[1] and RTSDT), were compared with the results of linear buckling analysis using the threedimensional elasticity theory [17], respectively, for S = 1 and S = 5.
Material properties (glass/epoxy) are:
, , , , , , , , . The direction of the fibers in all layers is in the peripheral direction (perpendicular to the axis of the cylinder). The main reason for the difference between the results of the present theory (RHOST) and the exact results of threedimensional elasticity is the hypotheses and approximations used in the highorder theory of the shell.
The effect of different mode numbers (n) on the critical buckling loads in the interaction mode of different loads S (combined loading) is shown in Fig. 3. As the number of environmental modes increases, the overall trend of critical loads decreases. However, the least critical load occurs in the fourth to eighth environmental modes. According to the figure, in the lower modes, the two theories are good compliance, but in the higher modes, the distance between the two theories increases.
Fig. 3. Changes of the critical value of buckling in combined loading in terms of the mode values of different shapes and different interaction coefficients
The effect of thickness to average radius ratio on critical buckling loads is shown in Fig. 4. As the thickness of the structure increases, the overall trend of critical loads decreases. However, the lowest critical load occurs at high thicknesses. In sum, the two theories coincide in all directions, and the differences between the new theory and the predefined theory of the highorder are not apparent except at higher thicknesses.
The effect of the load interaction parameter on the critical buckling load in the laminated layer mode is investigated, and the two buckling modes are shown in Fig. 5. As the effect of the load interaction (S) increases, the overall trend of the critical load moves toward a constant value. According to the figure, in different interaction effects, the two theories are completely coincident; for the first buckling mode, the critical load first increases and then reaches a constant value, and, for the second buckling mode, the critical load first decreases and then reaches a constant value.
Fig. 4. Changes of the critical buckling load in combined loading in terms of thickness to medium radius ratio
Fig. 5. Changes of the critical buckling load in combined loading during different load interactions (S)
As mentioned earlier, due to the presence of ribs and empty spaces in the geometry of the lattice cylinder, this nonuniformity in the stiffness and mass of the lattice shell must be entered into the calculations by a distribution function, which is done by the Heaviside step function. Then, by obtaining a stiffness matrix and a new mass for the lattice shell, natural frequencies were extracted.
It should be noted that due to the lack of similar research, the obtained results were compared with the results of the numerical solution. Before continuing the discussion about the results obtained from the present tables, some points about how to model and analyze the lattice cylinder in finite element software are mentioned. For finite element analysis, using Abaqus software, a lattice cylindrical shell model was modeled with a special method by creating parallel partitions and was prepared for dynamic and buckling analysis with simplysimply boundary conditions (Fig. 6). For buckling analysis, proper axial and peripheral loading is also required. Naturally, different dimensions and properties are selected for the model to provide different analyzes. Using a composite laminated layer, the thickness of the composite structure is applied to the structure. Then, using the meshing section, the structure meshed with the S8R element with an average of 7,980 elements, and finally, the present dynamic and buckling analyzes were performed for the structure.
In the present study, the main geometric properties of the cylindrical shell for the lattice code are described as follows:
(27) 

Fig. 6. Grid composite cylindrical shell modeled in Abaqus software
In the above cases, and are the number of cavities of the mentioned geometric design, in the longitudinal and the circumference direction of the cylinder, respectively. and are the thickness of the ribs present on the cylinder surface in the longitudinal and circumferential directions, respectively, and, then, and are the surface thickness of the cavities on the cylinder surface in both longitudinal and circumferential directions, the amount of which is determined by the rib size.
In the present section, the results of the static analysis on the lattice cylindrical shell are investigated. The geometric properties of the shell are similar to those of frequency analysis. In the following tables, the results of the present generalized theory (MRTSDT) are compared with the results of the present RTSDT theory and the results of the finite element method (FEM). In Table 3, the critical buckling load values of for the three corresponding buckling modes obtained from the present study (MRTSDT and RTSDT) are compared with the results of finite element analysis for different thickness to radius ratios (h/R). Material properties (glass/epoxy) are: , , , , , , , , .
By comparing the present results, a good approximation is established between the results and the finite element results, so that in the thinwalled shell and h/R equal to 0.01, the maximum difference between the present generalized theory and the finite element is 2.57% for the third buckling mode. If in the case of thickwalled shell due to the limitations of the theory and the results, the difference between the generalized theory and the finite element is 15.07% for the third buckling mode, it can be noted that with increasing the h/R ratio, the percentage of the difference between MRTSDT and FEM increases, and also with increasing the number of mode m, the percentage of difference usually decreases. It should be noted that in Table 3, the shape of the modes is plotted for h/R = 0.01. Table 3 presents the critical buckling load values for the first buckling mode and different orthotropic ratios for the whole structure. The difference between the results and the finite element results can be investigated due to large changes in the orthotropic ratio. In the isotropic state and , the difference between the results reached 8.08%, while in the state is 12.84%. It should be noted that the present difference is due to extensive changes in the stiffness of the structure in the present theory and is predictable. As the modulus of elasticity of the composite fibers increases, the stiffness of the structure increases significantly, and, as a result, the critical load of the structure increases significantly. In Table 5, the critical buckling load values for the three corresponding buckling modes obtained from the present study (MRTSDT and RTSDT) are compared with the results of finite element analysis for different orthotropic ratios for longitudinal and transverse ribs . is the modulus of elasticity of the circumferential rib, and is the modulus of the longitudinal rib. The direction of the fibers in all layers is in the peripheral direction (perpendicular to the axis of the cylinder). By comparing the present results, there is a good approximation between the results and the finite element results, so that in the case of the thinwalled shell, the maximum difference between the present generalized theory and the finite element is 2.57% for the third buckling mode. However, in higher , due to the limitations of the theory in terms of stiffness and the results, the difference between the generalized theory and the finite element is 15.53% for the first buckling mode, and is equal to 10.
Comparing the present results, there is a good approximation between the analytical results and the finite element results, such that in the case of the thinwalled shell and h/R equal to 0.001, the difference between the present generalized theory of MRTSDT and the finite element is 0% for all buckling modes. In the case of thickwalled cases, the difference between the generalized theory and the finite element is 10.88% for the third buckling mode.
In Table 6, the critical buckling load values for the three corresponding buckling modes from the present study (MRTSDT and RTSDT) are compared with the results of the FEM finite element analysis for different thickness to radius ratios (h/R). The number of cavities in the present circumferential and longitudinal lattice structure is 20 ( ) and 10 ( ).
Comparing the present results, there is a good approximation between the analytical results and the finite element results, such that in the case of the thinwalled shell and h/R equal to 0.001, the difference between the present generalized theory of MRTSDT and the finite element is 0% for all buckling modes. In the case of thickwalled cases, the difference between the generalized theory and the finite element is 6.42% for the third buckling mode. The reason for this can be explained by the higher density of the structure and its proximity to the integrated structure in the case of smaller networks. The stiffness of the lattice structure in the case of more cavities brings the results closer to an integrated structure in terms of accuracy.
The present MRTSDT theory is more accurate than the RTSDT theory, in which the effect of the trapezoidal term is not taken into account and also does not have higher terms in the definition of . This is especially evident in the comparison of theory with finite element simulations for cylindrical shells.
In the comparison between different highorder theories, the more accurate results of the present theory show a positive effect of the trapezoidal effect, as well as higherorder terms in the definition of .
Table 1. Comparison of the values of critical combined loads of axial ( ) and circumferential buckling ( ) with the number of corresponding buckling modes, in terms of S = 1
S = 1 

RTSDT 
MRTSDT 
L/b = 5 3D Elasticity [17] 

* Difference (%)


Difference (%)


Difference (%)


Difference (%)


(n, m) 


b/a 
5.34 
0.092 
5.40 
0.770 
4.384 
0.091 
4.409 
0.763 
(3, 1) 
0.088 
0.731 
1.03 
14.88 
0.106 
14.87 
0.536 
13.60 
0.105 
13.64 
0.530 
(2, 1) 
0.092 
0.466 
1.05 
13.24 
0.133 
13.23 
0.344 
11.30 
0.131 
11.30 
0.338 
(2, 1) 
0.117 
0.303 
1.1 
12.57 
0.176 
12.61 
0.310 
10.54 
0.173 
10.55 
0.304 
(2, 1) 
0.156 
0.275 
1.15 
12.04 
0.220 
11.99 
0.297 
9.99 
0.216 
9.976 
0.292 
(2, 1) 
0.196 
0.265 
1.2 
11.26 
0.261 
11.30 
0.289 
9.31 
0.257 
9.34 
0.284 
(2, 1) 
0.235 
0.26 
1.25 
* 100% difference × reference [17]/(reference [17]  MRTSDT)
Table 2. Comparison of the values of critical combined loads of axial ( ) and circumferential buckling ( ) with the number of corresponding buckling modes, in terms of S = 5
S = 5 

RTSDT 
MRTSDT 
L/b = 5 3D Elasisity [17] 

* Difference (%)


Difference (%)


Difference (%)


Difference (%)


(n, m) 


b/a 
13.71 
0.380 
13.72 
0.632 
13.07 
0.378 
13.07 
0.628 
(2, 1) 
0.334 
0.556 
1.03 
12.09 
0.335 
12.07 
0.337 
10.87 
0.331 
10.87 
0.334 
(2, 1) 
0.299 
0.301 
1.05 
9.96 
0.420 
9.99 
0.216 
8.13 
0.413 
8.13 
0.213 
(2, 1) 
0.382 
0.197 
1.1 
7.04 
0.506 
7.08 
0.178 
6.14 
0.502 
6.16 
0.176 
(2, 2) 
0.473 
0.166 
1.15 
6.45 
0.525 
6.44 
0.142 
5.40 
0.520 
5.37 
0.140 
(2, 2) 
0.494 
0.133 
1.2 
2.48 
0.514 
2.48 
0.113 
2.47 
0.514 
2.49 
0.113 
(1, 1) 
0.527 
0.116 
1.25 
* 100% difference × reference [17]/(reference [17]  MRTSDT)
Table 3. Comparison of critical buckling load values with the results of the finite element model (Abaqus) in terms of different thickness ratios (h/R) and the number of different buckling modes,
Critical buckling load ( ) * 

n = 3 
n = 2 
n = 1 


Difference (%) 
RTSDT 
MRTSDT 
FEM 
Difference (%)

RTSDT 
MRTSDT 
FEM 
Difference (%) 
RTSDT 
MRTSDT 
FEM 

2.57 
0.158 
0.1591 
0.163 
1.65 
0.785 
0.7861 
0.799 
2.34 
4.423 
4.4238 
4.530 
0.01 
5.21 
5.109 
5.1281 
5.410 
3.31 
9.003 
9.5493 
9.876 
2.38 
4.437 
4.4392 
4.449 
0.1 
8.39 
26.57 
29.423 
32.12 
6.14 
26.96 
28.518 
30.38 
5.39 
9.067 
9.1049 
9.624 
0.2 
12.10 
258.5 
261.34 
297.3 
9.38 
178.9 
180.67 
199.3 
9.31 
24.30 
30.026 
29.00 
0.5 
15.07 
995.5 
1012.6 
1192.3 
11.77 
649.9 
655.51 
743.0 
10.69 
59.00 
67.664 
72.42 
1 



Mode shape 
Table 4. Comparison of critical buckling load values with the results of the finite element model (Abaqus) for different orthotropic ratios ( ), L/R = 2, h/R = 0.01, m = n = 1
Critical buckling load ( ) * 

Difference (%) 
FEM 
MRTSDT 

8.08 
13.756 
12.644 
1 
8.01 
687.45 
635.57 
5 
10.27 
2246.5 
2015.7 
10 
12.84 
3785.1 
3298.8 
40 
Table 5. Comparison of critical buckling load values with the results of the finite element model (Abaqus) for different orthotropic ratios of ribs ( ) and the number of different buckling modes, L/R = 2
Critical buckling load ( ) * 

n =3 
n = 2 
n = 1 


Difference (%) 
MRTSDT 
FEM 
Difference (%) 
MRTSDT 
FEM 
Difference (%) 
MRTSDT 
FEM 

2.57 
0.1591 
0.1633 
1.65 
0.7861 
0.7993 
2.34 
4.4238 
4.53 
1 
5.22 
11.42 
12.05 
4.97 
8.02 
8.44 
5.43 
7.13 
7.54 
2 
8.07 
9.11 
9.91 
9.23 
5.21 
5.74 
10.54 
4.41 
4.93 
3 
12.91 
1.82 
2.09 
11.39 
1.71 
1.93 
12.57 
1.53 
1.75 
4 
13.86 
0.87 
1.01 
13.86 
0.87 
1.01 
15.53 
0.87 
1.03 
10 



Mode shape 
Table 6. Comparison of critical buckling load values with the results of the finite element model (Abaqus) for different thickness ratios (h/R) and the number of different buckling modes
Critical buckling load ( ) * 

n = 3 
n = 2 
n = 1 


Difference (%) 
RTSDT 
MRTSDT 
FEM 
Difference (%) 
RTSDT 
MRTSDT 
FEM 
Difference (%) 
RTSDT 
MRTSDT 
FEM 

0 
0.105 
0.105 
0.105 
0 
0.585 
0.585 
0.585 
0 
0.726 
0.726 
0.726 
0.001 
0 
1.073 
1.075 
1.075 
0 
5.90 
5.90 
5.90 
0 
7.23 
7.23 
7.23 
0.01 
7.90 
11.64 
11.89 
12.91 
0.003 
63.69 
63.90 
63.92 
0.005 
68.50 
68.57 
68.61 
0.1 
8.73 
68.25 
53.93 
59.09 
4.80 
342.44 
351.00 
368.7 
6.61 
219.89 
223.5 
239.3 
0.5 
10.88 
1238.8 
1132.7 
1271.1 
7.68 
40.81 
92.05 
99.71 
7.98 
66.48 
38.47 
41.81 
1 



Mode shape 
One of the differences between the present theory and Reddy’s theory is that in the present theory, the mass and stiffness matrices are symmetric, but in Reddy’s theory, the stiffness matrix is not completely symmetric.
As the h/R ratio increases, the critical buckling load increases, and the difference between the present theory and the finite element results also increases.
In the case of grid shells, it should be noted that by increasing the ratio of the cavity dimension to the dimensions of the whole shell, whether in singlecavity or multicavity mode, the present theory and finite element solution find more difference, indicating the higher accuracy of the present theory for integrated shells.
[1] N. Jaunky, N. F. Knight, D. R. Ambur, 1996. Formulation of an improved smeared stiffener theory for buckling analysis of gridstiffened composite panels, Composites Part B: Engineering, 27(5), 519526.
[2] S. Kidane, G. Li, J. Helms, S.S. Pang, E. Woldesenbet, 2003. Buckling load analysis of grid stiffened composite cylinders, Composites Part B: Engineering, 34(1), 19.
[3] E. Wodesenbet,. S. Kidane, S.S. Pang, 2003. Optimization for buckling loads of grid stiffened composite panels, Composite Structures, 60(2), 159169.
[4] M. Yazdani, H. Rahimi, A. A. Khatibi, S. Hamzeh, 2009. An experimental investigation into the buckling of GFRP stiffened shells under axial loading, Scientific Research and Essay, 9, 914920.
[5] G. H. Rahimi, M. Zandi, S. F. Rasouli, 2013. Analysis of the effect of stiffener profile on buckling strength in composite isogrid stiffened shell under axial loading, Aerospace Science and Technology, 24(1), 198203.
[6] M. Y. M.A.Ghasemi, S.M.Hoseini, 2013.Analysis of effective parameters on the buckling of grid stiffened composite shells based on first order shear deformation theory, Modares Mechanical Engineering, 13(10), 5161. (In Persian)
[7] Fadaei, M. and Kalantari, S., 2014.Stability analysis of weight optimum waffle cylindrical shells A new approach, Modares Mechanical Engineering, 14(14), 177184. (In Persian)
[8] Qatu, M.S., R.W. Sullivan, and W. Wang, 2010, Recent research advances on the dynamic of composite shells, Composite Structures, 93(1), 1431.
[9] Reddy, J. and C. Liu, 1985. A higherorder shear deformation theory of laminated elastic shells. International Journal of Engineering Science 23, 319330.
[10] Reddy, J.N., 2004.Mechanics of laminated composite plates and shells: theory and analysis. CRC press.
[11] Garg, A.K., R.K. Khare, and T. Kant, 2006. Higherorder closedform solutions for free vibration of laminated composite and sandwich shells. Journal of Sandwich Structures and Materials 8, 205235.
[12] Liew, K. and C. Lim, 1996. A higherorder theory for vibration of doubly curved shallow shells. Journal of applied mechanics 63, 587593.
[13] Bert, C.W., 1967.Structural theory for laminated anisotropic elastic shells. Journal of Composite Materials 1, 414423.
[14] Leissa, A. and J.D. Chang, 1996. Elastic deformation of thick, laminated composite shells. Composite structures 35, 153170.
[15] Qatu, A.S, 1967. Structural theory of laminated composite deep thick shells. Composite Materials, 1, 414423.
[16] Jafari. A.A, Bagheri. M, 2006. Free Vibration of Nonuniformly Ring Stiffened Cylindrical Shells Using Analytical, Exprimental and Numerical Methods. ThinWalled Structures Journal, 44, 8290.
[17] Kardomateas. G.A, Philobos. M.S, 1995. Buckling of thick orthotropic cylinderical shells under combined external pressure and axial compresion. AIAA Journal, 33, 19461953.
[18] Civalek, Ö. 2017."Buckling analysis of composite panels and shells with different material properties by discrete singular convolution (DSC) method." Composite Structures, 161, 93110.
[19] Ghahfarokhi, D. S. & G. Rahimi. 2018. "An analytical approach for global buckling of composite sandwich cylindrical shells with lattice cores." International Journal of Solids and Structures, 146, 6979.
[20] Mahani, R. B., et al. 2020. "Thermal buckling of laminated NanoComposite conical shell reinforced with graphene platelets." ThinWalled Structures, 155, 106913.
[1] Modified Reddy’s thirdorder shear deformation theory