A simple spline finite strip for buckling analysis of composite cylindrical panel with cutout

Abstract

In this study, a new sub-parametric strip element is developed to simulate the axially loaded composite cylindrical panel with arbitrary cutout. For this purpose, a code called SSFSM is developed in FORTRAN to analyze the buckling of panels. The first order shear deformation theory is used to form the strain-displacement relations. Spline and Lagrangian functions are used to derive element shape functions in the longitudinal and transverse directions, respectively. The computational cost of the SSFSM is decreased dramatically, as mapping functions of the strip element are very simple. The results obtained from the SSFSM are compared with those of the literature and the results obtained by ABAQUS to show the validity of the proposed approach. A parametric study is performed to show the capability of the SSFSM in calculating the panel buckling load. Results indicate that increasing the panel thickness and panel central angles cause an increase in panel buckling load. The cutout shape is an important factor influencing the panel buckling load. For instance, when the angle between the direction of big chord of the elliptical cutout and compressive load direction are 0 and 90 degrees, the panel buckling load reaches its minimum and maximum magnitude, respectively.

Keywords:
Composite; perforated cylindrical panel; spline finite strip; sub-parametric

1 INTRODUCTION

In the recent decades, plates and cylindrical composite panels with and without holes are used widely in the aerospace, civil and mechanical industries (Shuohui et al. 2016a-b, Nejati et al. 2017, Asadi and Qatu 2013, Khalfi et al. 2014, Khetir et al. 2017, Asadi et al. 2012). In this matter, buckling analysis of these structures under different load cases has been investigated by several researchers using different methods (Tiantang et al. 2017, Tiantang et al. 2016, Patel et al. 2006, Vuong and Dung 2017, Soltani et al. 2019, Asadi and Qatu 2013). Magnucki and Mackiewicz (2006) investigated the effects of panel boundary conditions on the elastic buckling of composite cylindrical panel. In their research, an analytical approach was used to calculate the buckling load of a panel without any cutout. Effects of rectangular and circular cutouts on the buckling load of laminated plate were investigated analytically by Masia et al (2005) and Ashrafi et al (2013), respectively. Although analytical approaches are simple to calculate the buckling load, but they have some limitations such as cutout shape, panel boundary conditions and the structure shape. To overcome these limitations, the Finite Element Method (FEM) was used extensively to calculate the buckling load of composite cylindrical panels. Anil et al (2007) evaluated the effect of rectangular cutout on the critical load of composite plate using FEM. The influence of circular cutout on the buckling load of laminated plate was investigated by Komur and Sonmez (Komur and Sonmez 2008). The effects of cutouts with various shapes (such as circular and elliptical shapes) on the buckling load of composite plates were also evaluated by FEM in the several researches (Arbelo et al. 2015, Shuohui et al. 2015, Thinh et al. 2016, Rajanna et al. 2016, Soni et al. 2013, Pascan et al. 2012).

Use of the Finite Strip Method (FSM) is simpler and more economical compared to the FEM, in order to analyze long thin-walled structures such as cylindrical panels (Mirzaei et al. 2015, Fazilati 2017, Balogh and Logo 2015, Rondal 1998). Moreover, while in the FEM the structure is discretized in transverse and longitudinal directions, it is only discretized in transverse direction in FSM (Rondal 1998). FSM can be divided into two categories. In the first one, called semi-analytical FSM, the longitudinal shape function of strip element is considered as harmonic function. Also, the transverse shape function of strip element is considered as polynomial function. Buckling load of laminated plates without any cutout subjected to different load cases was evaluated using semi-analytical FSM (Yao et al. 2015, Ovesy et al. 2015, Dawe et al. 1995, Ovesy and Fazilati 2012b). Semi- analytical FSM is very useful and efficient for buckling analysis of laminated plate when the boundary conditions of both ends of strip elements are simply supported (Wang et al. 1998, Dawe and Wang 1994, Cheung and Tham 1997). The longitudinal harmonic shape functions of strip element are complicated in case of more complex plate/panel boundary conditions (Wang et al. 1998, Dawe and Wang 1994, Cheung and Tham 1997). Spline FSM (SFSM) is used in the analysis of buckling load of plate/panel with different boundary conditions (Second category of the FSM) (Wang et al. 1998, Dawe and Wang 1994). The spline function is used as the shape function of strip element in the longitudinal direction. Au and Cheung (1996) investigated the free vibration and stability of shell without any cutout utilizing iso-parametric SFSM. In their work, First Order Shear Deformation plate Theory (FSDT) was used to model the shell. Fazilati and Ovesy (2013) investigated the effects of structure boundary conditions and material on the buckling load of laminated stiffened cylindrical panel with rectangular cutout. They developed a 3D model of cylindrical panel by using SFSM based on Higher order Shear Deformation Theory (HSDT). Finally, Eccher et al. (2008) evaluated the elastic buckling of perforated thin-walled structures such as cylindrical panels by iso-parametric SFSM. In their research, an iso-parametric spline finite strip element was developed based on the FSDT. The strip element shape functions of longitudinal and transverse directions have been considered as polynomial and B3 spline respectively (Eccher et al. 2008). Based on the literature, the order of mapping function of sub-parametric element is lower than that of iso-parametric element (Akin 2005) and therefore, equilibrium equations obtained in case of sub-parametric element have lower order than iso-parametric ones and could cause lower computational cost (Akin 2005).

A review of the literature reveals that there are two major limitations in buckling analysis of perforated laminated cylindrical panel. While the computational cost of sub-parametric element is lower than that of the iso-parametric element, there is no sign of considering sub-parametric SFSM in buckling analysis of perforated laminated cylindrical panel. Moreover, although the effects of single cutout on the buckling load of cylindrical panel are investigated in the literature, the effects of cutout with different shapes on the buckling load of composite laminated cylindrical panel are not compared with each other. Therefore, this paper is a response to these needs. In this regard, at first the shape and mapping functions of sub-parametric spline finite strip element are developed. Next, the displacement functions are derived. Then, the components of Jacobian matrix are calculated and based on the FSDT, the strain-displacement and stress-strain relationships are extended. According to the shape functions, the load vector of strip element is evaluated. Then, components of elastic stiffness and stability stiffness matrices of sub-parametric element are derived. Finally, according to the developed element, the elastic buckling analysis of laminated perforated cylindrical panel is carried out. To validate the results obtained in this paper, comparison are made between the results obtained in this study and those published in the literature and those obtained from commercial finite element software, namely ABAQUS (ABAQUS 2006).

2 Methodology

The cylindrical panel with arbitrary cutout is shown in Figure 1. As it is shown in Figure 1, R, B, $\varnothing$ and L represent the panel radius, width, central angle and length, respectively. The x, y and z axes are considered in the longitudinal, circumferential and radial directions of the panel, respectively. A local coordinate system is established for the panel using x, y and z axes. Based on the first order shear deformation theory and local coordinate system considered in Figure 1, five degrees of freedom ($u,v,w,{\psi }_x,{\psi }_y$) are considered for each node of the panel (Moradi et al. 2011). These degrees of freedom are shown in Figure 1. $u$ and$v$ are membrane (in-plane) displacement in x and y directions, w is the panel deflection in radial direction (out of plane), and ${\psi }_y$ and ${\psi }_x$ are the strip edge rotations about the x and y directions, respectively. According to the load vector and linear stiffness matrix of the cylindrical panel, the pre-buckling displacement vector of the panel is established. Based on the panel pre-buckling displacement vector, membrane force vector of the panel is calculated. Also, using the panel membrane force vector, the geometrical stiffness matrix is calculated. Finally, the buckling load of the cylindrical panel through the solution of linearized eigenvalue problem is calculated. A finite strip model of perforated cylindrical panel is shown in Figure 2. The process described above was carried out in a FORTAN code called Sub-parametric Spline Finite Strip Method (SSFSM) developed in this paper. All stages performed in this code are shown in Figure 3. The details of the process are described in the following sections. In this study, a sub-parametric strip element is developed to calculate the cylindrical panel axial buckling load. In this regard, some features of strip element such as shape functions, mapping functions, Jacobian matrix and stress-strain relationships are established. Based on these features, the stiffness matrix and load vector of the strip element are developed. The stiffness matrices and load vectors of strip elements are used in assemblage of global stiffness matrix and load vector of the cylindrical panel, respectively. In the tenth step of the algorithm, the calculation of ${\left[K_E\right]}^{-1}$ for large number of degrees of freedom is very time-consuming and in the current analysis, $\left[K_E\right]$ is decomposed by the modified Cholsky method as $\left[K_E\right]$=${\left[U\right]}^T\left[D\right]\left[U\right]$, in which [U] is an upper triangular matrix and [D] is a diagonal matrix (Rakotonirina 2013). After that, the pre-buckling deformations [d₀] can be extracted. Also, in the step 13^th, the linear eigenvalue problem is solved by the inverse iteration method (Van Ness 1969) which taking into account $\left[K_E\right]$ and $\left[K_G\right]$ as symmetric banded matrices. Several eigenvalues and their eigenvectors can be extracted using the aforementioned method.

Figure 1
View of cylindrical panel.

Figure 2
Finite strip model of cylindrical panel with a circular cutout.

Figure 3
SSFSM algorithm.

2.1 Shape functions of strip element

The $B_3$ spline (Moradi et al. 2011) and four-noded Lagrangian functions (Moradi et al. 2011) are used to estimate the strip element displacement functions in the longitudinal and transverse directions, respectively. As it is indicated in Figure 4, the $\eta$ coordinate is considered in the transverse direction (y direction) of the strip element. The shape functions of strip element in the transverse direction are presented in equation 1. In this equation, $L_j\left(\eta \right)$ is the Lagrangian function for jth nodal line in $\eta$ coordinate. $\eta$ is varied between -1 and +1. As it is presented in equation 1, the shape functions of strip element in the transverse direction are third degree polynomial function.

$\begin{array}{l}{L}_1\left(\eta \right)=-\frac{9}{16}\left({\eta }^3-{\eta }^2-\frac{1}{9}\eta +\frac{1}{9}\right)\\ L_2\left(\eta \right)=\frac{27}{16}\left({\eta }^3-\frac{1}{3}{\eta }^2-\eta +\frac{1}{3}\right)\\ L_3\left(\eta \right)=-\frac{27}{16}\left({\eta }^3+\frac{1}{3}{\eta }^2-\eta -\frac{1}{3}\right)\\ L_4\left(\eta \right)=\frac{9}{16}\left({\eta }^3+{\eta }^2-\frac{1}{9}\eta -\frac{1}{9}\right)\end{array}$ (1)

The C⁰ continuity is implied in the strip transverse direction (Eccher et al. 2008). Due to the use of this shape function, it is possible to impose any arbitrary boundary condition on the longitudinal strip edges. Four nodal lines are considered in the longitudinal direction of the strip element and m knots (node) are considered in each nodal line. In other word, each nodal line is divided into m sections. The $B_3$ spline function is used to derive the strip element displacement function in the longitudinal direction as shown in equation 2.

