Bi-Axial Buckling of Laminated Composite Plates Including Cutout and Additional Mass

In the presented paper, a study of bi-axial buckling of the laminated composite plate with mass variation through the cutout and additional mass is carried out using the improved shear deformation theory (ISDT). The ISDT mathematical model employs a cubic variation of thickness co-ordinates in the displacement field. A realistic parabolic distribution of transverse shear strains through the plate thickness is assumed and the use of shear correction factor is avoided. A C° finite element formulation of the mathematical model is developed to analyze the buckling behavior of laminated composite plate with cutout and additional mass. As no results based on ISDT for the considered problem of bi-axial buckling of the laminated composite plate with mass variation are available in the literature, the obtained results are validated with the data available for a laminated composite plate without cutout and additional mass. Novel results are obtained by varying geometry, boundary conditions and ply orientations.


Introduction
The composite materials proved to be much more efficient than the traditional materials that led to its wide use as a structural element around the world. Nowadays fiber reinforced polymer composites constitute the dominant materials applied in industry, e.g., the aerospace industry. Composites consist of a combination of two or more materials with different properties. They contain various fibers, which may be metallic or non-metallic, as well as grains of various materials and flakes. Distribution of these components may be homogenous or non-homogeneous. Out of all composite materials, laminated composites are very popular in industries. Excellent properties, such as high strength, high strength to weight ratio and low weight make them more popular in contemporary industrial environments [1]. They are characterized by high stiffness, damping and good directional properties. The fiber in composites is the main load-bearing constituent which is stronger and stiffer than other structural materials. For this reason, composites are used in structural elements, marine parts, aircrafts, etc.
Certain properties, such as static, dynamic and damping are very important [1,2]. In addition, there are notches in these laminated boards, because the cutouts can be used as doors and windows, presented numerical and experimental studies on buckling of skew laminates with circular cutouts under uniaxial compressive force.
The literature review indicates that few studies on biaxial buckling analysis of laminated composite panel with excision and concentrated mass were carried out and that significant part of the reported studies was based on the FSDT. Due to the complex mechanical behavior of such as structures, which, for example, exhibit complex modes of deformation, their theoretical analysis requires a slightly deeper approach.
In the presented work, an attempt was made to analyze the bi-axial buckling characteristics of laminated composite plates with cutouts and additional mass using a mathematical model (ISDT) which employs a cubic variation of thickness co-ordinate in displacement field. In the present study, an attempt was also made to incorporate shear along with bi-axial force for analyzing the buckling behavior of laminated plates.

Improved Shear Deformation Theory (ISDT)
In this deformation theory, transverse shear stress at the top and bottom of the laminate are taken as zero. It is assumed that the variation of transverse shear strains is realistic parabolic in shape and use of shear correction factor is hence avoided. The presented theory consists of a realistic cubic variation of in-plane displacement fields. The following equation for displacement fields is being adopted for the presented analysis: In the equation presented above, u, v and w represent the displacements of a point along the three directions (x, y, and z) respectively, whereas the associated midplane displacements are given by u 0 , v 0 , and w 0 respectively. θ x and θ y signifies the rotations of transverse normal in the x-z and y-z planes respectively.

of 23
When replaced, the values obtained in Equation (4) to Equation (1), we get: The linear strains may be represented in the form of linear displacement, as follows: ∂u/∂x ∂v/∂y ∂u/∂y + ∂v/∂x ∂u/∂z + ∂w/∂x ∂v/∂z + ∂w/∂y Using the values of displacement from Equation (5) to Equation (6), the following equation is obtained, The strains associated with Equation (8) are related to the generalized strains by means of the following expression: where {ε} = ε x ε y γ xy γ xz γ yz where For the typical lamina (kth), the constitutive relations with respect to the material axis can be expressed as: The stress-strain relationship with respect to global coordinate axis system (x, y, and z) for kth lamina can be expressed using the applied transformation coefficients as shown below: Integration of the stresses through the laminate thickness will help in obtaining the resultant forces and moments acting on the laminate, which is as follows:

Finite Element Formulations
In this paper, C • isoparametric elements having nine nodes, with seven unknowns per node (Figure 1) i.e., u 1 u 2 , u 3 , ψ 1, ψ 2 , w 1 and w 2 were used for the proposed finite element model. Thus, by following the standard procedure of FEM the element matrices are assembled which results in global stiffness matrices i.e., [K] and [M].

Finite Element Formulations
In this paper, C° isoparametric elements having nine nodes, with seven unknowns per node (Figure 1) i.e., were used for the proposed finite element model.  The generalized displacements included in the presented theory can be expressed as follows. (13) In the equation above, the shape function of the related node is represented by N i . Knowing the nodal unknown vector within an element, the mid-surface strains at any point in the plate can be expressed in the matrix form in terms of the global displacements as described below: In the equation shown above, [B i ] is used to express the differential operator matrix of the shape function.
For an element, the element stiffness matrix (say, eth), including the transverse shear effects, flexure and membrane can be given as: (15) [ In the above equations, [N] represents the matrix of shape function, [ρ] inertia matrix and |J| represents the determinant of a Jacobian matrix.
For all numerical integrations, a 3 × 3 Gaussian quadrature format was used. Then, the element matrices were grouped together to attain the global stiffness matrices [K], in accordance with the typical procedure of the FE method as considered by Bathe [54].
The C 0 FE formulation although being a 2D solutions allowed to obtain the results closer to 3D elastic solutions. This helped to reduce the complex solutions to a simpler form.
The geometric stiffness matrix [K Ge ] of an element can be developed using Equations (15) and (18) and it may be stated as: where, (20) and S i is in-plane stress components of the i-th layer.
Hence the final governing equation for buckling analysis may be written as: where {δ} is the nodal displacement vector, λ is the critical buckling load and [K], [K] G are the linear and geometric stiffness matrices, respectively.

Boundary Conditions
The boundary conditions mainly applied to the subsequent examples are, simply supported (SSSS) and clamped boundary conditions (CCCC) which include:

Engineering Properties
For all further investigations, unless mentioned otherwise, the composites with following properties were taken E 1 = 25E 2 , E 2 =1, G 12 = G 13 = 0.5E 2 , G 23 = 0.2E 2 , ρ = 1 and ν 12 = 0.25, after Chakrabarti and Sheikh [46] as the common standard values for these materials. The non-dimensional buckling load for composites were taken as λ = (λb 2 /E 2 h 3 ). The value of additional mass is given by the following equation:

Convergence and Comparison Studies
Convergence study was done to determine the required mesh size N × N at which the dimensionless buckling load values converged. It may be concluded from Table 1, that the values of non-dimensional buckling load converged at N = 20. Therefore, for all subsequent analyses mesh size of 20 × 20 (full) was taken into consideration. The numbering of edges in the plate is shown in Figure 2. Table 1. Convergence study of nondimensional buckling load λ = (λb 2 /E 2 h 3 ) with a/h for a simply supported cross-ply square laminated plate (0 • /90 • /90 • /0 • ) t with material properties:

Convergence and Comparison Studies
Convergence study was done to determine the required mesh size N × N at which the dimensionless buckling load values converged. It may be concluded from Table 1, that the values of non-dimensional buckling load converged at N = 20. Therefore, for all subsequent analyses mesh size of 20 × 20 (full) was taken into consideration. The numbering of edges in the plate is shown in Figure 2.
with a/h for a simply supported cross-ply square laminated plate (0°/90°/90°/0°)t with material properties: It is clear from the literature review that based upon present theory, for the buckling analysis of the laminated composite plates having mass variation in the form of cutout and concentrated mass no results are available. Therefore, in order to show the efficiency of the present FE model the obtained results were evaluated with the results published by Reddy and Phan [41], Pandit et al. [55], Liu et al. [56], and Singh et al. [57] based upon FSDT, HSDT, GRBF (Gaussian radial basic function) and MQRBF (Multiquadric radial basic function). Analysis of cross-ply (0°/90°/90°/0°) and (0°/90°/0°) square laminates under the effect of uni-axial compression was carried in this example and shown in Table 2. In this example, the analysis of a full plate was done with a/h ratio equal to 10. Table 2 shows, the utility of the present FE model in predicting the non-dimensional buckling load close to the analytical results from the previously It is clear from the literature review that based upon present theory, for the buckling analysis of the laminated composite plates having mass variation in the form of cutout and concentrated mass no results are available. Therefore, in order to show the efficiency of the present FE model the obtained results were evaluated with the results published by Reddy and Phan [41], Pandit et al. [55], Liu et al. [56], and Singh et al. [57] based upon FSDT, HSDT, GRBF (Gaussian radial basic function) and MQRBF (Multiquadric radial basic function).

References
Analysis of cross-ply (0 • /90 • /90 • /0 • ) and (0 • /90 • /0 • ) square laminates under the effect of uni-axial compression was carried in this example and shown in Table 2. In this example, the analysis of a full plate was done with a/h ratio equal to 10. Table 2 shows, the utility of the present FE model in predicting the non-dimensional buckling load close to the analytical results from the previously quoted literature. Table 2. Validation of nondimensional buckling load λ = (λb 2 /E 2 h 3 ) for a uni-axial buckling of simply supported cross-ply square laminated plate: G 12 = G 23 = 0.6E 2 , G 23 = 0.5E 2 , ν 12 = 0.25 and a/h = 10.

Lamination
Scheme In another problem, a square laminated composite plate having lamination scheme as (0 • /90 • /0 • ) was analyzed for different values of elastic modulus ratio i.e., E 1 /E 2 . The results were compared with the results of Vesconi and Dozio [30] as shown in Table 3. Table 3. Validation of nondimensional buckling load λ = (λb 2 /E 2 h 3 ) for a bi-axial buckling (N x /N y ) of simply supported (0 • /90 • /0 • ) cross-ply square laminated plate: The aspect ratios (a/h) taken into consideration for the study were 10, 20 and 50. Table 3 shows that for both ratios of N x /N y the results are in good accordance with those reported by Vescovini and Dozio [30]. Buckling load for two sides clamped laminates compared to experimental and FEM (ANSYS) results of Baba [58] are shown in Table 4.

Novel Results
After validating the presented FE model based upon the above-mentioned theory through comparison studies, new problems were worked out to analyze the effect of openings and additional mass on buckling of laminated composite plates. For the following examples of buckling analyses, various laminated composite plates with different lamination schemes and boundary conditions were considered. The lamination schemes adopted for the various problems are given in Table 5. The geometrical and material parameters adopted for the present analyses were assumed as defined in the previous section. The ply numbering scheme shown in Figure 3, is in such a way that the counting of the lamina was done from bottom to top. In the present examples, the laminated composite plates are considered having the additional mass (M = M ρha 2 = 0.5, 1 and 2) and square cutout (0.2a, 0.4a and 0.6a), concentrated at the center. Table 5. Plate lamination schemes.

Symbol
Lamination Scheme * were considered. The lamination schemes adopted for the various problems are given in Table 5. The geometrical and material parameters adopted for the present analyses were assumed as defined in the previous section. The ply numbering scheme shown in Figure 3, is in such a way that the counting    , 1 and 2), lamination schemes and boundary conditions. Further, the mode shapes of plates were also shown in Figures 4-6. In all the preceding problems the material properties were taken as defined in the previous section and the additional mass was concentrated at the central node.

