Buckling analysis of laminated composite plates using an efficient C0 FE model

Abstract

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

Buckling; laminated composites; Finite element analysis

Buckling analysis of laminated composite plates using an efficient C⁰ FE model

S. K. Singh^*; A. Chakrabarti

Department of Civil Engineering, Indian Institute of Technology, Roorkee-247 667, India Ph.+91(1332)285844, Fax.[phone]

ABSTRACT

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

Keywords: Buckling, laminated composites, Finite element analysis

1 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¹ continuity of the transverse displacement along the element edges at the time of their finite element implementation. However, a C¹ formulation is not encouraged in any practical applications. In this situation it is highly required to develop a C⁰ finite element formulation which will overcome the above C¹ continuity requirement of the theory (HZT).

Typically most of the C⁰ 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.

In view of the above an efficient C⁰ 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⁰ continuity of the present finite element model has been compensated in the stiffness matrix formulation by using penalty parameter approach. A nine noded C⁰ 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.

2 FORMULATION

The in-plane displacement fields (Figure 1) are taken as below:

where 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, , 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 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 (σ₃α|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, , are constants depending on material and geometric properties of individual layers,W,_γ is the derivatives of transverse displacement while γ = 1, 2 and = ψ_α 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 zterms and some terms related to material properties.

In present formulation based on higher order zigzag theory the in-plane displacement fields require C¹ continuity of the transverse displacement for its finite element implementation. In order to avoid the difficulties associated with C¹ 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⁰ continuous.

For the present study, a nine noded C⁰ continuous isoparametric element shown in Figure 2 with seven nodal unknowns per node(u₁, u₂, w, Ψ₁, Ψ₂, w₁, w₂) are used to develop the finite element formulation. The generalized displacements included in the present theory can be expressed as follows.

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

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

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

3.2 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: E₁/E₂ = open, G₁₂ = G₁₃ = 0.6E₂, G₂₃ = 0.5E₂, ν₁₂ =ν₁₃ = 0.25.

This problem is solved to assess the performance of the proposed C⁰ 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.

3.3 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₁/E₂ 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.

3.4 Cross-ply (0⁰/90⁰/90⁰/0⁰) laminate having different aspect ratio (a/b)

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.

3.5 Effect of Different boundary conditions

The influence of different boundary conditions on normalized critical uni-axial buckling load () 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⁰)₄ 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⁰)₄ plate under different boundary conditions are shown in Table 6. All these results are presented for the first time.

3.6 Cross-ply (0⁰/90⁰/0⁰) laminate having different modular ratio (E₁/E₂) under bi-axial loading

The effect of modular ratio on normalized critical bi-axial buckling load () 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₁/E₂ 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.

3.7 Anti-symmetric cross-ply (0⁰/90⁰/0⁰/90⁰) laminate

In the present example of anti-symmetric cross-ply (0⁰/90⁰/0⁰/90⁰) laminate (E₁/E₂ =40) firstly the effect of thickness ratio (a/h) on normalized critical ui-axial and bi-axial buckling load () 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₁/E₂) on normalized critical ui-axial and bi-axial buckling load () 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.

3.8 Angle-ply (45⁰/-45⁰/-45⁰/45⁰) laminate having different thickness ratio (a/h)

In this section the effect of thickness ratio (a/h) on normalized critical pure shear buckling load (λ_cr = λ×10^-3) is studied taking E₁/E₂ =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 F^{igure 4}.

4 CONCLUSIONS

In the present study, the buckling behaviour of laminated composite plates is studied using an efficient C⁰ 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.

Received 10 Apr 2012

In revised form 23 Apr 2012