${\varphi }_i\left(\xi \right)=\frac{1}{6}\left\{\begin{array}{c}\genfrac{}{}{0pt}{}{0\xi <{\xi }_{i-2}}{{(\xi -{\xi }_{i-2})}^3{\xi }_{i-2}\le \xi <{\xi }_{i-1}}\\ \genfrac{}{}{0pt}{}{1+3\left(\xi -{\xi }_{i-1}\right)+3{(\xi -{\xi }_{i-1})}^2-3{\left(\xi -{\xi }_{i-1}\right)}^3{\xi }_{i-1}\le \xi <{\xi }_i}{1+3\left({\xi }_{i+1}-\xi \right)+3{\left({\xi }_{i+1}-\xi \right)}^2-3{\left({\xi }_{i+1}-\xi \right)}^3{\xi }_i\le \xi <{\xi }_{i+1}}\\ \genfrac{}{}{0pt}{}{{({\xi }_{i+2}-\xi )}^3{\xi }_{i+1}\le \xi <{\xi }_{i+2}}{0\xi >{\xi }_{i+2}}\end{array}\right.$ (2)

where, ${\mathrm{\varphi }}_{\mathrm{i}}\left(\mathrm{\xi }\right)$ is the spline function for the ith knot in $\mathrm{\xi }$ coordinate, while $\mathrm{\xi }$ is varied between -1 to m+1. In Figure 5a, a ${\mathrm{B}}_3$ spline function is drawn. For each knot (node) considered on each nodal line, a spline function is assigned. In Figure 5b, all of the spline functions considered on a nodal line are drawn. Two virtual nodes are considered outside the panel (nodes -1 and m+1 in Figure 5b) for each nodal line to apply any arbitrary boundary condition on the transverse edges of the strip element (Moradi et al. 2011). The displacement functions of the strip element are derived based on the shape functions considered in the transverse and longitudinal directions of the element (see equation 3). The displacement function presented in the equation 3, is derived based on the FSDT (Khayat et al. 2016a). In equation 3, $\mathrm{Z}\left(\mathrm{\xi },\mathrm{\eta }\right)$ is the displacement of strip element at degree of freedom Z in coordinates $\mathrm{\xi }$ and $\mathrm{\eta }$, j is the number of nodal line and i is the number of knot (node) considered on the j^th nodal line,${\mathrm{L}}_j\left(\mathrm{\eta }\right)$ is the Lagrangian function for the j^th nodal line, ${\mathrm{\varphi }}_i\left(\mathrm{\xi }\right)$ is the spline function for i^th knot, and ${\mathrm{Z}}_{ij}$ is the unknown multiplier for the j^th nodal line, i^th knot and degree of freedom $\mathrm{Z}$. Z could be considered as $\mathrm{u},\mathrm{v},\mathrm{w},{\mathrm{\psi }}_{\mathrm{x}}$or ${\mathrm{\psi }}_{\mathrm{y}}$.

$\begin{array}{l}u\left(\xi ,\eta \right)={\displaystyle \sum }_{j=1}^4{\displaystyle \sum }_{i=-1}^{m+1}{L}_j\left(\eta \right)\times {\varphi }_i\left(\xi \right)\times u_{ij}\\ v\left(\xi ,\eta \right)={\displaystyle \sum }_{j=1}^4{\displaystyle \sum }_{i=-1}^{m+1}{L}_j\left(\eta \right)\times {\varphi }_i\left(\xi \right)\times v_{ij}\\ w\left(\xi ,\eta \right)={\displaystyle \sum }_{j=1}^4{\displaystyle \sum }_{i=-1}^{m+1}{L}_j\left(\eta \right)\times {\varphi }_i\left(\xi \right)\times w_{ij}\\ {\psi }_x\left(\xi ,\eta \right)={\displaystyle \sum }_{j=1}^4{\displaystyle \sum }_{i=-1}^{m+1}{L}_j\left(\eta \right)\times {\varphi }_i\left(\xi \right)\times {\psi }_x{}_{ij}\\ {\psi }_y\left(\xi ,\eta \right)={\displaystyle \sum }_{j=1}^4{\displaystyle \sum }_{i=-1}^{m+1}{L}_j\left(\eta \right)\times {\varphi }_i\left(\xi \right)\times {\psi }_y{}_{ij}\end{array}$ (3)

Figure 4
Lagrangian function in transverse direction.

The displacement function of the strip element is defined based on the $4\times (\mathrm{m}+3)$ knots. equation 3 can be expressed as equation 4 in which $\left\{\mathrm{\Delta }\right\}$ is strip element displacement vector, ${\mathrm{\Delta }}_{\mathrm{i}\mathrm{j}}$ is the unknown multiplier for ith knot of jth nodal line, ${\mathrm{\varphi }}_{\mathrm{i}}$ is the ith knot spline function and ${\mathrm{L}}_{\mathrm{j}}$ is the jth nodal line Lagrangian function. $\left\{\mathrm{\Delta }\right\}$ is comprised of $\left\{\mathrm{u}\right\}$, $\left\{\mathrm{v}\right\}$, $\left\{\mathrm{w}\right\}$, $\left\{{\mathrm{\psi }}_{\mathrm{x}}\right\}$ and $\left\{{\mathrm{\psi }}_{\mathrm{y}}\right\}$ (see equation 5).

Figure 5
Spline function in longitudinal direction (a) $B_3$ spline function (b) Spline functions for a nodal line.

$\left\{\mathrm{\Delta }\right\}=\sum _{i=-1}^{m+1}\sum _{j=1}^4{\mathrm{\Delta }}_{ij}{\varphi }_i{L}_j$ (4)

$\left\{\mathrm{\Delta }\right\}={\left\{\begin{array}{ccc}u& v& \begin{array}{ccc}w& {\psi }_x& {\psi }_y\end{array}\end{array}\right\}}^T$ (5)

2.2 Mapping functions

The sub-parametric strip element in the local Cartesian coordinate X and Y is shown in Figure 6a. This element is mapped into the parent element as it is shown in Figure 6b. The mapping functions are used to map the element from (X,Y) coordinate to the ($\mathrm{\xi },\mathrm{\eta }$) coordinate system. The strip element considered in this paper has two advantages. The longitudinal edges of the strip element are linear and parallel to each other and therefore, the width of the strip (b) is constant in the longitudinal direction of the strip element. On the other hand, the transverse edges of strip element (start and end edges of the element) are defined as curved function and consequently any arbitrary cutout could be considered in the cylindrical panel. According to these advantages, the mapping functions of the strip element are estimated using eight nodes on the start and end edges of the strip element. These nodes are shown in Figure 6. In fact, the mapping functions of the element are derived based on the eight nodes while the displacement functions of the element are estimated based on the $4\times (\mathrm{m}+3)$ nodes. In other word, this element is a sub-parametric element (Eccher et al. 2008). x is estimated with the x-coordinate of eight nodes as equation 6. In this equation, ${\mathrm{N}}_{\mathrm{i}}$ is the mapping function of the ith node and ${\mathrm{x}}_{\mathrm{i}}$ is the x-coordinate of the ith node.

$x=\sum _{i=1}^8{N}_i{x}_i$ (6)

Figure 6
Strip element (a) actual element (b) parent element.

The mapping functions ${\mathrm{N}}_1$ to ${\mathrm{N}}_8$ are written as equation 7. The magnitude of ${\mathrm{N}}_{\mathrm{i}}$ is 1 at the ith node and 0 at the other nodes. In equation 7, m is the number of knots considered on the nodal line.

$\begin{array}{l}{N}_1=\frac{9}{16m}\times \left(\xi -m\right)\times \left(\eta -1\right)\times \left({\eta }^2-\frac{1}{9}\right)\\ N_2=-\frac{9}{16m}\times \left(\xi -m\right)\times \left(\eta +1\right)\times \left({\eta }^2-\frac{1}{9}\right)\\ N_3=\frac{9}{16m}\times \xi \times \left(\eta +1\right)\times \left({\eta }^2-\frac{1}{9}\right)\\ N_4=-\frac{9}{16m}\times \xi \times \left(\eta -1\right)\times \left({\eta }^2-\frac{1}{9}\right)\\ N_5=-\frac{27}{16m}\times \left(\xi -m\right)\times \left({\eta }^2-1\right)\times \left(\eta -\frac{1}{3}\right)\\ N_6=\frac{27}{16m}\times \left(\xi -m\right)\times \left({\eta }^2-1\right)\times \left(\eta +\frac{1}{3}\right)\\ N_7=-\frac{27}{16m}\times \xi \times \left({\eta }^2-1\right)\times \left(\eta +\frac{1}{3}\right)\\ N_8=\frac{27}{16m}\times \xi \times \left({\eta }^2-1\right)\times \left(\eta -\frac{1}{3}\right)\end{array}$ (7)

${\mathrm{N}}_{\mathrm{i}}$ is simplified as the multiplication of two functions, ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\xi }}$ and ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\eta }}$ where ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\xi }}$ and ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\eta }}$ are functions of $\mathrm{\xi }$ and $\mathrm{\eta }$, respectively (see equation 8). The constant multiplier related to the m is considered as a part of ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\eta }}$.

$N_i=N_i^{\xi }\times N_i^{\eta }$ (8)

$N_i^{\eta }$ shape functions are cubic Lagrangian fuctions and are used to approximate boundary curve of each strip. So, each closed boundary of cutout is estimated by the several cubic polynomials. As it is illustrated in equation 6, x is proportional to $\mathrm{\xi }$, therefore $\frac{\partial \mathrm{x}}{\partial \mathrm{\xi }}$ and $\frac{\partial \mathrm{x}}{\partial \mathrm{\eta }}$ can be written as the following relations.

$\begin{array}{l}\frac{\partial x}{\partial \xi }={\displaystyle \sum }_{i=1}^8{x}_i\times N_i^{\eta }\\ \frac{\partial x}{\partial \eta }={\displaystyle \sum }_{i=1}^8{x}_i\times N_i^{\xi }\times \frac{\partial N_i^{\eta }}{\partial \eta }\end{array}$ (9)

Same as x, the y-coordinate is estimated by eight nodes laid down on the first and last strip edges. Therefore, y can be calculated by equation 10. In this equation, i is the node number, ${\mathrm{N}}_{\mathrm{i}}$ is mapping function of the i^th node and ${\mathrm{y}}_{\mathrm{i}}$ is the y-coordinate of the i^th node.

$y=\sum _{i=1}^8{N}_i{y}_i$ (10)

Based on the advantages of the strip elements mentioned in this section of the paper, y is simplified and calculated based on the node numbers 1, 2, 3 and 4 (four corner nodes of the element). Therefore, equation 10 can be rewritten as the following equation.

