Buckling analysis of laminated composite plates using an eﬃcient C 0 FE model

Buckling analysis of laminated composite plates is carried out by using an eﬃcient C 0 FE model developed based on higher order zigzag theory. In this model the ﬁrst derivatives of transverse displacement have been treated as independent variables to overcome the problem of C 1 continuity associated with the FE implementation of the plate theory. The C 0 continuity of the present FE model is compensated in the stiﬀness matrix calculations by using penalty parameter approach. Numerical results and comparison with other existing solutions show that the present model is very eﬃcient in predicting the buckling responses of laminated composites.


INTRODUCTION
Laminated composite plates are widely used in civil infrastructure systems due to their high strength to weight ratio and flexibility in design.One of the main failure mechanisms in composite plates is buckling.Accurate prediction of structural response characteristics is a challenging problem in the analysis of laminated composites due to the orthotropic structural behavior, the presence of various types of couplings and due to less thickness of the structural elements made of composites.Thus, an accurate buckling analysis of the laminated composite plates is an important part of the structural design.
Finite element method has been widely used for the buckling analysis of laminated composite plates.To predict the buckling load of composite plates, a number of plate theories have been proposed.The classical laminate plate theory (CLPT), which neglects the transverse shear deformation effect, yields acceptable results only for thin laminates [32].The structures designed based on CLPT theory may be unsafe because the CLPT overestimates the buckling load of the laminated composite plates.To take into account the effect of transverse shear deformation, the first-order shear deformation theory (FSDT) [1,28] has been used to predict the dynamic response of laminated composite structures.In the first-order shear deformation theory (FSDT) a shear correction factor is used to compensate for the assumed uniform transverse shear strain variations over the entire plate thickness.For a better representation of the transverse shear deformations, various higher-order shear deformation (HSDT) plate theories [5,15,16,30,36] were proposed.However, it has been observed and mentioned by many researchers [3] that increasing the number of terms in the in-plane displacement components does not improve the results and it is required to add the effect of inter-laminar transverse shear stress continuity in a multilayered composite plate problem.
A major development in this direction is due to Di Sciuva [34], Murakami [23], Liu et al. [20], and few others.They proposed zigzag plate theory where layer-wise theory is initially used to represent the in-plane displacements having piecewise linear variation across the thickness.The unknowns at the different interfaces are subsequently expressed in terms of those at the reference plane through satisfaction of transverse shear stress continuity at the layer interfaces.This theory is known as refined first order shear deformation theory (RFSDT).A further improvement in this direction is due to Di Sciuva [35], Bhaskar et al. [2], Cho et al. [8,16] and some other investigators.They included the condition of zero transverse shear stress at top and bottom of the plate in addition to the shear stress free conditions at the layer interfaces.This theory is known as higher order zigzag theory (HZT) or refined higher order shear deformation theory.
However, all these refined theories including the HSDT demand C 1 continuity of the transverse displacement along the element edges at the time of their finite element implementation.However, a C 1 formulation is not encouraged in any practical applications.In this situation it is highly required to develop a C 0 finite element formulation which will overcome the above C 1 continuity requirement of the theory (HZT).
Typically most of the C 0 refined theories [18,19,40] en-force the mentioned continuity assuming the displacement and transverse stresses as primary variable for multilayered composite plates employing the FE method and showed a good agreement with three dimensional elasticity solutions.On the other hand Ferreira et al. [39] combined the third order theory of Reddy [11] with a meshfree method based on the multi-quadric radial basis function approach to study the behavior of multilayered laminated composite beams and plates.Mondal et al. [29] analyzed the elastic stability behavior of simply supported anisotropic sandwich flat panels subjected to mechanical in-plane loads by using an analytical approach.Shufrin et al. [22] used a semi analytical approach for the buckling analysis of symmetrically laminated rectangular plates with general boundary conditions.Ganapathi et al. [37] investigated the free vibration characteristics of simply supported anisotropic composite laminates using analytical approach.Nali et al. [14] analyzed buckling behaviour of laminated plates subjected to combined biaxial/shear loading.Ferriera et al. [24] used meshfree method based on collocation with radial basis functions for buckling and vibration analysis of laminated composite plates.Ferriera et al. [12] used numerical method for the buckling analysis of laminated plates based on collocation with wavelets.Nguyen-Van et al. [10] presented buckling and vibration analysis of composite plate or shell structures of various shapes.Fiedler et al. [25] used refined higher-order shear deformation theory for buckling analysis of multi-layered plates subjected to unidirectional in-plane loads.Moreover, it has been also observed that in some cases the violation of inter-laminar shear stress continuity might lead to higher buckling load.
Latin American Journal of Solids and Structures 9(2012) 353 -365 In view of the above an efficient C 0 FE model has been developed based on a higher order zigzag theory (HZT) has been used to accurately predict the buckling load of laminated composite plates.The C 0 continuity of the present finite element model has been compensated in the stiffness matrix formulation by using penalty parameter approach.A nine noded C 0 continuous isoparametric finite element has been used for the development of the proposed finite element model.To assess the accuracy, numerical results obtained by using the proposed finite element model based on HZT have been compared with the results of 3D elasticity, exact and other finite element solutions available in the literature.