• [1] C.W. Bert and T.L.C. Chen. Effect of shear deformation on vibration of antisymmetric angle ply laminated rectangular plates. Int. J. Solids Struct, 14:465–473, 1978.
• [2] K. Bhaskar and T.K. Varadan. Refinement of higher order laminated plate theories. AIAA J, 27:1830–1831, 1989.
• [3] E. Carrera. Transverse normal stress effects in multilayered plates. J. Appl. Mech, 66:1004–1012, 1999.
• [4] A. Chakrabarti and A.H. Sheikh. Buckling of laminated composite plates by a new element based on higher order shear deformation theory. Mech. Composite Mater. Struct, 10(4):303–317, 2003.
• [5] A. Chakrabarti and A.H. Sheikh. A new triangular element to model inter-laminar shear stress continuous plate theory. Int. J. Numer. Meth. Engg, 60:1237–1257, 2004.
• [6] A. Chakrabarti and A.H. Sheikh. Buckling of laminated sandwich plates subjected to partial edge compression. Int. J.of Mech. sciences, 47:418–436, 2005.
• [7] M. Cho and R.R. Parmerter. An efficient higher order plate theory for laminated composites. Compos. Struct, 20:113–123, 1992.
• [8] M. Cho and R.R. Parmerter. Efficient higher order composite plate theory for general lamination configurations. AIAA J, 31(7):1299–1306, 1993.
• [9] M.E. Fares and A.M. Zenkour. Buckling and free vibration of non-homogeneous composite cross-ply laminated plates with various plate theories. Compos. Struct, 44:279–287, 1999.
• [10] A.J.M. Ferreira, L.M. Castro, C.M.C. Roque, J.N. Reddy, and S. Bertoluzza. Buckling analysis of laminated plates by wavelets. Comput. Struct, 89:626–630, 2011.
• [11] A.J.M. Ferreira, C.M.C. Roque, and P.A.L.S. Martins. Radial basis functions and higher order shear deformation theories in the analysis of laminated composite beams and plates. Compos. Struct, 66:287–293, 2004.
• [12] A.J.M. Ferreira, C.M.C. Roque, A.M.A. Neves, R.M.N. Jorge, C.M.M. Soares, and K.M. Liew. Buckling and vibration analysis of isotropic and laminated plates by radial basis functions. Composite Part B, 42:592–606, 2011.
• [13] L. Fiedler, W. Lacarbonara, and F. Vestroni. A generalized higher-order theory for buckling of thick multi-layered composite plates with normal and transverse shear strains. Compos. Struct, 92:3011–3019, 2010.
• [14] M. Ganapathi, Amit Kalyani, Bhaskar Mondal, and T. Prakash. Free vibration analysis of simply supported composite laminated panels. Compos. Struct, 90:100–103, 2009.
• [15] M. Ganapathi and D.P. Makhecha. Free vibration analysis of multi-layered composite laminates based on an accurate higher-order theory. Composites: Part B, 32:535–543, 2001.
• [16] T. Kant and K. Swaminathan. Free vibration of isotropic, orthotropic, and multilayer plates based on higher order refined theories. J. Sound Vib, 24:319–327, 2001.
• [17] A.A. Khdeir and L. Librescu. Analysis of symmetric cross-ply elastic plates using a higher-order theory. part ii: Buckling and free vibration. Compos. Struct, 9:259–277, 1988.
• [18] J. Lee. Thermally induced buckling of laminated composites by a layer-wise theory. Comput Struct, 65:917–922, 1997.
• [19] X.Y Li and Z. Liu. Generalized laminate theories based on double superposition hypothesis. Int. J. Numer. Methods Engg, 40:1197–1212, 1997.
• [20] D. Liu and X. Li. An overall view of laminate theories based on displacement hypothesis. J. Compos. Mater, 30:1539–1561, 1996.
• [21] L. Liu, L.P. Chua, and D.N. Ghista. Mesh-free radial basis function method for static, free vibration and buckling analysis of shear deformable composite laminates. Compos. Struct, 78:58–69, 2007.
• [22] Bhaskar Mondal, M. Ganapathi, and Amit Kalyani. On the elastic stability of simply supported anisotropic sandwich panels. Compos. Struct, 80:631–635, 2007.
• [23] H. Murakami. Laminated composite plate theory with improved in-plane responses. J. Appl. Mech, 53:661666, 196.
• [24] P. Nali, E. Carrera, and S. Lecca. Assessments of refined theories for buckling analysis of laminated plates. Compos. Struct, 93(2):456464, 2011.
• [25] H. Nguyen-Van, N. Mai-Duy, W. Karunasena, and T. Tran-Cong. Buckling and vibration analysis of laminated composite plate/shell structures via a smoothed flat shell element with in-plane rotations. Compos. Struct, 89:612–625, 2011.
• [26] A.K. Noor. Stability of multilayered composite plates. Fibre Sci. Technol, 8(2):81–89, 1975.
• [27] N.D. Phan and J.N. Reddy. Analyses of laminated composite plates using a higher-order shear deformation theory. Int. J. Numer. Methods Engg, 21:2201–2219, 1985.
• [28] J.N. Reddy. Free vibration of antisymmetric angle-ply laminated plates including transverse shear deformation by the finite element method. J. Sound Vib, 4:565–576, 1979.
• [29] J.N. Reddy. Cellular solids: structures & properties CRC Press, Boca Raton, FL, 2nd ed edition, 2004.
• [30] J.N Reddy and N.D. Phan. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher-order shear deformation theory. J. Sound Vib, 98:157–170, 1985.
• [31] J.N. Reddy and N.D. Phan. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher-order shear deformation theory. J. Sound Vibr, 98(2):157–170, 1985.
• [32] E. Reissner and Y. Stavasky. Bending and stretching of certain types of aelotropic elastic plates. J. Appl. Mech, 28:402–408, 1961.
• [33] D. Cook Robert, S. Malkus David, E. Michael Plesha, and J. Witt Robert. Concepts and Applications of Finite Element Analysis John Wiley & Sons Pte. Ltd, 4th edition edition, 2003.
• [34] M. Di Sciuva. A refined transverse shear deformation theory for multilayered anisotropic plates. Atti. Academia Scienze Torino, 118:279–295, 1984.
• [35] M. Di Scuiva. Multilayered anisotropic plate models with continuous interlaminar stress. Comput. Struct, 22(3):149167, 1992.
• [36] A.H. Sheikh and A. Chakrabarti. A new plate bending element based on higher order shear deformation theory for the analysis of composite plates. Finite Element Anal Des, 39:883–903, 2003.
• [37] I. Shufrin, O. Rabinovitch, and M. Eisenberger. Buckling of symmetrically laminated rectangular plates with general boundary conditions a semi analytical approach. Compos. Struct, 82:521–531, 2008.
• [38] P. Sundaresan, G. Singh, and G.V. Rao. Buckling of moderately thick rectangular composite plates subjected to partial edge compression. International Journal of Mechanical Sciences, 40(11):1105–1117, 1998.
• [39] W. Zhen and C. Wanji. Free vibration of laminated composite and sandwich plates using global-local higher-order theory. J. Sound Vib, 298:333–349, 2006.
• [40] W. Zhen and C. Wanji. Buckling analysis of angle-ply composite and sandwich plates by combination of geometric stiffness matrix. Comput. Mech, 39:839–848, 2007.

Publication Dates

• Publication in this collection
  17 Jan 2013
• Date of issue
  June 2012

History

• Received
  10 Apr 2012
• Reviewed
  23 Apr 2012