$y=\sum _{i=1}^8{N}_i{y}_i=N_1\times y_1+N_2\times y_2+N_3\times y_3+N_4\times y_4+\left(y_1+\frac{y_2-y_1}{3}\right)\times N_5+\left(y_1+\frac{2}{3}\left(y_2-y_1\right)\right)\times N_6+\left(y_4+\frac{2}{3}\left(y_3-y_4\right)\right)\times N_7+\left(y_4+\frac{y_3-y_4}{3}\right)\times N_8$ (11)

where, y-coordinate of nodes 5, 6, 7 and 8 is calculated based on the y-coordinate of node number 1, 2, 3 and 4. Inserting equations 7 and 8 into equation 11, it can be simplified to get the following relation.

$y=\sum _{i=1}^4{y}_i\times N_i^{\mathrm{*}}=y_1\times \frac{1}{2m}\times \left(\xi -m\right)\times \left(\eta -1\right)-y_2\times \frac{1}{2m}\times \left(\xi -m\right)\times \left(\eta +1\right)+y_3\times \frac{1}{2m}\times \xi \times \left(\eta +1\right)-y_4\times \frac{1}{2m}\times \xi \times \left(\eta -1\right)$ (12)

where, y-coordinate of strip element is estimated using four corner nodes of the element. In equation 12, ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{*}}$ is the modified mapping function of i^th node. Moreover, as the width of the strip element is constant in the longitudinal direction (${\mathrm{y}}_1={\mathrm{y}}_4,{\mathrm{y}}_3={\mathrm{y}}_2$), therefore equation 12 can be more simplified and rewritten as equation 13.

$y=\sum _{i=1}^2{N}_i^{\mathrm{*}\mathrm{*}}\times y_i=N_1^{\mathrm{*}\mathrm{*}}\times y_1+N_2^{\mathrm{*}\mathrm{*}}\times y_2=\left(\frac{1-\eta }{2}\right)y_1+\left(\frac{1+\eta }{2}\right)y_2$ (13)

In equation 13, ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{*}\mathrm{*}}$ is the simplified mapping function of strip element in y direction. As it is illustrated in this equation, y-coordinate of the element is estimated just with nodes number 1 and 2, and consequently y is only a function of $\mathrm{\eta }$. In other word, the simplified mapping function (${\mathrm{N}}_{\mathrm{i}}^{\mathrm{*}\mathrm{*}}$) in the y direction is only a function of $\mathrm{\eta }$. Thus, $\frac{\partial \mathrm{y}}{\partial \mathrm{\xi }}$ and $\frac{\partial \mathrm{y}}{\partial \mathrm{\eta }}$ of the strip element can be calculated by the following equations (see equation 14).

$\begin{array}{l}\frac{\partial y}{\partial \xi }=0\\ \frac{\partial y}{\partial \eta }=\frac{b}{2}\end{array}$ (14)

Therefore, the computational cost of calculation of Jacobian matrix is decreased dramatically due to the simplicity in evaluation of $\frac{\partial \mathrm{x}}{\partial \mathrm{\xi }}$, $\frac{\partial \mathrm{x}}{\partial \mathrm{\eta }}$, $\frac{\partial \mathrm{y}}{\partial \mathrm{\xi }}$ and $\frac{\partial \mathrm{y}}{\partial \mathrm{\eta }}$.

2.3 Jacobian matrix

The determinant of the Jacobian matrix ($\left|\mathrm{J}\right|$) is used to transform the integration from actual element to the parent element. It is calculated by equation 15 (Khayat et al. 2016a).

$\left|J\right|=\frac{\partial x}{\partial \xi }\times \frac{\partial y}{\partial \eta }-\frac{\partial x}{\partial \eta }\times \frac{\partial y}{\partial \xi }$ (15)

As it is shown in equation 15, the determinant of Jacobian matrix is related to the $\frac{\partial \mathrm{x}}{\partial \mathrm{\xi }}$, $\frac{\partial \mathrm{x}}{\partial \mathrm{\eta }}$, $\frac{\partial \mathrm{y}}{\partial \mathrm{\xi }}$ and $\frac{\partial \mathrm{y}}{\partial \mathrm{\eta }}$ terms which calculated in the previous section. Inserting equations 9 and 14 into equation 15, the determinant of the Jacobian matrix can be reached as equation 16.

$\left|J\right|=f\left(\eta \right)=\frac{b}{2}\sum _{i=1}^8{x}_i\times N_i^{\eta }$ (16)

where, $\mathrm{b}$ is the width of the strip element, i is node number, ${\mathrm{x}}_{\mathrm{i}}$ is the x-coordinate of the i^th node and ${\mathrm{N}}_{\mathrm{i}}^{\mathrm{\eta }}$ is the mapping function of ith node (see equation 9). As it is shown in equation 16, the determinant of Jacobian matrix is only a function of $\mathrm{\eta }$. Therefore, $\mathrm{\xi }$ and $\mathrm{\eta }$ coordinates can be separated in the integration process to calculate the load vector, linear stiffness matrix and geometrical stiffness matrix of strip element and result in decreasing the computational cost dramatically (Khayat et al. 2016a). Moreover, the partial derivatives with respect to x and y coordinates are presented in equation 17.

$\left\{\begin{array}{c}\frac{\partial }{\partial \mathrm{x}}\\ \frac{\partial }{\partial \mathrm{y}}\end{array}\right\}=\frac{1}{\left|J\right|}\times \left[\begin{array}{cc}\frac{\partial y}{\partial \eta }& -\frac{\partial y}{\partial \xi }\\ -\frac{\partial x}{\partial \eta }& \frac{\partial x}{\partial \xi }\end{array}\right]\left\{\begin{array}{c}\frac{\partial }{\partial \xi }\\ \frac{\partial }{\partial \eta }\end{array}\right\}$ (17)

For iso-parametric strip element used in the literature, approximation of geometry (x, y) follows the same rule as for displacement shape functions (Eccher et al. 2008). So, the geometry of iso-parametric strip element is approximated by following equations (Eccher et al. 2008).

$\begin{array}{l}x={\displaystyle \sum }_{i=-1}^{m+1}{\displaystyle \sum }_{j=1}^4{L}_j\left(\eta \right){\varphi }_i\left(\xi \right){\alpha }_{ij}\\ y={\displaystyle \sum }_{i=-1}^{m+1}{\displaystyle \sum }_{j=1}^4{L}_j\left(\eta \right){\varphi }_i\left(\xi \right){\beta }_{ij}\end{array}$ (18)

in which coefficients ${\alpha }_{ij}$ and ${\beta }_{ij}$ are calculated based on the geometry of each strip element. A full Jacobian matrix is derived using equation 18. Therefore, the Jacobian matrix determinant and derivatives of Jacobian matrix with respect to x and y are become dependent to $\xi$ and $\eta$. It leads to increase the computational cost of iso-parametric strip element (available element developed in the literature) compared to that of the sub-parametric element (developed in this study).

2.4 Strain- displacement relationships

In this paper, the strain- displacement relationship is assumed to be Sanders-type (Moradi et al. 2011, Khayat et al. 2016a, Khayat et al. 2016b). As it is illustrated in equation 19, the cylindrical panel strain vector $\left\{\mathrm{\epsilon }\right\}$ is comprised of linear part (${\left\{\mathrm{\epsilon }\right\}}^{\mathrm{l}}$) and nonlinear part (${\left\{\mathrm{\epsilon }\right\}}^{\mathrm{n}\mathrm{l}}$).

$\left\{\epsilon \right\}={\left\{\epsilon \right\}}^l+{\left\{\epsilon \right\}}^{nl}$ (19)

${\left\{\epsilon \right\}}^l=\left\{\begin{array}{c}\begin{array}{c}{\epsilon }_{xx}^l\\ {\epsilon }_{yy}^l\\ {\gamma }_{xy}^l\end{array}\\ k_{xx}^l\\ \begin{array}{c}{k}_{yy}^l\\ k_{xy}^l\\ \begin{array}{c}{\gamma }_{xz}^l\\ {\gamma }_{yz}^l\end{array}\end{array}\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial \mathrm{v}}{\partial \mathrm{y}}+\frac{w}{R}\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\end{array}\\ \frac{\partial {\psi }_x}{\partial x}\\ \begin{array}{c}\frac{\partial {\psi }_y}{\partial y}\\ \frac{\partial {\psi }_x}{\partial y}+\frac{\partial {\psi }_y}{\partial x}\\ \begin{array}{c}\frac{\partial w}{\partial x}+{\psi }_x\\ \frac{\partial w}{\partial y}+{\psi }_y-\frac{v}{R}\end{array}\end{array}\end{array}\right\}$ (20)

The components of linear part of strain vector (${\left\{\mathrm{\epsilon }\right\}}^{\mathrm{l}}$) are presented in equation 20 where, in this equation, ${\mathrm{\epsilon }}_{\mathrm{x}\mathrm{x}}^{\mathrm{l}}$ , ${\mathrm{\epsilon }}_{\mathrm{y}\mathrm{y}}^{\mathrm{l}}$, ${\mathrm{\gamma }}_{\mathrm{x}\mathrm{y}}^{\mathrm{l}}$, ${\mathrm{k}}_{\mathrm{x}\mathrm{x}}^{\mathrm{l}}$, ${\mathrm{k}}_{\mathrm{y}\mathrm{y}}^{\mathrm{l}}$, ${\mathrm{k}}_{\mathrm{x}\mathrm{y}}^{\mathrm{l}}$, ${\mathrm{\gamma }}_{\mathrm{z}\mathrm{x}}^{\mathrm{l}}$ and ${\mathrm{\gamma }}_{\mathrm{z}\mathrm{y}}^{\mathrm{l}}$ are the linear in-plane strain in x direction, in-plane linear strain in y direction, linear shear strain in xy plane, linear curvature strain about the x axis, linear curvature strain about the y axis, linear curvature strain in xy plane, linear shear strain in zx plane and linear shear strain in zy plane respectively. The components of nonlinear part of strain vector (${\left\{\mathrm{\epsilon }\right\}}^{\mathrm{n}\mathrm{l}}$) are presented in equation 21 (Wang et al. 1998), where ${\mathrm{\epsilon }}_{\mathrm{x}\mathrm{x}}^{\mathrm{n}\mathrm{l}}$, ${\mathrm{\epsilon }}_{\mathrm{y}\mathrm{y}}^{\mathrm{n}\mathrm{l}}$ and ${\mathrm{\gamma }}_{\mathrm{x}\mathrm{y}}^{\mathrm{n}\mathrm{l}}$ are the nonlinear in-plane strain in x direction, nonlinear in-plane strain in y direction and nonlinear shear strain in xy plane, respectively.