FORMULATION
The in-plane displacement fields (Figure 1) are taken as below: where u 0 α denotes the in-plane displacement of any point on the mid plane, n u and n l represent number of upper and lower layers respectively, S k α , T k α are the slopes of k -th layer corresponding to the upper and lower layers.ξ α , φ α are the higher order unknown terms, H(Z−Z k ), H(Z−P k ) are unit step functions and the subscript α represents the co-ordinate directions [α=1, 2 i.e.
x, y in this case ] respectively and The stres-sstrain relationship of a lamina, say k -th lamina having any fiber orientation with respect to the structural axes system (x-y-z ) may be expressed as: The rigidity matrix Qk can be evaluated by material properties and fibre orientation following usual techniques followed in case of laminated composites and {ε} is the strain field vector of size (5×1) at the reference plane (i.e., at the mid plane).Now by utilizing the transverse shear stress free boundary conditions (σ 3 α|z = ±h/2 = 0) at the top and bottom surface of the plate the components ξ α and φ α may be expressed as: and Similarly by imposing the transverse shear stress continuity conditions at the layer interfaces the following expressions for and are obtained: where, a k αγ , b k αγ , c k αγ , d k αγ , are constants depending on material and geometric properties of individual layers,W, γ is the derivatives of transverse displacement while γ = 1, 2 and S 0 α = ψ α is the rotation of normal at the mid plane about co-ordinate directions [α =1, 2 i.e., x, y in this case].
By using equations (2 -7) the strain field vector can be evaluated as below, where {ε} is the strain vector of size (17×1) at the reference plane (i.e., at the mid plane) where is the matrix of size (5×17) containing z terms and some terms related to material properties.
In present formulation based on higher order zigzag theory the in-plane displacement fields require C 1 continuity of the transverse displacement for its finite element implementation.In order to avoid the difficulties associated withC 1 continuity, the derivatives of w with respect to x and y are expressed as follows.
which helps to define all the variables including and as C 0 continuous.For the present study, a nine noded C 0 continuous isoparametric element shown in Figure 2 with seven nodal unknowns per node(u 1 , u 2 , w, Ψ 1 , Ψ 2 , w 1 , w 2 ) are used to develop the finite element formulation.The generalized displacements included in the present theory can be expressed as follows.
Latin American Journal of Solids and Structures 9(2012) 353 -365 where N i is the shape function for the nine noded isoparametric element [33].Employing Hamilton's principle, the dynamic equation of equilibrium for a finite element can be formed.Accordingly the element stiffness matrix can be written as below in which B is the strain matrix and Q i is the transformed material constant matrix.where In a similar manner the geometric strain vector may be expressed as, With the matrix [G] in the above equation, the geometric stiffness matrix of an element can be derived and may be written as where [S k ] is the in-plane stress components of the k -th layer Combining all the element geometric and elastic stiffness matrices, the dynamic equilibrium equation can be written as in which [K] and [K G ] are the global elastic and geometric stiffness matrix, where[∂] is the mode shape vector of buckling and λ is the critical buckling load parameter respectively.