Laminated Composite Plates with Additional Mass
Many novel problems were solved and shown in Tables 6-8 with variation in values of aspect ratio (a/h), nature of loading (N x /N y & N xy ), additional mass M (0.05, 1 and 2), lamination schemes and boundary conditions. Further, the mode shapes of plates were also shown in Figures 4-6. In all the preceding problems the material properties were taken as defined in the previous section and the additional mass was concentrated at the central node.   dimensional buckling load parameter increases with a decrease in thickness. It may be also observed that for laminated plates having a/h ratio greater than 5, the CCCC boundary condition has the greater values of non-dimensional buckling loads (this is due to the highest stiffness of CCCC boundary conditions) whereas the lowest for CFCF boundary condition. Whereas for an a/h ratio equal to 5, any loading condition has almost the same nondimensional buckling load for all boundary conditions except for CFCF boundary condition. The nondimensional buckling load was also found to have the least variation for increasing values of additional mass for a particular aspect ratio and applied loading conditions. In Tables 7a,b, and 9a,b, the other problems were solved for the composite plates having lamination scheme B and C with different boundary conditions. In the present problem, different values of aspect ratio (a/h = 100, 20 and 5) and Nx/Ny (Nx/Ny = 0, 1 and 2) were taken in consideration. The values of shear forces (Nxy) were taken as 0, 1 and 2. In the tables, all trends of variations for non-dimensional buckling load with respect to different loading and boundary conditions are similar to the ones presented in Table 6a,b. It may be concluded from the tables that for any boundary condition different than CFCF and CCFF as the number of plies are increased the values of non-dimensional buckling load increases along with any particular value of aspect ratio, loading condition and additional mass. The mode shapes for the buckled laminated composite plates are shown in Figures 4-6. Figure  4 shows mode shapes for composite plates having lamination scheme A and clamped at all edges. The value of additional mass and a/h ratio was taken as 0.05 and 100, respectively. The different conditions of uni-axial and bi-axial loading with or without shear were taken into consideration. Figure 4 shows that mode shapes for bi-axial loading are different from uni-axial buckling whereas mode shapes for different values of bi-axial loading remain the same, only the values are changed. Figure 5 shows mode shapes for different lamination schemes keeping other conditions like a/h ratio, nature of the load, the values of additional mass and boundary conditions are constant. It can be seen that on changing the lamination scheme, the mode shape changes; this may be due to the fact that a different orientation of fiber in plies affects the strength of laminate and thereby changing the mode of deformation.         In Table 6a,b, the non-dimensional buckling loads for composite plates having lamination scheme A and different boundary conditions are shown. In the present problem, different values of aspect ratio (a/h = 100, 20 and 5) and N x /N y (N x /N y = 0, 1 and 2) are taken in consideration. The values of shear forces (N xy ) are taken as 0, 1 and 2. It is observed from the table that the non-dimensional buckling loads are minimum for lower values of a/h and higher values of biaxial (N x /N y ) and shear (N xy ) forces. The buckling load on plates is proportional to h 2 and buckling load should decrease with a decrease in thickness. But the presented results are in terms of the non-dimensional buckling load parameter (λ = N x a 2 /(E 2 h 3 )) which is inversely proportional to cubed thickness. Hence, the non-dimensional buckling load parameter increases with a decrease in thickness. It may be also observed that for laminated plates having a/h ratio greater than 5, the CCCC boundary condition has the greater values of non-dimensional buckling loads (this is due to the highest stiffness of CCCC boundary conditions) whereas the lowest for CFCF boundary condition. Whereas for an a/h ratio equal to 5, any loading condition has almost the same nondimensional buckling load for all boundary conditions except for CFCF boundary condition. The nondimensional buckling load was also found to have the least variation for increasing values of additional mass for a particular aspect ratio and applied loading conditions. In Table 7a,b, and Table 9a,b, the other problems were solved for the composite plates having lamination scheme B and C with different boundary conditions. In the present problem, different values of aspect ratio (a/h = 100, 20 and 5) and N x /N y (N x /N y = 0, 1 and 2) were taken in consideration. The values of shear forces (N xy ) were taken as 0, 1 and 2. In the tables, all trends of variations for non-dimensional buckling load with respect to different loading and boundary conditions are similar to the ones presented in Table 6a,b. It may be concluded from the tables that for any boundary condition different than CFCF and CCFF as the number of plies are increased the values of non-dimensional buckling load increases along with any particular value of aspect ratio, loading condition and additional mass. Table 9. (a) Non-dimensional buckling load λ = (λb 2 /E 2 h 3 ) for a bi-axial buckling of a composite plate having lamination scheme A and central square cutout; (b) Non-dimensional buckling load λ = (λb 2 /E 2 h 3 ) for a bi-axial buckling of a composite plate having lamination scheme A and central square cutout.   Figure 4 shows mode shapes for composite plates having lamination scheme A and clamped at all edges. The value of additional mass and a/h ratio was taken as 0.05 and 100, respectively. The different conditions of uni-axial and bi-axial loading with or without shear were taken into consideration. Figure 4 shows that mode shapes for bi-axial loading are different from uni-axial buckling whereas mode shapes for different values of bi-axial loading remain the same, only the values are changed. Figure 5 shows mode shapes for different lamination schemes keeping other conditions like a/h ratio, nature of the load, the values of additional mass and boundary conditions are constant. It can be seen that on changing the lamination scheme, the mode shape changes; this may be due to the fact that a different orientation of fiber in plies affects the strength of laminate and thereby changing the mode of deformation. Figure 6 presents the variation of mode shape for different boundary conditions. The figure indicates that mode shapes are affected by variation in the boundary condition. The least distorted mode shape is observed for SSSS boundary condition.

