Buckling Analysis of Unidirectional Polymer Matrix Composite Plates Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 42 Buckling Analysis of Unidirectional Polymer Matrix Composite Plates Dr. Jawad Kadhim Uleiwi Material Engineering. Dept. / University of Technology (Received 14 November 2005; accepted 4 April 2006) Abstract:- This study deals with the estimation of critical load of unidirectional polymer matrix composite plates by using experimental and finite element techniques at different fiber angles and fiber volume fraction of the composite plate. Buckling analysis illustrated that the critical load decreases in nonlinear relationship with the increase of the fiber angle and that it increases with the increase of the fiber volume fraction. The results show that the maximum value of the critical load is (629.54 N/m) at ( = 0) and (Vf = 40 %) for the finite element method, while the minimum value of the critical load is (49 N/m) at ( = 90) and (Vf = 10 %) for the experimental results. The results also indicated that the maximum difference between the finite element analysis and experimental work is about (11 % ) at (  = 0) and (Vf = 40 %). 1.Introduction Many aeroengine components consist of thin-walled plates under loads which can potentially cause failure by buckling. These loads comprise combinations of compressive, shear and bending force [1]. The word "buckling" evokes an image of failure of structures which has been compressed in some way. Pictures and perhaps sound come to mind of sudden, catastrophic collapse involving very large deformations [2]. When the magnitude of the load on a structure is such that the equilibrium is changed from stable to unstable the load is called a critical load ( or buckling load ) [3]. The nature of the buckle pattern in plate depends not only upon the type of the applied loading but also upon the manner in which the edges are supported, shape (dimension) of the problem and material properties. The composite plate used in this study was a unidirectional polymer matrix composite plate made from glass fiber – epoxy matrix composite with constant length to width ratio (L/w = 2) at different fiber angles and fiber volume fraction. The aim of this work is to investigate the buckling analysis of the Al-khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol.2,No.2,pp 42-55 (2006) Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 43 composite plate by using the finite element analysis and the experimental work. An enormous amount of effort has been exerted on the subject of buckling since the beginning of the 20th century mainly the work performed by NACA (National Advisory Committee for Aeronautics). Rais-Rohani, M and Marcellier, P. [4,5] studied the effects of elastic edge conditions, aspect ratio, and face sheet ply pattern on the buckling loads and natural frequencies of anisotropic rectangular sandwich plates by using Rayleigh – Ritz method. Timothy L.C. Chen and Charles W. Bert [6] used the theory for buckling of laminated, anisotropic, thin plates to study buckling of the rectangular plates, laminated of composite material and subject to uniaxial compressive loading for optimal design plates laminated of glass/epoxy, boron/epoxy, and carbon/epoxy composite materials. Hyonny Kim and Keith T. Kedward [7] studied the global buckling and local sublaminate buckling of the laminated plate containing Zone of delamination by using Rayleigh – Ritz method and finite element models. Timoshenko, Stephen P. and James M. Gere [8] explained the approximate method to obtain the results of buckling stresses which are always higher than the exact solution. Yasui [9] presented what appears to be the first in-depth parametric study of buckling behavior of laminated composite plates with a central circular cutout that has been obtained by using the finite element method. Rouse [10] presented an experimental investigation of the buckling and postbuckling behavior of square graphite-epoxy and graphite- thermoplastic plates loaded in shear. 2.Theoretical Analysis Using lamination theory, the strains in structural laminate resulting from a given applied stress from information on the basic lamina which are: glass and epoxy moduli, poisson's ratios, densities, and volume fractions of fibers, in addition to geometric data for the composite plate. The constitute relations may be established for a lamina by assuming that the lamina is a homogeneous orthotropic material in a state of plane stress. The geometry of a single filamentary lamina is shown in figure (1) [11]. Hook's law for a composite material with orthotropic properties in a plane stress state is:                                 12 2 1 33 2221 1211 12 2 1 00 0 0       Q QQ QQ (1) The components of the reduced stiffness matrix, [Q]. are:   21121111 1/  EQ (2)   211211212112 1/   EQQ (3) 1233 GQ  (4) And the longitudinal and lateral elastic modulus is predicted by the rule of mixtures formula. mmff VEVEE  1 (5) fmmf mf VEVE EE E    2 (6) And the major poisson's ratio is given by:- mmff VV   12 (7) And the minor poisson's ratio 21 can be estimated from 1 2 1221 E E   (8) Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 44 For a given fiber angle ( ) the stress/strain relationship                                 xy y x xy y x QQQ QQQ QQQ       333231 232221 131211 (9)     4 22 22 3312 4 1111 sincossin 22cos   Q QQQQ     4 22 22 3312 4 1122 coscossin 22sin   Q QQQQ     44 12 22 3322112112 cossincossin 4   Q QQQQQ   )cos(sincossin 22 44 33 22 3312221133    Q QQQQQ       cossin2cos sin2 3 332212 3 3312113113   QQQ QQQQQ       3 332212 3 3312113223 cossin2cos sin2   QQQ QQQQQ Also the stress/strain equation in the ply can be written in reference to the geometric midplane as follows:                                                        xy y x xy y x xy y x k k k QQQ QQQ QQQ Z QQQ QQQ QQQ 333231 232221 131211 333231 232221 131211          (10) To find the force resultants acting on the total laminate, Equation (10) is integrated over the total thickness in the z direction.                                        n k z z xy y x k t t xy y x xy y x dzdz N N N k k 1 2/ 2/ 1       (11)                                        n k z z xy y x k t t xy y x xy y x dzzdzz M M M k k 1 2/ 2/ 1       (12) After integration and rearranging equations (11 and 12 ), the constitutive equations can be obtained as follows:                                                                     xy y x xy y x xy y x xy y x k k k DDDBBB DDDBBB DDDBBB BBBAAA BBBAAA BBBAAA M M M N N N       333231333231 232221232221 131211131211 333231333231 232221232221 131211131211 (13) Where:         n k kkijij zzQA 1 1         n k kkijij zzQB 1 2 1 2 2 1         n k kkijij zzQD 1 3 1 3 3 1 3.Experimental Work The experimental part was carried out in the field to determine experimentally the critical load on many specimens of the composite material of unidirectional glass fiber – epoxy matrix composite plates with different fiber angles and volume fraction of fiber as following: Vf = 10 %, 20%, 30% and 40%  = 0, 15, 30, 45, 60, 75, 90 The dimension of the test specimen has a length of (200 mm) and width of (100 mm) and a thickness of (1mm). Figure (2) represents the test machine with a test specimen of the Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 45 composite plate with some accessories like dial-gauge and load indicator to measure the lateral deflection of the composite plate. 4.Finite element Analysis The finite element analysis carried out as part of this work was performed using the ANSYS package in the buckling analysis of composite plate to determine the critical load at which the structure becomes unstable. Element Selected This element is used to modal the composite plate. The element has six degrees of freedom at each node translation in the x,y and z directions and rotations about the nodal x,y and z-axes. The geometry, node locations, and the coordinate system for this element are shown in figure (3). The element is defined by eight nodes four thicknesses and orthotropic material properties [12]. 5.Boundary Conditions And Mesh Generation There are two main types of boundary condition (displacement and load boundary condition). The composite plate is a simply – supported plate from the two opposite edge and free from other edges and the load is applied axially on the two opposite edges. As for the mesh generation of the composite plate, it is necessary to discretize it into a sufficient number of elements in order to obtain a reasonable accuracy. Figure (4) represents the mesh generation and boundary conditions of the composite plate. 6.Results and Discussions The results obtained from the experimental work and the finite element analyses of the buckling analysis of the composite plate are discussed here. Figure (5) shows the lateral deflection contours and the value of the critical load of the unidirectional composite plate at (Vf = 10 % ) It is clear from these figures that the buckling happens at the midline of the plane and the values of (Pcr = 195.83 N/m, 139.82 N/m, 73.111 N/m and 55.142 N/m ) at fiber angle are (  = 0, 30, 60 and 90) respectively. Figure (6) shows the relationship between the critical load and the fiber angle of the composite plate at different fiber volume fraction and at ( L/w = 2)for the finite element analysis. It is clear from this figure that the critical load decreases in nonlinear relationship with the increase of the fiber angle in different manners depending on the volume fraction. The maximum value of ( Pcr = 629.54 N/m ) at ( = 0 and Vf = 40 %) while the minimum value of ( Pcr = 55.142 N/m ) at ( = 90 and Vf = 10 %). Figure (7) shows the relationship between the critical load and fiber volume fraction of the composite plate at different fiber angles for the finite element method. It can be seen that the critical load increases with the increase of the fiber volume fraction in different manners depending on the fiber angle, where the behavior of curve at ( = 0 ) is the greatest one. At this angle the maximum value of (Pcr = 629.54 N/m) at (Vf = 40 %) while the minimum value of (Pcr = 195.83 N/m) at (Vf = 10 %). Figure (8) shows the relationship between the critical load and the fiber angle of the composite plate at different fiber volume fractions and at ( L/w = 2)for the experimental work. It is clear from this figure that the critical load decreases in nonlinear relationship with the increase of the fiber angle in different manners depending on the volume fraction. The maximum value of ( Pcr = 560 N/m ) at ( = 0 and Vf = 40 %) while the minimum value of Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 46 ( Pcr = 49 N/m ) at ( = 90 and Vf = 10 %). Figure (9) shows the relationship between the critical load and the fiber volume fraction of the composite plate at different fiber angle for experimental work. It can be seen that the critical load increases with the increase of the fiber volume fraction in different manners depending on the fiber angle, where the behavior of the curve at ( = 0 ) is the greatest one. At this angle the maximum value of (Pcr = 560 N/m) at (Vf = 40 %) while the minimum value of (Pcr = 173 N/m) at (Vf = 10 %). Figure (10) shows the relationship between the critical load and the fiber angle at ( Vf = 40 %) and at (L/w = 2) for the experimental and finite element analysis. It is clear from this figure that the critical load decreases in nonlinear relationship with the increase of the fiber angle. It has also been found that the maximum difference between the finite element analysis results and the experimental results ( 11 %) at ( = 0). The reason for this difference is attributed to the condition of preparing the specimens. Figure (11) shows the relationship between the critical load and the fiber volume fraction at ( = 0) and at (L/w = 2) for the experimental and finite element analysis. It is clear from this figure that the critical load increases in linear relationship with the increase of the fiber volume fraction and it has also been found that the maximum difference between finite element analysis results and experimental results ( = 11 %) at ( Vf = 40 %). 7.Conclusion The main conclusions of the buckling analysis of the composite plates for finite element analysis and the experimental work are:- (1-) Critical load decreases in nonlinear relationship with the increase of the fiber angle in different manners depending on the fiber volume fraction and the maximum value of (Pcr = 629.54 N/m) at ( = 0) and (Vf = 40 %) for finite element analysis. (2-) Critical load increases with the increase of the fiber volume fraction and the minimum value ( Pcr = 55.142 N/m ) at ( = 90) and (Vf = 10 %) for finite element analysis. (3-) The behavior of the critical load with fiber angle and with fiber volume fraction for experimental work is similar to the behavior of finite element analysis but here the maximum value of (Pcr = 560 N/m) at ( = 0) and (Vf = 40 %) the minimum value of ( Pcr = 49 N/m ) at ( = 90) and (Vf = 10 %). (4-) Maximum difference between the results of finite element analysis and experimental work results ( = 11 %) at ( = 0) and ( Vf = 40 %). 8.References [1-] C.A. Featherston and C Ruiz, " Buckling of Flat Plates under Bending and Shear ", Proc Instn Mech Engrs, Vol.212, Part C, (1998). [2-] David Bushnell, " Buckling of Shells - Pitfall for Designers", AIAA80-0665R, Vol.19, No.9, September, (1981). [3-] N.G. Lyenger and S.Gupta, " Programming Methods in Structure Design", (1980). Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 47 [4-] Rais-Rohani, M and Marcellier, P.," Buckling and Vibration Analysis of Composite Sandwich Plates with Elastic Rotational Edge Restraints". AIAA Journal, Vol.37, No. 5, (1999). [5-] Marcellier, P. and Rais-Rohani, M., " Free Vibration and Buckling Analysis of Anisotropic Sandwich Plates with Edges Elastically Restrained Against Rotation", Proceedings of the 39th AIAA/ ASME/ ASCE/ AHS/ ASC Structures, Structural Dynamics, and Materials Conference, Long Beach, CA, April 20-23, 1998. AIAA-98-2084-CP. [6-] Timothy L.C. Chen and Charles W. Bert " Design of Composite Material Plates for Maximum Uniaxial Compressive Buckling Loads", Mechanical and Nuclear Engineering, The University of Oklahoma, Norman, Oklahoma, (1976). [7-] Hyonny Kim and Keith T. Kedward," A Method for Modeling the Local and Global Buckling of Delaminated Composite Plate ", Composite Structures, Vol.44, No.1, (1999). [8-] Timoshenko, Stephen P. and James M. Gere, " Theory of Elastic Stability", 2nd Edition, McGraw- Hill Book Company, New York, (1961). [9-] Yasui, Yoshiaki," The Buckling of Rectangular Composite Plates with Cutout under Uniaxial and Biaxial Compression", Proceedings the 8th International Conference on Composite Materials, July, (1991). [10-] Rouse, Marshall," Effect of Cutouts or Low-Speed Impact Damage on the Postbuckling Behavior of Composite Plates Loaded in Shear", Proceedings of the 31st AIAA, ASME, ASCE, AHS, ASC, Structures, Structural Dynamics and Materials Conference, April, (1990). [11-] R.M., Jones, " Mechanics of Composite Material", McGraw- Hill, New York, (1975). [12-] Kohnke P.," ANSYS Theory Reference Release 5.4", ANSYS, Inc., (1975). Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 48 Aij Extensional stiffness terms (N/m). Bij Coupling stiffness terms (N). Dij Bending stiffness terms (N.m). E11 Young's modulus in the fiber direction (N/m2). E22 Young's modulus in the matrix direction (N/m2). Ef Young's modulus of the fiber (N/m2). Em Young's modulus of the matrix (N/m2). G12 Shear Modulus along the principal axis (N/m2). kx, ky, & kxy Curvatures of the midplane (m). Mx Resultant moment intensity about the y-axis per unit width of the laminate cross- section (N.m). My Resultant moment intensity about the x-axis per unit width of the laminate cross- section (N.m). Mxy Resultant twisting moment intensity about the xy- axis/unit width of the laminate cross-section (N.m). Nx Resultant in-plane force intensity in the x-axis/unit width of the laminate cross- section (N). Ny Resultant in-plane force intensity in the y-axis/unit width of the laminate cross- section (N). Nxy Resultant in-plane shearing force intensity in the xy- axis/unit width of the laminate cross-section (N). ij Q Reduced stiffness terms (N/m2). ij Q Transformed reduced stiffness terms (N/m2).  Q Reduced stiffness matrix.  Q Transformed reduced stiffness matrix. Vf Volume fraction of fiber (%). Vm Volume fraction of matrix (%). z Coordinate measured from the midplane (m).  xyyx  ,, Midplane strains.  Strain vector in ply.  Stress vector in ply. i Stress components (N/m 2). 12 Shear strength along the principal axis (N/m2) 12 Shear strain along the principal axis.  Fiber angle (). 12,21 Major and Minor poisson's ratio. Notation Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 49 Figure (1): Location of Material (Principal) Coordinate System on the Lamina [11]. Figure ( 2 ): Test Apparatus and Test Specimen Load Indicator Dial Gauge Unidirectional Composite Plate Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 50 Figure ( 3 ): SHELL93 8-Node Structural Shell [12]. Figure ( 4 ): Mesh Generation and Boundary Conditions of the Composite Plate. Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 51 Figure ( 5 ): Deflection Contours of the Composite Plates at (Vf=10 %) and at Different Value of Fiber Angle.  = 0  = 30  = 60  = 90 Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 52 Figure ( 6 ): Relationship Between Critical Load and Fiber Angle for Finite Element Analysis at Different Volume Fraction. Figure ( 7 ): Relationship Between Critical Load and Fiber Volume Fraction for Finite Element Analysis at Different Fiber Angle. 10 20 30 40 Fiber Volume Fraction (%) 0 50 100 150 200 250 300 350 400 450 500 550 600 650 C ri ti c a l L o a d ( N /m ) F.E.M.  = 0  = 45  = 90 0 15 30 45 60 75 90 Fiber Angle (Deg.) 0 100 200 300 400 500 600 700 C ri ti c a l L o a d ( N /m ) F.E.M. Vf = 40 % Vf = 30 % Vf = 20 % Vf = 10 % Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 53 Figure ( 8 ): Relationship Between Critical Load and Fiber Angle for Experimental at Different Volume Fraction. Figure ( 9 ): Relationship Between Critical Load and Fiber Volume Fraction for Experimental at Different Fiber Angle. 0 15 30 45 60 75 90 Fiber Angle (Deg.) 0 100 200 300 400 500 600 C ri ti c a l L o a d ( N /m ) Experimental Vf = 40 % Vf = 30 % Vf = 20 % Vf = 10 % 10 20 30 40 Fiber Volume Fraction (%) 0 50 100 150 200 250 300 350 400 450 500 550 600 C ri ti c a l L o a d ( N /m ) Experimental  = 0  = 45  = 90 Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 54 Figure ( 10 ): Relationship Between Critical Load and Fiber Angle for Finite Element Analysis Experimental at Fiber Volume Fraction( = 40 %). Figure (11): Relationship Between Critical Load and Fiber Volume Fraction for Finite Element Analysis and Experimental at Fiber Angle (= 0). 0 15 30 45 60 75 90 Fiber Angle (Deg.) 0 100 200 300 400 500 600 700 C ri ti c a l L o a d ( N /m ) F.E.M Experimental 10 20 30 40 Fiber Volume Fraction (%) 0 50 100 150 200 250 300 350 400 450 500 550 600 650 C ri ti c a l L o a d ( N /m ) F.E.M. Experimental Dr. Jawad Kadhim Uleiwi /Al-khwarizmi Engineering Journal ,Vol.2, No. 2 PP 42-55 (2006) 55 تحليل األنبعاج للصفائح المركبة األحادية األتجاه ذات األساس البوليمري د. جواد كاظم عليوي قسم هندسة المواد / الجامعة التكنولوجية :الخالصة تهتم هذه الدراسة بتقيم الحمل الحرج للصفائح المركبة األحادية األتجاه ذات األسااس الباوليمري بتساتخدام طريقااة تقنيااات العناصاار المحااددو انااد ماويااا ليااح مختلفااة وكساار حجماا مختلااح لتلاا الطريقااة العمليااة و الصفائح. اذ يبين تحليل األنبعاج بان الحمل الحرج يقل بعالقاة ييار خطياة ماا مياادو ماوياة اللياح ويامداد ماا مياادو الكسر الحجم لليح. ( = 0 ( اناد N/m rcP 629.5 =بينت النتائج بان اقصى قيمة للحمل الحرج لكل وحدو طاول ( اناد N/m crP 49 =( لطريقة العناصر المحددو. بينما اقل قيمة للحمل الحرج الحارج fV % 40 =و = 90  10 =( و % fV.للطريقة العملية ) = ( اند % 11 = وكذل بينت النتائج بان اقصى فرق بين طريقة العناصر المحددو والطريقة العملية 0 40 =( و % fV .)