RESULTS AND DISCUSSION
In this section some problems on buckling of laminated composite plates are solved using the present finite element model considering different features such as boundary conditions, ply orientations, thickness ratio and aspect ratio.The results obtained by using the proposed FE model based on HZT is first validated with some published results and then many new results are also generated.

Simply supported square isotropic plates
In this example simply supported square isotropic (ν= 0.3) plate subjected to uni-axial loading is considered.The analysis is carried out for different thickness ratio (a/h = 100, 10).The critical buckling load obtained by using present FE model obtained are presented with the results of Chakrabarti and Sheikh [6] based on refined higher order shear deformation theory (RFSDT), Sundaresan et al. [38] based on first order shear deformation theory and that of general purpose finite element program MARC [38] in Table 1.It is observed that the present results are somewhat lesser than Chakrabarti and Sheikh [6], Sundaresan et al. [38] and general purpose finite element program MARC [38].It may be due to efficient modeling of the transverse shear deformation followed in the present finite element model.Chakrabarti and Sheikh [6] 4.0 Sundersan et al. [38] 4.0 MARC [38] 4.0 10 Present(12×12 1 ) 3.768 Chakrabarti and Sheikh [6] 3.782

Cross-ply multilayered (0/90/..) composite plate simply supported at all the edges
In this example the buckling of a laminated (0 ○ /90 ○ /90 ○ /0 ○ ) square composite plate is considered.The material properties of the individual layers are given by: This problem is solved to assess the performance of the proposed C 0 finite element model.The convergence and comparison of the normalized uni-axial critical buckling load obtained by using the proposed element is shown in Table 2.The thickness ratio (a/h) is varied (100 to 5) to cover the range of very thin to thick plates.It is found that with refining meshes, results obtained by the present FE model converge closely to the exact solution [26] and reasonably close to the 3D elasticity solution given by Noor [13].The present finite element model based on higher order zigzag theory gives the more accurate results compared to the results given by Fiedler et al. [33] based on a similar theory, results of Van et al. [25] based on first order shear deformation theory and results of Ferreira et al. [10] based on first order shear deformation theory respectively.Some new results are also presented in Table 2.
In the next example, the effect of number of layers and degree of orthotropy of individual layers on normalized uni-axial critical buckling load has been studied.The numbers of layers are varied from 3 to 5 and thickness ratio (a/h) is considered as 10 where h is the overall thickness of the plate.The material properties are taken same as in the previous section.The results are given in Table 3.It is to be noted that the present results are somewhat lesser than 3D elasticity solution of Noor [13] as compared to Fiedler et al. [33], Ferreira et al. [10] and Reddy and Phan [27] respectively.It has been observed in the previous example that the present finite element model based on higher order zigzag theory gives more accurate results compared to the results based on various other theories.It may be due to accurate modeling of the transverse shear deformation followed in the present finite element model to incorporate the exact flexible pattern of the structures.Due to inclusion of shear deformation in a refined manner the stiffness of the structure becomes less and therefore, the buckling load calculated by the present model gives lesser values than those obtained by other theories.

Cross-ply (0 ○ /90 ○ ) laminate having different thickness ratio (a/h)
The effect of span to thickness ratio (a/h) on the uni-axial critical buckling load for simply supported two layer cross-ply (0 ○ /90 ○ ) square plate is studied in this example.The material properties are same as in section 3.2 where E 1 /E 2 is taken as 40.The results obtained by the present finite element model are shown in Table 4 along with those obtained by Van et al. [25], Reddy and Phan [27], and Chakrabarti and Sheikh [31] respectively.The same trend as explained before is also observed in the results presented in Table 4.The present results based on higher order zigzag theory are on the lower side compared to the results base on other theories as expected.In this example the effect of aspect ratio (a/b) on the uni-axial critical buckling load for simply supported four layer cross-ply (0 ○ /90 ○ /90 ○ /0 ○ ) square plate is studied.The material properties are same as in the previous example.The results obtained by using the proposed model are shown in Table 5 along with those obtained by Fiedler et al. [13].The results are in good agreement due the fact that the theory used in the present finite element model and in the analytical approach adopted by Fiedler et al. [13] is nearly same.Some new results are also presented in Table 5.