Laminated Composite Plates with Central Cutout
Many novel problems were solved and shown in Tables 9-11 taking different values of aspect ratio (a/h), nature of loading (N x /N y & N xy ), cutout sizes (0.2a × 0.2a, 0.4a × 0.4a and 0.6a × 0.6a), lamination schemes and boundary conditions. Further, the mode shapes were also drawn and shown in Figures 7-9. In all the preceding problems the material properties were assumed as defined in the previous section and a square cutout is taken at the center of the plate  In Table 9a,b, the non-dimensional buckling loads for the composite plates having lamination scheme A and different boundary conditions are shown. In the present problem, the different values of aspect ratio (a/h = 100, 20 and 5) and N x /N y (N x /N y = 0, 1 and 2) were taken in consideration. The values of shear forces (N xy ) were taken as 0, 1 and 2. The table shows that for lower values of aspect ratio the non-dimensional buckling load is minimum. The values of loads are also found to be lower for higher values of biaxial (N x /N y ) and shear (N xy ) forces. It may be also observed, that for laminated plates having the a/h ratio greater than 5, CCCC boundary conditions have the greater values of the non-dimensional buckling loads whereas the lowest is for CFCF boundary conditions. Whereas for a/h ratio equal to 5, any loading condition has almost the same nondimensional buckling load for all the boundary conditions except for CFCF boundary condition. The non-dimensional buckling load is found to increase along with the size of the central cutout area.
In Tables 10 and 11, other solved problems are presented for composite plates having lamination scheme B and C with different boundary conditions. In the present problem, different values of aspect ratio (a/h = 100, 20 and 5) and N x /N y (N x /N y = 0, 1 and 2) were taken in consideration. The values of shear forces (N xy ) were taken as 0, 1 and 2. In the tables, all trends of variations for non-dimensional buckling load with respect to different loading and boundary conditions were similar to the ones presented in Table 9. The values presented in the discussed tables show that for any boundary condition the values of non-dimensional buckling load increase along with the number of plies for any particular value of aspect ratio, loading condition and cutout size.         The mode shapes for the central area of buckled laminated composite plates are shown in Figures 7-9. The plane section area with no deformations represents the cutout in the plate. As there is no material at the cutout in the plate, it does not move in the mode shape. Figure 8 shows mode shapes for central area of composite plates clamped at all edges and having lamination scheme A. The value of a central square cutout and a/h ratio were taken as 0.4a and 100, respectively. The different conditions of uni-axial and bi-axial loading with or without shear were taken into consideration. From Figure 7, it can be seen that mode shapes for bi-axial loading are different from uni-axial buckling, whereas the mode shapes for different values of bi-axial loading remain same, only the values are different. In Table 9a,b, the non-dimensional buckling loads for the composite plates having lamination scheme A and different boundary conditions are shown. In the present problem, the different values of aspect ratio (a/h = 100, 20 and 5) and Nx/Ny (Nx/Ny = 0, 1 and 2) were taken in consideration. The values of shear forces (Nxy) were taken as 0, 1 and 2. The table shows that for lower values of aspect ratio the non-dimensional buckling load is minimum. The values of loads are also found to be lower for higher values of biaxial (Nx/Ny) and shear (Nxy) forces. It may be also observed, that for laminated plates having the a/h ratio greater than 5, CCCC boundary conditions have the greater values of the non-dimensional buckling loads whereas the lowest is for CFCF boundary conditions. Whereas for a/h ratio equal to 5, any loading condition has almost the same nondimensional buckling load for all the boundary conditions except for CFCF boundary condition. The non-dimensional buckling load is found to increase along with the size of the central cutout area.
In Tables 10 and 11, other solved problems are presented for composite plates having lamination scheme B and C with different boundary conditions. In the present problem, different values of aspect ratio (a/h = 100, 20 and 5) and Nx/Ny (Nx/Ny = 0, 1 and 2) were taken in consideration. The values of shear forces (Nxy) were taken as 0, 1 and 2. In the tables, all trends of variations for non-dimensional buckling load with respect to different loading and boundary conditions were similar to the ones presented in Table 9. The values presented in the discussed tables show that for any boundary condition the values of non-dimensional buckling load increase along with the number of plies for any particular value of aspect ratio, loading condition and cutout size.
The mode shapes for the central area of buckled laminated composite plates are shown in Figures 7-9. The plane section area with no deformations represents the cutout in the plate. As there is no material at the cutout in the plate, it does not move in the mode shape. Figure 8 shows mode shapes for central area of composite plates clamped at all edges and having lamination scheme A. The value of a central square cutout and a/h ratio were taken as 0.4a and 100, respectively. The different conditions of uni-axial and bi-axial loading with or without shear were taken into consideration. From Figure 7, it can be seen that mode shapes for bi-axial loading are different from uni-axial buckling, whereas the mode shapes for different values of bi-axial loading remain same, only the values are different. Figure 8 shows the mode shapes for different lamination schemes obtained preserving other conditions, like a/h ratio, loading condition, cutout size and boundary conditions constant. It is visible that changes in lamination scheme trigger changes in the mode shapes, this may be due to the fact that the different orientation of fiber in plies affects the strength of laminate and thereby changing the mode of deformation. Figure 9 shows the variation of mode shape for different boundary conditions keeping other parameters, like cutout size, aspect ratio and nature of loading, constant. The figure indicates that mode shapes are affected by variation in the boundary conditions.  Figure 8 shows the mode shapes for different lamination schemes obtained preserving other conditions, like a/h ratio, loading condition, cutout size and boundary conditions constant. It is visible that changes in lamination scheme trigger changes in the mode shapes, this may be due to the fact that the different orientation of fiber in plies affects the strength of laminate and thereby changing the mode of deformation. Figure 9 shows the variation of mode shape for different boundary conditions keeping other parameters, like cutout size, aspect ratio and nature of loading, constant. The figure indicates that mode shapes are affected by variation in the boundary conditions.

Conclusions
In this paper, using the presented ISDT formulation and C • finite element model a computer code for Uni-axial and Bi-axial buckling analysis of laminated composite plates were developed. The proposed FE model based on the presented theory was analyzed for the first time. Along with this theory, the transverse shear stress continuity was also incorporated at the interface of each layer in addition to zero transverse shear stress at the top and bottom of the plate. The performed analyses showed that the obtained results are much improved over the other existing models (FSDT and HSDT). Various novel problems with different geometrical properties, loading, and boundary conditions and ply orientation were analyzed for the laminated composite plates. The following conclusions were drawn from the presented study: