Buckling Analysis of Composite Grid-Stiffened Cylindrical Shells Using a Generalized Equivalent Single Layer Theory

1. M. Panahiha ^a, A. Davar ^a*, J. Eskandari Jam ^a, M. Heydari Beni ^a,
2. M. Eskandari Shahraki^b

^a Faculty of Materials and Manufacturing Technologies, Malek Ashtar University of Technology, Tehran, 1774-15875, Iran

^b Department of Aerospace Engineering, Ferdowsi University of Mashhad, Mashhad, 9177948974, Iran

1.  Introduction   

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 Rayleigh-Ritz 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 three-dimensional 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 quasi-static 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 cross-section 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 first-order 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 Nano-Composite conical shells reinforced with graphene platelets. They studied the effects of the volume fraction of GPLs, semi vertex angle, length-to-thickness 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 three-layer 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.

2.    Derivation of Equilibrium Equations Based on a High-Order Theory

In recent decades, many researchers have analyzed composite shells using the theory of three-dimensional 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 two-dimensional and three-dimensional high-order shell engineering theories.

2.1. Hypotheses

To analyze cylindrical shells and derive equilibrium equations, the following hypotheses are first considered:

• Thin-walled shell hypotheses and the first love approximation are used.
• Dimensions of the shell (length and diameter) are limited.
• The damping effect is ignored.
• The material is in the linear elastic range.
• The effect of the nonlinear terms in the equations is not considered.

2.2. Determining the displacement components

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].

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 in-plane displacements, and w0 is the out-of-plane 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 x-axes, respectively. Parameters , , , and  are high-order terms in the Taylor expansion and represent the transverse deformation modes of the shell and are defined as follows:

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:

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]:

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

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.

3.  Defining Strain-Displacement Relations

By defining strains from the theory of linear elasticity for cylindrical shells, the general strain-displacement relations in the cylindrical coordinate system taking into account the effect of the trapezoidal term are as follows [13, 14]:

In the present theory, the definition of the     strain is more complete than in the present Reddy’s third-order 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 higher-order shell theory (RHOST) theory. The reason is that in the present theory,  is defined instead of    .

4.  Definition of Stress-Strain Relationships in Composites and Stress Resultants

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 three-dimensional stress-strain 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:

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

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]:

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 single-layer axes (1,2,3) to the multi-layer reference axes. The final relationships are as follows [10]:

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]:

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 x-axis. 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):

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

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:

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.

5.  Governing Equations

The extraction of governing equations was done using Hamilton’s principle and the principle of virtual work.

6.  Hamilton’s principle

The analytical form of Hamilton’s principle can be expressed as follows [14]:

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

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

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

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:

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

In the following, each of the terms , and are calculated.

7.  Finite Element Modeling and Analysis of Results

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 supported-simply 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.

8.  Results of Buckling Analysis

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 simply-simply boundary conditions of the shell are considered. Load interaction parameter S is defined as follows [16]:

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]:

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 three-dimensional 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 three-dimensional elasticity is the hypotheses and approximations used in the high-order 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 high-order 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)

9.  Validation of Lattice Cylindrical Shell Results

As mentioned earlier, due to the presence of ribs and empty spaces in the geometry of the lattice cylinder, this non-uniformity 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 simply-simply 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:

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.

10. Results of Buckling Analysis

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 thin-walled 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 thick-walled 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 thin-walled 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 thin-walled 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 thick-walled 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 thin-walled 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 thick-walled 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.

11. Conclusions

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 high-order theories, the more accurate results of the present theory show a positive effect of the trapezoidal effect, as well as higher-order 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

* 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

* 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,

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

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

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

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 single-cavity or multi-cavity mode, the present theory and finite element solution find more difference, indicating the higher accuracy of the present theory for integrated shells.

References

[1]          N. Jaunky, N. F. Knight, D. R. Ambur, 1996. Formu-lation of an improved smeared stiffener theory for buckling analysis of grid-stiffened composite pan-els, Composites Part B: Engineering, 27(5), 519-526.

[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), 1-9.

[3]          E. Wodesenbet,. S. Kidane, S.-S. Pang, 2003. Opti-mization for buckling loads of grid stiffened com-posite panels, Composite Structures, 60(2), 159-169.

[4]          M. Yazdani, H. Rahimi, A. A. Khatibi, S. Hamzeh, 2009. An experimental investigation into the buck-ling of GFRP stiffened shells under axial loading, Scientific Research and Essay, 9, 914-920.

[5]          G. H. Rahimi, M. Zandi, S. F. Rasouli, 2013. Analy-sis of the effect of stiffener profile on buckling strength in composite isogrid stiffened shell under axial loading, Aerospace Science and Technology, 24(1), 198-203.

[6]          M. Y. M.A.Ghasemi, S.M.Hoseini, 2013.Analysis of effective parameters on the buckling of grid stiff-ened composite shells based on first order shear deformation theory, Modares Mechanical Engi-neering, 13(10), 51-61. (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), 177-184. (In Persian)

[8]          Qatu, M.S., R.W. Sullivan, and W. Wang, 2010, Re-cent research advances on the dynamic of compo-site shells, Composite Structures, 93(1), 14-31.

[9]          Reddy, J. and C. Liu, 1985. A higher-order shear deformation theory of laminated elastic shells. In-ternational Journal of Engineering Science 23, 319-330.

[10]        Reddy, J.N., 2004.Mechanics of laminated com-posite plates and shells: theory and analysis. CRC press.

[11]        Garg, A.K., R.K. Khare, and T. Kant, 2006. Higher-order closed-form solutions for free vibration of laminated composite and sandwich shells. Journal of Sandwich Structures and Materials 8, 205-235.

[12]        Liew, K. and C. Lim, 1996. A higher-order theory for vibration of doubly curved shallow shells. Journal of applied mechanics 63, 587-593.

[13]        Bert, C.W., 1967.Structural theory for laminated anisotropic elastic shells. Journal of Composite Materials 1, 414-423.

[14]        Leissa, A. and J.-D. Chang, 1996. Elastic defor-mation of thick, laminated composite shells. Com-posite structures 35, 153-170.

[15]        Qatu, A.-S, 1967. Structural theory of laminated composite deep thick shells. Composite Materials, 1, 414-423.

[16]        Jafari. A.A, Bagheri. M, 2006. Free Vibration of Non-uniformly Ring Stiffened Cylindrical Shells Using Analytical, Exprimental and Numerical Methods. Thin-Walled Structures Journal, 44, 82-90.

[17]        Kardomateas. G.A, Philobos. M.S, 1995. Buckling of thick orthotropic cylinderical shells under combined external pressure and axial compre-sion. AIAA Journal, 33, 1946-1953.  

[18]        Civalek, Ö. 2017."Buckling analysis of composite panels and shells with different material proper-ties by discrete singular convolution (DSC) meth-od." Composite Structures, 161, 93-110.

[19]        Ghahfarokhi, D. S. & G. Rahimi. 2018. "An analyti-cal approach for global buckling of composite sandwich cylindrical shells with lattice cores." International Journal of Solids and Structures, 146, 69-79.

[20]        Mahani, R. B., et al. 2020. "Thermal buckling of laminated Nano-Composite conical shell rein-forced with graphene platelets." Thin-Walled Structures, 155, 106913.

[1] Modified Reddy’s third-order shear deformation theory