Effect of Different boundary conditions
The influence of different boundary conditions on normalized critical uni-axial buckling load (λ cr = λ a 2 E 2 h 3 ) is investigated in the present section.The plate is considered simply supported (S) along the edges parallel to the y axis while the other edges have simply supported (S), clamped (C) or free (F) boundary conditions.The notation SSCF, for example refers to the simply supported boundary conditions of the two edges parallel to the y axis and the clamped and free conditions for the two edges parallel to the x axis.The eight layer (0 ○ /90 ○ ) 4 square plate is analyzed with material properties same as in previous example while the thickness ratio (a/h) has been taken as 10 and 100.The normalized critical uni-axial buckling loads of the eight layer (0 ○ /90 ○ ) 4 plate under different boundary conditions are shown in Table 6.All these results are presented for the first time.The effect of modular ratio on normalized critical bi-axial buckling load (λ cr = λ a 2 E 2 h 3 ) for simply supported cross-ply [0 ○ /90 ○ /0 ○ ] square plate with thickness ratio (a/h = 10) is studied in this example.The material properties are same as in section 3.2.The results obtained by the present element are shown in Table 7 along with other published results.The results obtained by the proposed finite element model agree reasonably well with the results obtained by Van et al.,FSDT [25], Ferreira et al., FSDT [10] and Khedir and Librescu, HSDT [26] for lower E 1 /E 2 however with higher value of this ratio i.e., with increased orthotrpic behavior the present results are on the lower side as expected and as explained earlier.No pre-buckling state is considered since other results are also based on the same assumption.In the present example of anti-symmetric cross-ply (0 ○ /90 ○ /0 ○ /90 ○ ) laminate (E 1 /E 2 =40) firstly the effect of thickness ratio (a/h) on normalized critical ui-axial and bi-axial buckling load (λ cr = λ a 2 E 2 h 3 ) with simply supported boundary conditions is studied.The material properties are same as in section 3.2.The comparison between the normalized critical ui-axial and bi-axial buckling load is shown in Figure 2. The comparison shows expected trend.
Next the effect of modular ratio (E 1 /E 2 ) on normalized critical ui-axial and bi-axial buckling load (λ cr = λ a 2 E 2 h 3 ) is studied with thickness ratio (a/h) = 10.The material properties are same as in section 3.1.The comparison between the normalized critical ui-axial and bi-axial buckling load is shown in Figure 3.In this section the effect of thickness ratio (a/h) on normalized critical pure shear buckling load (λ cr = λ × 10 −3 ) is studied taking E 1 /E 2 =40 for an angle-ply (45 ○ /-45 ○ /-45 ○ /45 ○ ) laminate with simply supported boundary conditions.The material properties are same as in section 3.2.The variation of the normalized critical pure shear buckling load with different thickness ratio (a/h) is shown in Figure 4.

CONCLUSIONS
In the present study, the buckling behaviour of laminated composite plates is studied using an efficient C 0 finite element model developed based on higher order zigzag theory.Many problems are solved covering different features of laminated composites such as boundary conditions, ply orientations aspect ratio, thickness ratio and loading etc. Numerical results on buckling responses for different problems show that the proposed FE model for laminated composite plates is capable in predicting results very close to the exact solutions.Many new results are also generated.As such the model may be recommended for wide use in research and industrial applications.

Figure 1
Figure 1 Geometry of laminated composite plate

Figure 4
Figure 4 Variation of normalized critical shear buckling load with different thickness ratio (a/h)

Table 1
Normalized Critical buckling loads (λcr) for square isotropic plate

Table 3
Normalized critical uni-axial buckling loads (λcr) with various modular ratios for simply supported cross-ply square plate with thickness ratio (a/h = 10)