${\left\{\epsilon \right\}}^{nl}=\left\{\begin{array}{c}\begin{array}{c}{\epsilon }_{xx}^{nl}\\ {\epsilon }_{yy}^{nl}\\ {\gamma }_{xy}^{nl}\end{array}\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{c}\frac{1}{2}\left[{\left(\frac{\partial w}{\partial x}\right)}^2\right]\\ \frac{1}{2}\left[{\left(\frac{\partial w}{\partial y}-\frac{v}{R}\right)}^2\right]\\ \frac{\partial w}{\partial x}\left(\frac{\partial w}{\partial y}-\frac{v}{R}\right)\end{array}\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right\}$ (21)

The relationship between generalized linear stress vector ($\left\{{\mathrm{\sigma }}^{\mathrm{l}}\right\}$) and linear strain vector (${\left\{\mathrm{\epsilon }\right\}}^{\mathrm{l}}$) is presented in the following equation.

$\left\{\sigma \right\}=\left[M\right]{\left\{\epsilon \right\}}^l$ (22)

where the components of generalized linear stress vector (i.e. $\left\{{\mathrm{\sigma }}^{\mathrm{l}}\right\}$) is shown by the equation 23.

$\left\{\sigma \right\}={\left\{\begin{array}{ccc}{N}_{xx}^{}& N_{yy}^{}& \begin{array}{ccc}{N}_{xy}^{}& M_{xx}^{}& \begin{array}{ccc}{M}_{yy}^{}& M_{xy}^{}& \begin{array}{cc}{Q}_{xz}^{}& Q_{yz}^{}\end{array}\end{array}\end{array}\end{array}\right\}}^T$ (23)

In this equation, ${\mathrm{N}}_{\mathrm{x}\mathrm{x}}^{}$, ${\mathrm{N}}_{\mathrm{y}\mathrm{y}}^{}$, ${\mathrm{N}}_{\mathrm{x}\mathrm{y}}^{}$, ${\mathrm{M}}_{\mathrm{x}\mathrm{x}}^{}$, ${\mathrm{M}}_{\mathrm{y}\mathrm{y}}^{}$, ${\mathrm{M}}_{\mathrm{x}\mathrm{y}}^{}$, ${\mathrm{Q}}_{\mathrm{x}\mathrm{z}}^{}$ and ${\mathrm{Q}}_{\mathrm{y}\mathrm{z}}^{}$ are the membrane stress in x direction, membrane stress in y direction, shear membrane stress in xy plane, moment stress about y axis, moment stress about x axis, torsional stress in xy plane, transverse shear stress in xz plane and transverse shear stress in yz plane, respectively (Moradi et al. 2011).

$M=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{A}_{11}& A_{12}\\ A_{12}& A_{22}\end{array}& \begin{array}{cc}{A}_{16}& B_{11}\\ A_{26}& B_{12}\end{array}\\ \begin{array}{cc}{A}_{16}& A_{26}\\ B_{11}& B_{12}\end{array}& \begin{array}{cc}{A}_{66}& B_{16}\\ B_{16}& D_{11}\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{B}_{12}& B_{16}\\ B_{22}& B_{26}\end{array}& \begin{array}{cc}{0}_{}& 0_{}\\ 0_{}& 0_{}\end{array}\\ \begin{array}{cc}{B}_{26}& B_{66}\\ D_{12}& D_{16}\end{array}& \begin{array}{cc}{0}_{}& 0_{}\\ 0_{}& 0_{}\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{cc}{B}_{12}& B_{22}\\ B_{16}& B_{26}\end{array}& \begin{array}{cc}{B}_{26}& D_{12}\\ B_{66}& D_{16}\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{D}_{22}& D_{26}\\ D_{26}& D_{66}\end{array}& \begin{array}{cc}{0}_{}& 0_{}\\ 0_{}& 0_{}\end{array}\\ \begin{array}{cc}{0}_{}& 0_{}\\ 0_{}& 0_{}\end{array}& \begin{array}{cc}{A}_{44}& A_{45}\\ A_{45}& A_{55}\end{array}\end{array}\end{array}\right]$ (24)

In equation 22, $\left[\mathrm{M}\right]$ is the material constitutive matrix which is expanded by equation 24. where, ${\mathrm{A}}_{\mathrm{i}\mathrm{j}}$, ${\mathrm{B}}_{\mathrm{i}\mathrm{j}}$ and ${\mathrm{D}}_{\mathrm{i}\mathrm{j}}$ are the extensional, bending, and bending-extensional coupling rigidities respectively and can be calculated from the following relations (Moradi et al. 2011).

$\left(A_{ij},B_{ij},D_{ij}\right)=\sum _{k=1}^N\underset{t_{k-1}}{\overset{t_k}{\int }}{\left({\stackrel{-}{Q}}_{ij}\right)}_k\left(1,Z,Z^2\right)dZ(i,j=\mathrm{1,2},6)$ (25)

Note that, ${\mathrm{A}}_{44}$, ${\mathrm{A}}_{55}$ and ${\mathrm{A}}_{45}$ are the shear stiffness constants (Moradi et al. 2011). In equation 25, N is the number of layers, and ${\mathrm{t}}_{\mathrm{k}-1}$ and ${\mathrm{t}}_{\mathrm{k}}$ are the depth of lower and upper surfaces of the k^th layer, respectively. Moreover, ${\left({\stackrel{-}{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}\right)}_{\mathrm{k}}$ are the transformed elastic constants which are obtained from equation 26 (Eccher et al. 2008, Moradi et al. 2011, Khayat et al. 2016a).

$\left[\begin{array}{c}\begin{array}{c}{\stackrel{-}{Q}}_{11}\\ {\stackrel{-}{Q}}_{12}\end{array}\\ \begin{array}{c}{\stackrel{-}{Q}}_{22}\\ {\stackrel{-}{Q}}_{16}\end{array}\\ \begin{array}{c}{\stackrel{-}{Q}}_{26}\\ {\stackrel{-}{Q}}_{66}\end{array}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}{c}^4\\ c^2{s}^2\end{array}& \begin{array}{c}2c^2{s}^2\\ c^4+s^4\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{s}^4\\ c^2{s}^2\end{array}& \begin{array}{c}4c^2{s}^2\\ -4c^2{s}^2\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{c}{s}^4\\ c^{3}s\end{array}& \begin{array}{c}{2c}^2{s}^2\\ c{s}^3-c^{3}s\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}{c}^4\\ -c{s}^3\end{array}& \begin{array}{c}4c^2{s}^2\\ -2cs(c^2-s^2)\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{c}{c}^3{s}^3\\ c^2{s}^2\end{array}& \begin{array}{c}{c}^{3}s-c{s}^3\\ -2c^2{s}^2\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}-c^{3}s\\ c^2{s}^2\end{array}& \begin{array}{c}2cs(c^2-s^2)\\ {(c^2-s^2)}^2\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{Q}_{11}\\ Q_{12}\end{array}\\ \begin{array}{c}{Q}_{22}\\ Q_{66}\end{array}\end{array}\right]$ (26)

where, c and s stand for the $\cos \mathrm{\theta }$ and $\sin \mathrm{\theta }$, in which $\mathrm{\theta }$ is the lamination angle for each layer. ${\mathrm{Q}}_{\mathrm{i}\mathrm{j}}$ are rigidities of an orthotropic layer in its material directions (1 and 2) and are calculated by the following equation.

$\begin{array}{l}{Q}_{11}=\frac{E_{11}}{\left(1-v_{12}{v}_{21}\right)}\\ Q_{22}=\frac{E_{22}}{\left(1-v_{12}{v}_{21}\right)}\\ Q_{12}=\frac{v_{21}{E}_{11}}{\left(1-v_{12}{v}_{21}\right)}\\ Q_{66}=G_{12}\end{array}$ (27)

In equation 27, ${\mathrm{E}}_{11}$, ${\mathrm{E}}_{22}$, ${\mathrm{G}}_{12}$, ${\mathrm{v}}_{12}$ and ${\mathrm{v}}_{21}$ are material properties which are well-known in the literature (Eccher et al. 2008, Moradi et al. 2011, Khayat et al. 2016a).

2.5 Strip element linear stiffness matrix

The linear stiffness matrix components are calculated based on the internal virtual work (${\mathrm{\delta }\mathrm{W}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\mathrm{L}}$) done by components of the stress vector ($\left\{\mathrm{\sigma }\right\}$) on the strip element volume. The internal virtual work (${\mathrm{\delta }\mathrm{W}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\mathrm{L}}$) is written as equation 28.

${\delta W}_{int}^L=\underset{-1}{\overset{+1}{\int }}\underset{0}{\overset{m_i}{\int }}\left[\left(N_{xx}^{}\delta {\epsilon }_{xx}^l\right)+\left(N_{yy}^{}\delta {\epsilon }_{yy}^l\right)+\left(N_{xy}^{}\delta {\gamma }_{xy}^l\right)+\left(M_{xx}^{}\delta k_{xx}^l\right)+\left(M_{yy}^{}\delta k_{yy}^l\right)+\left(M_{xy}^{}\delta k_{xy}^l\right)+\left(Q_{xz}^{}\delta {\gamma }_{xz}^l\right)+\left(Q_{yz}^{}\delta {\gamma }_{yz}^l\right)\right]\left|J\right|d\xi d\eta$ (28)

where ${\mathrm{m}}_{\mathrm{i}}$ is the number of ith element knots on a nodal line. The integrations on $\mathrm{\xi }$ and $\mathrm{\eta }$ are separated from each other because of the simplicity of Jacobian matrix determinant (see equation 16). So, the one-dimensional Gauss quadrature integration method is applied separately in each direction to compute the equation 28 (Eccher et al. 2008). In the iso-parametric element presented in the literature (Eccher et al. 2008), the determinant of Jacobian is complicated and the integrations on $\mathrm{\xi }$ and $\mathrm{\eta }$ are dependent to each other. This prevents of decomposition of surface integrals encountered in the calculation of stiffness matrix. It proves that the computational cost of sub-parametric strip element developed in this study is lower than that of the iso-parametric strip element due to the simplicity of its Jacobian matrix. According to the internal virtual work, the linear stiffness matrix of the strip element ($\left[{\mathrm{K}}_{\mathrm{E}}\right]$) is obtained from equation 29.

${\delta W}_{int}^L={\left\{\mathrm{\delta }\mathrm{\Delta }\right\}}^T\left[K_E\right]\left\{\mathrm{\Delta }\right\}$ (29)

in which, $\left\{\mathrm{\Delta }\right\}$ is the strip element displacement vector (see equations 4 and 5) and ${\left\{\mathrm{\delta }∆\right\}}^{\mathrm{T}}$is transpose of the virtual displacement vector of the panel.

2.6 Strip element force vector

There are two types of strip element based on the panel loading condition and location of the cutout on the panel. In Figure 7, a circular cutout with free edges is considered in the center of the panel. The figure shows two types of strip elements, extended strip element (i.e. i^th element) and truncated strip element (i.e. j^th element). The length of the extended strip element is equal to the panel length. The end edges of extended strip element are under compressive distributed load with magnitude of q. The truncated strip element is between the panel transverse edge and cutout edge; therefore, its length is less than the length of the panel. Also as can be seen from Figure 7, for the truncated strip elements only one transverse edge is under compression load of q in u direction. The external virtual work (${\mathrm{\delta }\mathrm{W}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{q}}$) done by compressive load applied on transverse edges of the cylindrical panel is written to obtain the load vector of the panel by the following relation (Khayat et al. 2016a):

${\delta W}_{ext}^q=\sum _{i=1}^{N_e}\underset{0}{\overset{b_i}{\int }}\left(q_{i(x=0)}\times \delta u_{i(x=0)}-q_{i(x=l_i)}\times \delta u_{i(x=l_i)}\right)dy$ (30)

In equation 30, ${\mathrm{N}}_{\mathrm{e}}$ is the total number of strip elements, ${\mathrm{b}}_{\mathrm{i}}$ is the width of ith strip element,${\mathrm{l}}_{\mathrm{i}}$ is the length of ith strip element, ${\mathrm{q}}_{\mathrm{i}}$ is the distributed compressive load on the element edges in direction u and $\mathrm{\delta }{\mathrm{u}}_{\mathrm{i}}$ is the virtual displacement. For the truncated strip elements, one of the ${\mathrm{q}}_{\mathrm{i}(\mathrm{x}={\mathrm{l}}_{\mathrm{i}})}$ or ${\mathrm{q}}_{\mathrm{i}(\mathrm{x}=0)}$ is set to zero. Also due to the costraint on u in $x=l_i$, equation 30 can be summerized as equation 31.

${\delta W}_{ext}^q=\sum _{i=1}^{N_e}\underset{0}{\overset{b_i}{\int }}{q}_{i(x=0)}\times \delta u_{i(x=0)}dy$ (31)

According to Figure 1, y-coordinate is along circumferential direction of the strip and can be related to the natural coordinate $\mathrm{\eta }$ as:

$\mathrm{d}\mathrm{y}=\frac{b_i}{2}d\eta$ (32)

Virtual displacement $\mathrm{\delta }{\mathrm{u}}_{\mathrm{i}(\mathrm{x}=0)}$ can be rewritten as equation 33 based on equation 3.

$\delta u_{i(x=0)}=\sum _{s=1}^4\sum _{r=-1}^1{L}_s{\varphi }_r\times \delta u_{sr}$ (33)

It should be noted that, three spline functions (${\varphi }_{-1}$, ${\varphi }_0$ and ${\varphi }_1$) are contributed to the external virtual work (${\delta W}_{ext}^q$) (Eccher et al. 2008). Substituting equations 32 and 33 into equation 31 leads to following equation.

${\delta W}_{ext}^q=\sum _{i=1}^{N_e}\frac{b_i}{2}\sum _{s=1}^4\left({\varphi }_{-1}\delta u_{-1s}+{\varphi }_0\delta u_{0s}+{\varphi }_1\delta u_{1s}\right)\underset{-1}{\overset{+1}{\int }}{L}_s\times q_i\times d\eta$ (34)

Finally, equation 30 is simplified as equation 35 to obtain the load vector of cylindrical panel.

${\delta W}_{ext}^q={\left\{\delta \mathrm{\Delta }\right\}}^T\left\{F\right\}$ (35)

where, $\left\{\mathrm{F}\right\}$ and ${\left\{\mathrm{\delta }\mathrm{\Delta }\right\}}^{\mathrm{T}}$ are the load vector and transpose of the virtual displacement of the cylindrical panel, respectively.

Figure 7
Cylindrical panel loading.

2.7 Strip geometrical stiffness matrix

An initial deflection is induced in the cylindrical panel due to the distributed compressive load applied on the edges of the panel. This pre-buckling displacement vector is calculated by the following equation.

$\left\{{∆}_0\right\}=\left\{F\right\}{\left[K_E\right]}^{-1}$ (36)

where, $\left\{{∆}_0\right\}$ is the pre-buckling displacement vector, $\left\{\mathrm{F}\right\}$ is the panel load vector and $\left[{\mathrm{K}}_{\mathrm{E}}\right]$ is the linear stiffness matrix of the cylindrical panel. Based on the initial displacement vector of the panel, the pre-buckling generalized linear stress vector $\left\{{\mathrm{\sigma }}_0\right\}$ is obtained by equation 22. In the linearized buckling analysis it is assumed that the load-deflection relation is linear up to the onset of buckling perturbation. The stresses induced in this stage will remain constant and are independent of buckling deformations. In the geometrical nonlinear analysis which is another way to obtain buckling load of structures, the stress vector $\left\{{\mathrm{\sigma }}_0\right\}$ depend on the linear and nonlinear strains. This makes the problem to be nonlinear in terms of displacements. This pre-buckling bending deformations have significant effects on the calculation of the critical axial load in shells. The components of the initial generalized linear stress vector $\left\{{\mathrm{\sigma }}_0\right\}$ are shown in equation 37.

$\left\{{\sigma }_0\right\}={\left\{\begin{array}{ccc}{N}_{xx}^0& N_{yy}^0& \begin{array}{ccc}{N}_{xy}^0& M_{xx}^0& \begin{array}{ccc}{M}_{yy}^0& M_{xy}^0& \begin{array}{cc}{Q}_{xz}^0& Q_{yz}^0\end{array}\end{array}\end{array}\end{array}\right\}}^T$ (37)

From the internal virtual work done by the initial linear in-plane stresses (${\mathrm{\delta }\mathrm{W}}_{\mathrm{i}\mathrm{n}\mathrm{t}}^{\mathrm{n}\mathrm{l}}$) (Khayat et al. 2016a), the geometrical stiffness matrix of the panel is obtained (see equation 38).

${\delta W}_{int}^{nl}=\sum _{i=1}^{N_e}\underset{-\frac{b_i}{2}}{\overset{\frac{b_i}{2}}{\int }}\underset{0}{\overset{l_i}{\int }}\left[\left(N_{xx}^0\delta {\epsilon }_{xx}^{nl}\right)+\left(N_{yy}^0\delta {\epsilon }_{yy}^{nl}\right)+\left(N_{xy}^0\delta {\gamma }_{xy}^{nl}\right)\right]dxdy$ (38)

In equation 38, ${\mathrm{N}}_{\mathrm{x}\mathrm{x}}^0$ is the initial membrane stress in x direction, ${\mathrm{N}}_{\mathrm{y}\mathrm{y}}^0$ is the initial membrane stress in y direction, ${\mathrm{N}}_{\mathrm{x}\mathrm{y}}^0$ is the initial shear membrane stress in xy plane, ${\mathrm{\epsilon }}_{\mathrm{x}\mathrm{x}}^{\mathrm{n}\mathrm{l}}$ is the nonlinear in-plane strain in x direction, ${\mathrm{\epsilon }}_{\mathrm{y}\mathrm{y}}^{\mathrm{n}\mathrm{l}}$ is the nonlinear in-plane strain in y direction and ${\mathrm{\gamma }}_{\mathrm{x}\mathrm{y}}^{\mathrm{n}\mathrm{l}}$ is the nonlinear shear strain in xy plane. equation 38 is rewritten and simplified as equation 39.

${\delta W}_{int}^{nl}={\left\{\mathrm{\delta }\mathrm{\Delta }\right\}}^T\left[K_G\right]\left\{\mathrm{\Delta }\right\}$ (39)

where, $\left[{\mathrm{K}}_{\mathrm{G}}\right]$ is the geometrical stiffness matrix, $\left\{\mathrm{\Delta }\right\}$ is the displacement vector and $\left\{\delta ∆\right\}$ is the virtual displacement vector of the strip element.

2.8 Buckling analysis

The linear stiffness matrix ($\left[{\mathrm{K}}_{\mathrm{E}}\right]$) and geometrical stiffness matrix ($\left[{\mathrm{K}}_{\mathrm{G}}\right]$) of the cylindrical panel are obtained by assembling linear and geometrical stiffness matrices of all strip elements. To calculate the buckling load of the panel, the eigenvalue problem presented in equation 40 is solved by the inverse iteration method (Khayat et al. 2016b). In this equation, $\mathrm{\lambda }$ is the eigenvalue and $\mathrm{\varnothing }$ is its associated eigenvector.

$(K_E-\lambda K_G)\left\{\mathrm{\varnothing }\right\}=0$ (40)

3 Verification

Based on the formulation presented in the previous sections, a computer program namely SSFSM, was written in FORTRAN to compute the buckling load of a cylindrical panel having cutouts with different shapes. In order to check the validity of the presented formulation, several case studies are carried out and the results obtained from SSFSM are compared with the results available in the literature and the results obtained from the ABAQUS commercial finite element software (ABAQUS 2006). The details are as follows.

3.1 Buckling load of cylindrical panel without cutout

Au and Cheung (1996) calculated the buckling load of a square cylindrical panel without any cutout using spline finite strip method. In their research, the length of the panel, panel thickness, modulus of elasticity and poison ratio were considered as 600 mm, 5 mm, $70\mathrm{k}\mathrm{N}/{\mathrm{m}\mathrm{m}}^2$ and 0.3, respectively. Moreover, all edges of the panel were simply supported. A view of the cylindrical panel considered in the Au and Cheung study is shown in Figure 8. In their research, a uniform compressive load was applied at the transverse edges of the panel. The cylindrical panel was simulated in ABAQUS using 3600 S8R5 elements and roughly 11000 degrees of freedom (ABAQUS 2006). S8R5 is an eight nodded iso-parametric shell element. In this element, the stiffness matrix is obtained using reduced integraton. A cylindrical coordinate system is defined in the ABAQUS to define/apply the boundary conditions on the panel edges. A unique uniform load is applied on the curved edge of the panel using the option named as "shell edges load" in the ABAUQS. "Buckle analysis" is performed to calculate the panel buckling load in the ABAQUS. Their panel was modeled in SSFSM with 6 strip elements, 20 sections and roughly 1900 degrees of freedom. A view of the cylindrical panel modeled in the SSFSM and ABAQUS is shown in Figure 8. The cylindrical panel buckling load for different central angle (0, 0.1, 0.2 and 0.3 radian) is calculated by SSFSM and ABAQUS. In Figure 9, the buckling load obtained from the SSFSM is compared with those of Au and Cheung (1996) and ABAQUS. As it is illustrated in the figure, as the central angle of the panel is increased, the buckling load is increased. The figure also shows that there is a good agreement between results obtained from the SSFSM and those of the Au and Cheung (1996) and ABAQUS.

Moreover, in order to verify the convergence ability of the SSFSM in calculating the panel buckling load, the buckling load of the panel for different number of strip elements and knots are calculated and presented in Table 1. As it is indicated in the table, 6 strip elements and 15 knots are sufficient for convergence of the panel buckling load.

Figure 8
SSFSM and ABAQUS models of cylindrical panel (a) ABAQUS (b) SSFSM.

Figure 9
Cylindrical panel buckling load obtained from Au and Cheung, SSFSM and ABAQUS

Table 1
Convergence of SSFSM buckling loads

Number of strips		2	3	4	5	6	7	8
Number of knots	5	183.2950	183.0030	182.7109	182.7028	182.6947	182.7677	182.8407
	10	183.2098	182.9117	182.6136	182.5649	182.5162	182.5933	182.6704
	15	183.1247	182.8204	182.5162	182.5081	182.5000	182.5000	182.5000
	20	183.0395	182.7779	182.5162	182.5081	182.5000	182.5000	182.5000
	25	182.9543	182.7353	182.5162	182.5081	182.5000	182.5000	182.5000

In order to evaluate the computational cost of SSFSM, a comparsion is performed between execution time of SSFSM (sub-parametric element) and the same code with iso-parametric element. In this regrad, the panel described in this section (Figure 8) is discretized by six sub-parametric strip elements and the buckling load of the panel is calculated. Different numbers of konts are considered to simulate the panel. The execution time for both codes is measured for different degrees of freedom. The solution time of code is dependent on the number of degrees of freedom. For instance, when degrees of freedom is 2040, the solution time of SSFSM is about 148 seconds for a 2.8 GHz CPU. Figure 10 shows the computational cost of SSFSM against the iso-parametric code. In this figure, the precentage of reduction in execution time for sub-parametric element is drawn against the number of degrees of freedom. As shown in Figure 10, the execution time of SSFSM is significantly lower than that of the iso-parametric one. Two factors are responsible for drastical reduction of SSFSM computational cost. The simplicity of Jacobian determinant in the SSFSM causes decomposition of surface integral encountered in the linear elastic and geometric stiffness matrices. The number of Gauss points requried in the SSFSM to accuratly integrate the stiffness matrices is lower than those of the iso-parametric. As the degrees of freedom is increased from 675 to 2040, the precentage of reduction in the execution time is reduced from 62 to 57.5 percent. This reduction is related to the calculation of eigenvalue problem and as the degrees of freedom increases, the time to solve the eigenvalue problem increases.

Figure 10
Computational cost of sub-parametric strip element developed in this study.

3.2 Buckling load of cylindrical panel with cutout

A composite cylindrical panel with layup configurations of ${\left(\mathrm{\theta }/-\mathrm{\theta }\right)}_{2\mathrm{s}}$ and ${\left(\mathrm{\theta }/-\mathrm{\theta }\right)}_4$ is simulated in ABAQUS and SSFSM. The angle layer $\mathrm{\theta }$ is considered as 0, 15, 30, 45, 60, 75 and 90 degrees. The material properties of panel are presented in Table 2. The panel length, width and radius are considered 180 mm, 180 mm and 1800 mm, respectively. The ratio of thickness to length of the panel is 1/45. The panel comprises of eight layers and all edges are simply supported. The panel is under uniform compressive load on the transverse edges. A cutout with a square shape and area of 1296 ${\mathrm{m}\mathrm{m}}^2$ was considered in the center of the panel.

Figure 11
View of cylindrical panel with square cutout (a) ABAQUS (b) SSFSM.

Figure 12
Panel buckling load obtained from ABAQUS and SSFSM for layup ${\left(\theta /-\theta \right)}_4$.

The cylindrical panel was modeled in the ABAQUS with 700 of S8R5 elements type. Moreover, the panel was simulated in the SSFSM with 10 strip elements and 20 knots. ABAQUS and SSFSM models are shown in Figure 11. The buckling loads of the panel calculated from the SSFSM and ABAQUS are shown in Figures 12 and 13. As seen from Figures 12 and 13, the maximum and minimum buckling load of the panel is about 490 and 221 N/mm for layup angle 45 and 90 degree, respectively. The maximum difference between results obtained from the SSFSM and those of the ABAQUS is about 4 percent, which shows that there is a good agreement between the result obtained from the SSFSM and the ABAQUS.

Figure 13
Panel buckling load obtained from ABAQUS and SSFSM for layup ${\left(\theta /-\theta \right)}_{2s}$.

4 Evaluation of SSFSM capability in calculation of panel axial buckling load

A parametric study was carried out to clarify the capability of the SSFSM in calculating the buckling load of the composite cylindrical panel with cutout. The effects of panel boundary condition, cutout shape, cutout area, panel thickness, central angle of the panel and layup angle configuration on buckling load of panel were investigated. The details are as follows.

4.1 Effect of cutout shape and panel thickness on panel buckling load

In this section, the effects of cutout shape and panel thickness are investigated on the buckling load of the panel. To this end, the effect of several cutout shapes such as square, diamond, circular and elliptical on the buckling load of the panel are investigated. All cutouts are considered in the center of the panel and their cutout area are considered to be 1296 ${\mathrm{m}\mathrm{m}}^2$ (the ratio of cutout area to panel area is 4%). The ratio of big chord to small chord of the elliptical cutout is kept at 0.5, $2/3$, 1, 1.5 and 2. The panel is simply supported on its edges and the in-plane displacements are considered to be free on all of its edges. The edges of the cutout are considered to be free of any mechanical constraint. The boundary of cutout coincides with boundaries of small strips and is approximated by eight cubic polynomials, and subsequently, 24 nodes are contributed in approximating the boundary of circular/elliptical cutout. The panel width, length and central angle are 180 mm, 180 mm and 0.3 radians, respectively. The properties of the panel material are presented in Table 2.

Table 2
Mechanical properties of panel material

${\mathbf{E}}_{11}$	${\mathbf{E}}_{22}$	${\mathbf{G}}_{12}$	${\mathbf{G}}_{13}$	${\mathbf{G}}_{23}$	$\mathbf{v}$
$181\frac{\mathrm{k}\mathrm{N}}{{\mathrm{m}\mathrm{m}}^2}$	$10.3\frac{\mathrm{k}\mathrm{N}}{{\mathrm{m}\mathrm{m}}^2}$	$7.17\frac{\mathrm{k}\mathrm{N}}{{\mathrm{m}\mathrm{m}}^2}$	$7.17\frac{\mathrm{k}\mathrm{N}}{{\mathrm{m}\mathrm{m}}^2}$	$6.21\frac{\mathrm{k}\mathrm{N}}{{\mathrm{m}\mathrm{m}}^2}$	$0.28$

The parametric study is performed for two cross-ply layups configuration of ${(0/90)}_4$ and ${(0/90)}_{2\mathrm{s}}$. Three types of panel thickness to length ratio of 1/25, 1/50 and 1/100 are used in the parametric study. The cylindrical panel is modeled in the SSFSM using 12 four-nodded strip elements and 20 knots (roughly 3700 degrees of freedom). Figure 14 shows the configuration of the panel with different cutout shapes as modeled in SSFSM. The buckling loads of the panels computed from SSFSM are tabulated in Table 3.

The results presented in Table 3 indicate that the buckling load of the panel increases as its thickness is increased. In case of the panel without any cutout, when its thickness is getting larger by 2 and 4 times, the buckling load of the panel is increased by 3.7 and 26 times, respectively. It means that the buckling load of the panel is changed nonlinearly in terms of its thickness. In all cases, buckling load of the panel with cutout is less than that of the panel without any cutout. Although all cutouts areas are the same, the cutout shape has an important effect on the buckling load.

Figure 14
Cylindrical panel simulated by SSFSM (a) Without cutout (b) Square cutout (c) Diamond cutout (d) Elliptical cutout with ratio of 0.5 (e) Elliptical cutout with ratio of 2/3 (f) Circular cutout (g) Elliptical cutout with ratio of 1.5 (h) Elliptical cutout with ratio of 2.

When the elliptical cutout rotates about the normal axis to the center of the panel, the buckling load varies, dramatically. In all cases, when big chord of the elliptical cutout is normal to the compressive load direction (elliptical cutout with ratio 2), the cutout has negligible effect on the buckling load of the panel compared with other cutout shapes. For these kinds of cutouts, the panel buckling load decreases by 17%, 11% and 3% for thick, medium and thin panels respectively. It indicates that when the panel thickness is decreased, the effect of cutout on changes of the panel buckling load is decreased. In most cases, when the big chord of the elliptical cutout is parallel to compressive load direction (elliptical cutout with ratio 0.5), the maximum reduction is seen in the panel buckling load. The buckling load of the panel with diamond cutout is always larger than that of the panel with square cutout. Moreover, the buckling load of the panel with square cutout is always less than that of the panel with circular cutout. As it is illustrated in Table 3, the layup configuration has an interesting influence on the buckling load of the panel. In the most cases, the buckling load of the panel with layup configuration of ${[0/90]}_{2\mathrm{s}}$ is more than that of ${[0/90]}_4$ configuration. But as an interesting point, when the elliptical cutouts with ratio 1.5 and 2 are considered in the center of the panel, the buckling load for ${[0/90]}_4$ layup configuration is more than that of the ${[0/90]}_{2\mathrm{s}}$.

4.2 Effects of cutout size, layup, panel central angle and boundary conditions on the buckling load

In this section, the effects of cutout area, layup configuration and panel central angle are investigated on the buckling load of the panel. A circular cutout is considered at the center of the panel (see Figure 14(f)). The panel length, width and the ratio of panel thickness to its length are considered as 180 mm, 180 mm and 1/50, respectively. The panel is composed of eight layers and its material properties are presented in Table 2. The cutout area is considered as 5%, 10%, 15%, 20% and 25% of panel area. Furthermore, the central angle of the panel is considered as 0.1, 0.3 and 0.5 radian. The layup arrangements of ${\left[0\right]}_8$, ${\left[90\right]}_8$, ${[0/90]}_{2\mathrm{s}}$, ${[0/90]}_4$, ${\left[45\right]}_8$, ${[45/-45]}_4$, ${[30/-30]}_4$ and ${[0/45/90/-45]}_{\mathrm{s}}$ are used in the study. At the first step, the boundary conditions of the panel and cutout are the same as those used in previous section. The buckling loads of the simply supported panel calculated from the SSFSM are presented in Table 4.

Table 3
Panel buckling load for different cutout shapes.

Panel thickness Thick (1/25)	Cutout shape	Buckling load (N/mm)
		${(0/90)}_4$	${(0/90)}_{2\mathbf{s}}$
Thick (1/25)	Without cutout	2201.668	2244.298
	Square	1792.830	1822.665
	Diamond	1712.729	1730.092
	Elliptical (a/b=0.5)	1640.561	1733.972
	Elliptical (a/b=2/3)	1697.209	1762.393
	Circle	1772.675	1798.089
	Elliptical (a/b=1.5)	1847.753	1829.905
	Elliptical (a/b=2)	1912.064	1859.781
Medium (1/50)	Without cutout	363.776	368.892
	Square	313.892	317.070
	Diamond	304.338	305.599
	Elliptical (a/b=0.5)	286.228	297.150
	Elliptical (a/b=2/3)	297.237	304.590
	Circle	310.934	313.378
	Elliptical (a/b=1.5)	323.922	321.119
	Elliptical (a/b=2)	334.446	327.491
Thin (1/100)	Without cutout	84.535	86.0832
	Square	79.550	80.271
	Diamond	78.240	78.357
	Elliptical (a/b=0.5)	72.711	74.263
	Elliptical (a/b=2/3)	75.844	76.921
	Circle	79.298	79.802
	Elliptical (a/b=1.5)	82.169	82.062
	Elliptical (a/b=2)	84.198	83.566

The results presented in Table 4 indicate that the panel buckling load is decreased as the cutout area is increased or central angle of the panel is decreased. For instance, in case of ${\left[0\right]}_8$ layup and central angle of 0.1 radian, when the area contribution is changed from 5% to 25% (cutout area is increased 5 times), the buckling load of the panel is decreased from 191 to 110.7 N/m (the panel buckling load is decreased 42%). On the other hand, for the same layup and cutout area ratio of 25%, when the central angle is changed from 0.5 to 0.1 radian (central angle is decreased 5 times), the buckling load of the panel is decreased from 151 to 110.7 N/m (the panel buckling load is decreased 27%). This shows that the variation of cutout area has more effect on the buckling load than the variation of the panel central angle. The results obtained in this section indicate that the panel layup configuration is one of the most important factors on its buckling load. The panel buckling load is reached at the maximum magnitude of 707.5 N/m for cutout area ratio of 5%, central panel angle of 0.5 radian and ${[0/45/90/-45]}_{\mathrm{s}}$ layup configuration. While, it is reached at 639.6 and 513.3 N/m in ${[45/-45]}_4$ and ${[30/-30]}_4$ layup configuration, respectively, the minimum buckling load is reached at 110.6 N/m for ${\left[0\right]}_8$ layup configuration. As an interesting point, in case of ${\left[90\right]}_8$ layup and cutout area ratio of 25%, when the central angle is varied from 0.1 to 0.3, the panel buckling load is changed from 127.5 to 127.6 N/m which shows only a slight variation.

Figure 15 shows the buckling modes of laminated composite cylindrical panel having 15% central circular cutout with different layups and central angles. Panel stacking sequence is considered as [90]₈, [45]₈, [45/-45]₄ and [0/45/90/-45]_s. The central angle of the panel is 0.1 and 0.5 radians. For uni-directional layups [90]₈ and [45]₈ wavelength in the fiber direction is more than that of fiber transverse direction. This is pronounced with increasing in the panel central angle. The local shape of modes in these cases justifies reduction of buckling loads. The buckling mode shapes of multi-directional laminates such as [45/-45]₄ and [0/45/90/-45]_s behave globally and therefore endure more axial buckling loads.

To show the capability of SSFSM to consider the effects of different cylindrical boundary conditions on the buckling load, the buckling load of the panel was calculated using SSFSM when the loading edges of the panel were fixed. The panel properties and geometry were considered as mentioned at the beginning of this section. All layups presented in Table 4 were considered in the analysis. The central angle of the panel was considered as 0.1 and 0.5 radian. The area ratio of circular cutout was considered as 5%, 10%, 15%, 20% and 25%. The buckling load of the simply supported cylindrical panel (simple panel) is compared with that of fix supported cylindrical panel (fixed panel) in Figures 16 to 23.

Table 4
Buckling load of cylindrical panel with circular cutout.

Layup	Panel central angle (Radian)	buckling load (N/mm)
		Cutout area to panel area ratio
		$5\mathbf{\%}$	$10\mathbf{\%}$	15%	20%	25%
${\left[0\right]}_8$	0.1	191.565	157.045	136.396	121.755	110.602
	0.3	224.302	185.779	159.564	140.111	124.960
	0.5	291.893	240.619	202.801	173.774	150.949
${\left[90\right]}_8$	0.1	152.005	139.993	133.221	129.046	127.449
	0.3	200.351	164.023	144.246	133.015	127.633
	0.5	269.667	193.089	157.496	138.601	128.390
${[0/90]}_{2\mathbf{s}}$	0.1	211.942	184.828	169.256	158.427	149.793
	0.3	307.965	275.370	250.331	228.820	209.289
	0.5	493.606	441.490	390.148	343.196	301.987
${[0/90]}_4$	0.1	207.756	181.447	166.355	155.683	146.882
	0.3	305.682	273.351	248.185	226.400	206.613
	0.5	491.235	438.266	384.229	335.748	294.745
${\left[45\right]}_8$	0.1	184.912	171.255	164.816	157.359	148.421
	0.3	237.307	207.731	188.777	172.980	159.064
	0.5	334.709	263.502	219.285	190.414	170.618
${[45/-45]}_4$	0.1	386.753	370.087	368.392	366.912	361.698
	0.3	474.016	439.697	419.616	405.847	392.931
	0.5	639.636	563.598	503.161	465.010	439.207
${[30/-30]}_4$	0.1	339.907	317.775	307.730	298.518	287.405
	0.3	397.086	366.716	345.626	327.320	309.259
	0.5	513.295	460.208	416.627	380.318	348.782
${[0/45/90/-45]}_{\mathbf{s}}$	0.1	271.319	247.538	235.994	227.448	218.481
	0.3	423.066	378.692	342.118	310.317	281.378
	0.5	707.511	603.327	502.818	420.998	358.210

The buckling load of the panel with layup ${\left[0\right]}_8$ is shown in Figure 16. As shown in Figure 16, the buckling load of a fixed panel with layup ${\left[0\right]}_8$ is greater than that of simple panel. As the cutout area ratio is increased or central angle of the panel is decreased, the buckling load is decreased. The maximum buckling load of fixed panel is about 840 N/mm. In the fixed panel, the effect of central angle on the buckling load of panel decreases dramatically when cutout area ratio is increased. In the simple panel with central angle of 0.1 radians, the area of the cutout has a negligible effect on the buckling load. The buckling load of the panel with layup ${\left[90\right]}_8$ is presented in Figure 17. As an interesting point, the boundary condition of the panel has not considerable effect on the variation of buckling load of the panel. For instance, the maximum buckling load of the simple panel is 1% lower than that of the fixed panel. This is due to the fact that the lamination angle is perpendicular to the loading direction. The central angle of the panel has very significant effect on the buckling load. As the central angle decreases, the buckling load decreases. As the cutout area is increased, the buckling load of the panels (fixed and simple) with central angle of 0.5 radians is decreased. On the other hand, the cutout area has a slightly effect on the variation of panel buckling load with central angle of 0.1 radians. For instance, when the cutout area ratio is changed from 5 to 25%, the buckling load of the simple and fixed panels with central angle of 0.1 radians is changed about 15%. Moreover, when the cutout area ratio is reached to 25%, the effect of central angle of the panel on the buckling load of the panel is negligible. For instance, the buckling loads of fixed panel with central angle of 0.5 radian and simple panel with central angle of 0.1 radian are about 133 and 127 N/mm, respectively (in presence of a circular hole with cutout area ratio of 25%).

Figure 15
Cylindrical panel mode shape (a) layup [90]₈ and central angle of 0.1 (b) layup [90]₈ and central angle of 0.5 (c) layup [45]₈ and central angle of 0.1 (d) layup [45]₈ and central angle of 0.5 (e) layup [45/-45]₄ and central angle of 0.1 (f) layup [45/-45]₄ and central angle of 0.5 (g) layup [0/45/90/-45]_s and central angle of 0.1 (h) layup [0/45/90/-45]_s and central angle of 0.5.

The effect of boundary condition on the buckling load the panel with layup ${[0/90]}_4$ is shown in Figure 18. As illustrated in this figure, as the cutout ratio increases, the buckling load of the panel decreases. However, when the central angle of the panel is considered as 0.1 radians, the reduction is negligible for any boundary condition. Figure 19 shows the effect of boundary condition on the buckling load for layup ${[0/90]}_{2\mathrm{s}}$. The trend of variation of buckling load in this figure is same as Figure 18. The buckling loads of the fixed and simple panels with layups ${\left[45\right]}_8$ and ${[45/-45]}_4$ are presented in Figures 20 and 21, respectively. As an interesting point depicted in Figure 20, when the central angle of the fixed and simple panels is low (0.1 radians), the cutout area has not significant effect on the buckling load of the panel. In both panels (fixed and simple), considerable difference is seen between the buckling loads of the panel with central angles of 0.1 and 0.5 radian. However, as the cutout area increases, the difference is decreased, dramatically. As shown in Figure 19 (panel with layup ${[45/-45]}_4$), the cutout area has not a significant effect on the buckling load of the panel with central angle of 0.1 radians (fixed and simple panels). The buckling load of the panel is decreased when the central angle of the panel is 0.5 radians. When the panel layup is ${[30/-30]}_4$ (Figure 22), the trend of variation of the buckling load is same as Figure 20 (layup ${\left[45\right]}_8$). The panel with Layup ${[0/45/90/-45]}_{\mathrm{s}}$ is reached to the maximum and minimum in comparison with other layups (Figure 23). The trend of variation of the buckling load of the panel is same as Figure 20.

Figure 16
Effect of boundary conditions on the buckling load for layup ${\left[0\right]}_8$.

Figure 17
Effect of boundary conditions on the buckling load for layup ${\left[90\right]}_8$.

Figure 18
Effect of boundary conditions on the buckling load for layup ${[0/90]}_4$.

Figure 19
Effect of boundary conditions on the buckling load for layup ${[0/90]}_{2s}$.

Figure 20
Effect of boundary conditions on the buckling load for layup ${\left[45\right]}_8$.

Figure 21
Effect of boundary conditions on the buckling load for layup ${[45/-45]}_4$.

Figure 22
Effect of boundary conditions on the buckling load for layup ${[30/-30]}_4$.

Figure 23
Effect of boundary conditions on the buckling load for layup ${[0/45/90/-45]}_s$.

5 CONCLUSION

In this paper a sub-parametric spline finite strip element is developed to study the buckling analysis of perforated composite cylindrical panel. To this end, shape, mapping and displacement functions, Jacobian matrix, strain-displacement relations based on first order shear deformation theory, load vector, and linear and geometrical stiffness matrices of the strip element were derived. The obtained eigenvalue problem was solved by inverse subspace iteration to calculate the buckling load of the cylindrical panel. All of these steps were carried out in a FORTRAN code called SSFSM. The results obtained by the SSFSM were compared with those of the published literature and those obtained from ABAQUS finite element software to check the validity of the results.

A parametric study was performed to show the influences of cutout shape and area, panel thickness, central angle of the panel, boundary conditions of the panel and layup scheme on the axial buckling load of the composite cylindrical panel. The cutouts with the shape of circle, elliptic, square and diamond were considered in the center of the panel. The ratio of big to small chords of the elliptical cutout and the ratio of circular cutout area to panel area were varied. In addition eight layup configurations (angle-ply and cross-ply), three panel thicknesses (thick, medium and thin) and three panel central angles were considered in the parametric study. The following conclusions can be extracted based on the current study.

• - The simplicity of the Jacobian matrix in the current sub-parametric study causes the surface integrals encountered in the derivation of stiffness matrices turn into multiplication of two single variable integrals.

• - Another advantage of current sub-parametric formulation over iso-parametric one returns to the reduction in the required number of Gauss points to accurately estimate the domain integrals.

• - The reduction in execution time of current sub-parametric code against the standard iso-parametric formulation is about 60 percent. With increasing in panel degrees of freedom due to increasing in eigenvalue solvers contribution, this reduction is decreased slightly.

• - Increasing in the cutout chord prependicular to the loading direction while unaltering cutout area increases the axial buckling load of cylindrical panel.

• - Among the layups considered in the parametric study, quasi-isotropic laminate [0/45/90/-45]_s and uni-directional [0]₈ laminate carry the maximum and minimum axial buckling loads, respectively.

• - Fixed or simply supported curved loaded edges of the panel have not any effect on the axial buckling load of the panel for [90]₈ layup while they have the most effect in the case of [0]₈ arrangement.

• - Axial buckling mode of uni-directional laminates such as [90]₈ and [45]₈ has local expansion in the fibre direction and this is pronounced with increasing in the panel curvature. Unlike the uni-directional, in the multi-directional laminates like [45/-45]₄ and [0/45/90/-45]_s global expansion is illustrated.

References

• ABAQUS, (2006), User’s manual version 6.5. Hibbit, Karlsson & Sorensen, USA.
• Akin, J.E. (2005). Finite Element Analysis with Error Estimators, Elsevier publisher, New York, USA.
• Anil, V., Upadhyay, C.S. and Iyengar, N.G.R., (2007). Stability analysis of composite laminate with and without rectangular cutout under biaxial loading. Composite Structures 80: 92-104.
• Arbelo, M.A., Herrmann, A., Castro, S.G.P., Khakimova, R., Zimmermann, R., Degenhardt, R., (2015). Investigation of Buckling Behavior of Composite Shell Structures with Cutouts. Applied Composite Materials 22: 623-636.
• Asadi, E., Qatu, M.S., (2013). Free vibration of thick laminated cylindrical shells with different boundary conditions using general differential quadrature. Journal of Vibration and Control 19: 356-366.
• Asadi, E., Wang, W., Qatu, M.S., (2012). Static and vibration analyses of thick deep laminated cylindrical shells using 3D and various shear deformation theories. Composite Structures 94: 494-500.
• Ashrafi, H., Asemi, K., Shariyat, M., (2013). A three dimensional boundary element stress and bending analysis of transversely/longitudinally graded plates with circular cutouts under biaxial loading. European Journal of Mechanics A Solids 42: 344-357.
• Au, F.T.K., Cheung, Y.K. (1996). Free vibration and stability analysis of shells by the isoparametric spline finite strip method. Thin-Walled Structures 24: 53-82.
• Balogh, B., Logo, J. (2015). Optimization of curved plated structures with the finite strip and finite element methods. Civil Engeering Series 15: 1-10.
• Cheung, Y.K., Tham, L.G. (1997). The Finite Strip Method, CRC Press, New York, USA.
• Dawe, D.J., Wang, S., Lam, S.S.E., (1995). Finite strip analysis of imperfect laminated plates under end shortening and normal pressure. International Journal for Numerical Methods in Engineering 38: 4193-4205.
• Dawe, D.J., Wang, S., (1994). Buckling of composite plates and plate structures using the spline finite strip method. Composites Engineering 4: 1099-1117.
• Eccher, G., Rasmussen, K.J.R., Zandonini, R., (2008). Elastic buckling analysis of perforated thin-walled structures by the isoparametric spline finite strip method. Thin-Walled Structures 46: 165 - 191.
• Fazilati, J., (2017). Stability analysis of variable stiffness composite laminated plates with delamination using Spline FSM. Latin American Journal of Solids and Structures 14: 528- 543.
• Fazilati, J., Ovesy, H.R., (2013). Parametric instability of laminated longitudinally stiffened curved panels with cutout using higher order FSM. Composite Structures 95: 691-696.
• Khalfi, Y., Houari, M.S.A., Tounsi, A., (2014). A refined and simple shear deformation theory for thermal buckling of solar functionally graded plates on elastic foundation. International Journal of Computational Methods 11: 1350077.
• Khayat, M., Poorveis, D., Moradi, S., Hemmati, M. (2016a). Buckling of thick deep laminated composite shell of revolution under follower forces. Structural Engineering and Mechanics 58: 59-91.
• Khayat, M., Poorveis, D., Moradi, S., (2016b). Buckling analysis of laminated composite cylindrical shell subjected to lateral displacement-dependent pressure using semi analytical finite strip method. Steel Composite Structures 22: 45-59.
• Khetir, H., Bouiadjra, M.B., Houari, M.S.A., Tounsi, A., Mahmoud, S.R., (2017). A new nonlocal trigonometric shear deformation theory for thermal buckling analysis of embedded nanosize FG plates. Structural Engineering and Mechanics 64: 391-402.
• Komur, M.A., Sonmez, M., (2008). Elastic buckling of rectangular plates under linearly varying in-plane normal load with a circular cutout. Mechanics Research Communucations 35: 361-371.
• Magnucki, K., Mackiewicz, M., (2006). Elastic buckling of an axially compressed cylindrical panel with three edges simply supported and one edge free. Thin-Walled Structures 44: 387-392.
• Masia, U., Avalos, D.R., Laura, P.A.A., (2005). Displacement amplitudes and flexural moments for a rectangular plate with a rectangular cutout under a uniformly distributed static load. Journal of Sound and Vibration 280: 433-442.
• Mirzaei, S., Azhari, M., Bondarabady Rahimi, H.A., (2015). On the use of finite strip method for buckling anaylsis of moderately thick plate by refined plate theory and using new types of functions. Latin American Journal of Solids and Structures 12: 561- 582.
• Moradi, S., Poorveis, D., Khajehdezfuly, A., (2011). Geometrically nonlinear analysis of anisotropic laminated cylindrical panels with cut-out using spline finite strip method. Proceeding of conference on the Advances in Structural Engineering and Mechanics, Seoul, South Korea.
• Nejati, M., Dimitri, R., Tornabene, F., Hossein Yas, M. (2017). Thermal buckling of nanocomposite stiffened cylindrical shells reinforced by Functionally Graded Wavy Carbon Nano-Tubes with temperature-dependent properties. Applied Science 7: 1-24.
• Ovesy, H.R., Totounferoush, A., Ghannadpour, S.A.M., (2015). Dynamic buckling analysis of delaminated composite plates using semi-analytical finite strip method. Journal of Sound and Vibration 343: 131-143.
• Ovesy, H.R., Fazilati, J., (2012b). Parametric instability analysis of moderately thick FGM cylindrical panels using FSM. Computers and Structures 108: 135-143.
• Pascan, O.Z., Wei-hong, Z., Ponthot, J.P., (2012). Mechanical buckling analysis of composite panels with/without cutouts. International Journal of Plant Engineering Management 17: 65-76.
• Patel, B.P., Singh, S., Nath, Y., (2006). Stability and nonlinear dynamic behaviour of cross-ply laminated heated cylindrical shells. Latin American Journal of Solids and Structures 3: 245- 261.
• Rakotonirina, C., (2013). On the Cholesky method. Journal of Interdisciplinary Mathematics 12:875-882.
• Rajanna, T., Banerhee, S., Desai, M.Y., (2016). Effects of partial loading and fiber configuration on vibration and buckling characteristics of stiffened composite plates. Latin American Journal of Solids and Structures 13: 854- 879.
• Rondal, J., (1998). Coupled Instabilities in Metal Structures: Theoretical and Design Aspects, Springer publisher, New York, USA.
• Shuohui, Y., Tiantang, Y., Tinh, Q.B., Peng, L., Sohichi, H., (2016a). Buckling and vibration extended isogeometric analysis of imperfect graded Reissner-Mindlin plates with internal defects using NURBS and level sets. Computers & Structures 177: 23-38.
• Shuohui, Y., Tinh, Q.B., Shifeng, T., Satoyuki, T., Sohichi, H., (2016b). NURBS-based isogeometric analysis of buckling and free vibration problems for laminated composites plates with complicated cutouts using a new simple FSDT theory and level set method. Thin-Walled Structures 101: 141- 156.
• Shuohui, Y., Tiantang, Y., Tinh, Q.B., Shifeng, X., Sohichi, H., (2015). A cutout isogeometric analysis for thin laminated composite plates using level sets. Composite Structures 127: 152- 164.
• Soltani, K., Bessaim Houari, M.S.A., Kaci, A., Benguediab, M., Tounsi, A., Alhodaly, M., (2019). A novel hyperbolic shear deformation theory for the mechanical buckling analysis of advanced composite plates resting on elastic foundations. Steel and Composite structures 30:13-29.
• Soni, G., Singh, R., Mitra, M., (2013). Buckling behavior of composite laminates (with and without cutouts) subjected to non-uniform in-plane loads. International Journal of Structural Stability and Dynamics 13: 1-20.
• Tiantang, Y., Shuohui, Y. Tinh, Q.B., Chen, L., Nuttawit W., (2017). Buckling isogeometric analysis of functionally graded plates under combined thermal and mechanical loads. Composite Structures 162: 54-69.
• Thinh, T.I., Tu., T.M., Quoc, T.H., Long, N.V., (2016). Vibration and buckling analysis of functionally graded plates using new eight-uknown higher order shear deformation theory. Latin American Journal of Solids and Structures 13: 456- 477.
• Tiantang, Y., Tinh, Q.B., Duc, H.D., Wud, C.T., Thom, V.D., Satoyuki, T., (2016). On the thermal buckling analysis of functionally graded plates with internal defects using extended isogeometric analysis. Composite Structures 136: 684- 695.
• Van Ness, J., (1969). Inverse iteration method for finding eigenvectors. IEEE Transactions on Automatic Control 14: 63-66.
• Vuong, P.M., Dung, D.V., (2017). Nonlinear analysis on buckling and postbuckling FGM imperfect cylindrical shells filled inside by elastic foundations in thermal environment using TSDT. Latin American Journal of Solids and Structures 14: 950- 977.
• Wang, S., Chen, J., Dawe, D.J., (1998). Linear transient analysis of rectangular laminates using spline finite strips. Composites Structures 41: 57-66.
• Yao, L.K., He, B., Zhang, Y., Zhou, W., (2015). Semi-analytical finite strip transfer matrix method for buckling analysis of rectangular thin plates. Mathematical Problems in Engineering 15: 1-11.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  30 Sept 2019
• Date of issue
  2019

History

• Received
  01 Mar 2019
• Reviewed
  19 Aug 2019
• Accepted
  02 Sept 2019
• Published
  06 Sept 2019