3D Buckling Analysis of Multidelaminated Composite Specimens

©Freund Publishing House Ltd. Science and Engineering of Composite Materials 14, 181-188 (2007) 3D Buckling Analysis of Multidelaminated Composite Specimens D. Tumino1, F. Cappello, D. Rocco Dipartimento di Meccanica, Università degli Studi di Palermo, 90128, Palermo, Italy. ABSTRACT The behaviour of thin composite laminates (unidirectional, cross-ply and angle-ply) under compressive loads has been examined in cases where multiple delaminations are present. The problem is solved using the Finite Element Method (FEM) both with linear analyses, based on the eigenvalues research problem, and with nonlinear analyses, based on incremental-iterative procedures. In particular, the role of the delamination length, of the angle of the plies and of the stacking sequence on the critical load is investigated. Results are compared with those found in literature derived from experimental or numerical 2D analyses. Key Words: composite delaminations, FEM analysis. laminates, buckling, 1. INTRODUCTION The use of laminated composite materials is, nowadays, a well-established practice in several engineering fields, particularly in aeronautics where the high specific stiffness and strength of such materials allows consistent reduction in weight for the finished structures if compared with the ones made with traditional materials. Long-fibre reinforced plastics are largely used and have been studied in details /1/. Due to imperfections caused by the manufacturing processes or to slow-speed impacts or to stress concentrations (like 1 those in correspondence of the free edges) these laminated materials can exhibit partial or total separation between the laminae and a resulting loss of stiffness and strength. The buckling phenomenon (the mechanical instability of a thin structure subjected to compressive loads) can lead to out-of-service of the structure although the level of the internal stress is not critical. This is the reason why a correct design of thin components must take into account the instability problems and, when dealing with composite materials, with the presence of delaminations. A plane composite with delaminations subjected to a uniaxial compressive load can exhibit a global, a local or a mixed buckling (Fig. 1) depending on the position and on the length of delaminations. When global buckling occurs the whole structure is deformed (Fig. 1.a) but there is no evident separation between each delaminated part (sublaminate). In the case of local buckling, only one (or few) sublaminate shows large deflections while the rest of the structure remains almost undeformed (Fig. 1.b). In the mixed case the whole laminate buckles and each sublaminate can undergo different curvatures (Fig. 1.c,1.d). When the deflection of the sublaminate subjected to local buckling is free to occur we talk about unconstrained buckling, while constrained buckling occurs when the local buckling of one (or more) sublaminates is prevented by adjacent sublaminates. The aim of the buckling analysis is to determine the critical load in correspondence with the neutral equilibrium between external and internal forces under Phone: +39916657124, Fax: +3991484334, Email: tumino@dima.unipa.it Vol. 14, No. 3, 2007 large deflection conditions, and the related deformed shape of the structure. Fig. 1: Mode of buckling: (a) global, (b) local, (c,d) mixed. Some authors use analytical solutions for this problem. Shu /2/ solved the case of a beam with two centred delaminations; sublaminates are considered as beams linked to each other. Chen /3/ has solved the case with a single delamination; sublaminates are treated as plates and the solution is achieved using energetic methods; the possibility for the crack to grow is considered by evaluating the strain energy release rate /4/. All analytical approaches apply to single or double delaminated specimens and in the case of unidirectional laminates. When dealing with a greater number of delaminations or with different stacking sequences, alternative methods must be used, as experiments or numerical analyses. With FEM the buckling problem can be solved adopting a linear and a nonlinear approach. Linear buckling furnishes information when the structure has reached the neutral equilibrium between external loads and elastic reactions (the bifurcation point), under the hypotheses that the load-displacements diagram is linear as load increases and large displacements occur just at the moment of buckling. In /5/ linear buckling analyses of multi-delaminated composite plates are performed, using shell elements with Mindlin formulation; interpenetration between sublaminates is avoided using penalty functions. Non-linear buckling with FEM allows one to consider the presence of large deflections and contact constraint between sublaminates; this means that not only the bifurcation point but also pre- and postbuckling phases can be studied without any restrictive hypothesis. In /6/ and /7/ the case of a single 182 3D Buckling Analysis of Multidelaminated Composite Specimens delamination with different lengths is analysed in unidirectional and cross-ply laminates either as throughthe-width and as embedded circular flaws, modelling the structure with shell elements. Shell elements are also used in /8/ for cases of multi-delaminations in unidirectional laminates. In /9/ and /10/ composite laminates, with fibres oriented along the load direction, are studied in the presence of multiple delaminations variable in number and length using a model meshed with 2D plane elements; a comparison between linear and non-linear method has also been done. In a previous work /11/ the effects of number, position and the length of delaminations and stacking sequence in unidirectional and cross-ply specimens were examined by performing 2D plane FE analyses. In the present work angle-ply laminates are also considered with the aim to understand the influence of the misalignment of laminae on the critical load of a multidelaminated specimen. The use of 3D models is required by such problems. Results obtained both with linear and nonlinear numerical approaches are validated with analytical and numerical results from the literature. 2. NUMERICAL PROCEDURE 2.1. Linear and nonlinear approach A brief explanation on the numerical approach used in the analyses is summarized here; more details can be found in /11/. The linear approach gives information about the structure at the time corresponding to the neutral equilibrium between external forces and internal reactions /11//12/ by solving the equation: > K @ ^d ` > K0 @  O > KV @ ^d ` ^0` (1) where > K @ is the global stiffness matrix, > K 0 @ is the stiffness matrix under small displacements, > KV @ is the geometric stiffness matrix correspondent to a reference load, O is the eigenvalue (the factor to be multiplied to the reference load to obtain the critical value) and ^d ` is the vector of displacements. Assuming that ^d ` is not equal to the zero vector ^0` , solutions of equation D. Tumino, F. Cappello and D. Rocco Science and Engineering of Composite Material (1) can be found for different values of O by solving the following expression: det > K 0 @  O > KV @ 0 (2) Equation (2) is an nth-order polynomial in the unknown O that admits n different roots. Due to the fact that we research the smallest value of the load that causes instability, we are interested only in the smallest value of O, together with the related eigenvector ^d ` that represents the deformed shape of the structure at the instant of buckling. The nonlinear approach consists in performing a standard incremental-iterative analysis. With this approach several nonlinearities that characterise the problem can be modelled: for our case, for example, large displacements and contact between adjacent sublaminates. The entire evolution of the phenomenon can be observed by adopting a nonlinear approach, not only the moment of buckling but also the pre- and postbuckling phases. By monitoring the internal reactions of the structure as a whole and of each portion of it, it is possible to evaluate the critical load as the one where the tangent to the load/displacement curve becomes horizontal. The linear approach has a great advantage, compared with the nonlinear approach, that solution can be achieved in a very short computational time. On the contrary, no information can be carried out regarding pre and post-buckling and, in particular cases, interpenetration can occur between adjacent laminae. As explained later, this can cause an underestimation of the load-carrying capability of the buckled specimen. 2.2. Analysed model and procedure setup Assume that the specimen is clamped at the two ends and that delaminations are through the width and centred with respect to the length (see Fig. 2). It is also assumed that the length of the delaminations remains constant during the analysis. Geometric and elastic properties are taken from /13/. In particular, the entities that do not vary are: height H=1.776 mm, length L=90 mm, width B=18 mm. Specimens are made of n=16 laminae of a graphite/epoxy material, which elastic properties are reported in Table 1. To analyse the problem, 3D FE models have been realised. The mesh is obtained using hexahedral structural 8-noded elements from the element library of the code ANSYS®. To avoid interpenetration between adjacent sublaminates during the buckling process, surface-to-surface contact elements have been added in the delaminated areas. The total number of elements is a compromise between accuracy of results and computational effort: preliminary analyses have been performed in order to search for an optimal value for the element dimension. Figure 3 shows the variation of the eigenvalue O with the total number of nodes of the structure for a unidirectional undelaminated specimen. From this graph one can observe that an approximate number of 5·104 nodes gives satisfactory results. The model is such that any stacking sequence can be simulated and the angle of misalignment of the laminae can vary from 0° to 90°. Delaminations can be inserted between each pair of adjacent laminae. Fig. 2: Analysed model. Table 1 Elastic properties of the graphite/epoxy material, from /13/. E11 130 GPa E22 = E33 10 GPa Q12 = Q13 0.31 Q23 0.5 G12 = G13 4.85 GPa G23 3.62 GPa 183 Vol. 14, No. 3, 2007 3D Buckling Analysis of Multidelaminated Composite Specimens O 34 33.5 33 4000 24000 44000 64000 84000 number of nodes Fig. 3: Influence of the number of nodes on the solution. 3. NUMERICAL RESULTS 3.1. Comparison between 2D and 3D analyses A preliminary activity has been accomplished to compare numerical results given by the 3D model with those obtained in /11/ using a 2D model. Due to the fact that for unidirectional specimen a very good agreement was demonstrated between 2D numerical results and analytical results by Shu /2/, attention has been now focused mainly on the behaviour of cross-ply specimen as [04//904]s with double delamination (the symbol “//” identifies the position of the delaminations in the laminate). As will be repeated in the following, we assume that delaminations are located at the interface between sublaminates with a different angle of misalignment. Numerical 3D analyses of buckling on the [04//904]s specimen lead, for the undelaminated material, to a buckling load Pc=4550 N, calculated both with linear and nonlinear approach. The same results was found in /11/ with 2D analyses. Figure 4 shows the variation with the relative length a/L of the normalised buckling load P/Pc for the specimen with two equal-length delamination. In this diagram three different zones can be distinguished: the first for 0 d a / L  0.166 where only global buckling appears, the second for 0.166  a / L  0.3 which characterised the mixed buckling, and then the third zone for 0.3 d a / L  1 184 where only local buckling occurs. The type of buckling that occurs can be recognized exactly by analysing the curvatures of each sublaminate and of the overall structure at the moment of buckling. One can observe that the main loss of load capability due to the instability appears in the mixed buckling phase. There is a very good agreement between results reported in Fig. 4 relative to 3D analyses. The presence of delamination at the interface between laminae with different angles has a strong influence on the buckling behaviour of the specimen, also for very small value of a/L. As the delamination length increases, the critical load decreases to negligible value (for a/L=1 we calculate P/Pc=4%). The same agreement is appreciable between results of 2D analyses from /11/ and results from 3D analyses in Table 2. Table 2 Normalised buckling load for the [04//904]s laminate. a/L 2D - 3D - linear 2D 3D - Nonlinear 0 1.000 1.000 1.000 1.000 0.1 0.998 0.990 1.000 0.990 0.167 0.923 0.925 0.900 0.943 0.2 0.769 0.759 0.777 0.776 0.3 0.376 0.379 0.371 0.374 0.4 0.219 0.218 0.221 0.218 0.5 0.143 0.144 0.144 0.144 0.7 0.075 0.075 0.072 0.069 1 0.038 0.040 0.043 0.040 1.2 1 0.8 3D - Linear P/Pc 34.5 3D - Nonlinear 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 a/L Fig. 4: Normalised buckling load versus relative crack length for the [04//904]s laminate. D. Tumino, F. Cappello and D. Rocco Science and Engineering of Composite Materials 3.2. 3D simulations for angle-ply specimens (a) 5000 4500 P [N] Three types of angle-ply laminates have been studied using 3D models: x AP1: Symmetrical laminates [04//T4]s with two equal-length delaminations (h1/H=0.25, h2/H=0.75) for T=15°, 30°, 45°, 90°. x AP2: Symmetrical balanced laminates [04//+T2//-T2]s with four equal-length delaminations (h1/H=0.25, h2/H=0.375, h3/H=0.625, h4/H=0.75) for T=30°, 45°, 60°, 90°. x AP3: Symmetrical laminates [02//+T2//-T2//902]s with six equal-length delaminations (h1/H=0.125, h2/H=0.25, h3/H=0.375, h4/H=0.625, h5/H=0.75, h6/H=0.875) for T=30°, 45°. At first numerical analyses have been performed to calculate buckling loads of undelaminated specimens for the three configurations stated above. Linear and nonlinear analyses give very similar results, therefore, for reason of clarity, only results from linear analyses are reported. Figure 5 shows the absolute buckling load for the three configurations versus the angle of misalignment of the laminae from 0° to 90°. For configurations AP1 and AP2 we notice that the critical load decreases with T very slowly up to T =45°; then for 45°<T <90° the load is constant. Moreover, there is no preference between AP1 and AP2 because absolute values of the load are about the same. The configuration AP3 exhibits smaller values of loads, with respect to AP1 and AP2 (as expected because of a lower number of 0° laminae on the surface), and a great influence of T. 5500 AP1 AP2 4000 AP3 3500 3000 0 15 30 45 60 75 90 T>q@ Fig. 5: Absolute buckling load versus angle of misalignment. As previously stated, linear analyses can lead to incompatibilities during the buckling phenomenon, in particular for high values of a/L. Figure 6 shows the two typical shapes of buckled specimens resulting from linear analyses for the case [04//454]s with a/L=0.3. In Fig. 6.a the map of displacements in the z-direction is shown for the S-shaped buckling mode, and in figure 6.b for the thin-film buckling mode. Both of deformed configurations exhibit interpenetration between sublaminates and are related to different eigenvalues. It is not easy to predict which one of the buckling modes is associated to the smallest eigenvalue. On the other hand, results from nonlinear analyses give, in most cases, the thin-film mode as the buckled shape associated to the smallest eigenvalue. (b) Fig. 6: Shape of the buckled specimen: (a) S-shaped, (b) thin-film. 185 Vol. 14, No. 3, 2007 3D Buckling Analysis of Multidelaminated Composite Specimens 1.2 1 0.8 P/Pc The diagrams of normalised buckling load versus relative crack length for the configuration [04//T4]s are plotted in Figures 7 and 8 for different values of T. As can be clearly noticed, there is no relevant difference when comparing results from linear and nonlinear approach. Results for T=30, 45 and 60 are similar, only the curve for T=15° shows higher values of normalised buckling load. For this configuration global buckling occurs for a/L<0.133, mixed buckling for 0.133<a/L<0.3 and local buckling for a/L>0.3. 30° 45° 0.6 60° 90° 0.4 0.2 0 1.2 0 0.2 0.4 0.6 0.8 1 a/L 1 P/Pc 0.8 Fig. 9: Normalised buckling load versus relative crack length for the [04//+T2//-T2]s laminate: results from linear analyses. 15° 30° 0.6 45° 90° 0.4 1.2 0.2 1 0 0 0.2 0.4 0.6 0.8 1 0.8 Fig. 7: Normalised buckling load versus relative crack length for the [04//T4]s laminate: results from linear analyses. P/Pc a/L 30° 45° 0.6 60° 90° 0.4 0.2 1.2 0 0 1 0.2 0.4 0.6 0.8 1 a/L 0.8 P/Pc Fig. 10: 15° 30° 0.6 45° 90° 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 a/L Fig. 8: Normalised buckling load versus relative crack length for the [04//T4]s laminate: results from nonlinear analyses. 186 Normalised buckling load versus relative crack length for the [04//+T2//-T2]s laminate: results from nonlinear analyses. The diagrams of normalised buckling load versus relative crack length for the configuration [04//+T2//-T2]s are plotted in Figures 9 and 10 for different values of T. In this case the differences between linear and nonlinear approach are not negligible. Table 3 summarises the numerical values obtained with the two approaches. For different values of T we notice, from figure 9, that the curves from the linear approach are different; in particular as T increases the load decreases for the same D. Tumino, F. Cappello and D. Rocco Science and Engineering of Composite Materials 1.2 1 0.8 P/Pc value of a/L. This behaviour is in contradiction with the one shown by Figure 10 relative to nonlinear analyses. In this case there is no influence of T on the buckling response as the coincidence of the four curves confirms. The reason of this discrepancy between the two approaches can be motivated because of the incompatibility of the deformed structure led by the linear approach. In fact sublaminates at +T are those that buckle first, and their deflection is free when the structure is analysed with the linear approach. When nonlinear analyses are performed, deflection of sublaminates at +T is prevented by contact elements at the interfaces with external sublaminates at 0°. Then higher levels of load are needed to cause deflection of external and internal sublaminates. 30° 0.6 45° 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 a/L Fig. 11: Normalised buckling load versus relative crack length for the [02//+T2//-T2//902]s laminate: results from nonlinear analyses. Table 3 Normalised buckling load for the [04//+T2//-T2]s laminate. T = 30° a/L Linear 0 1.0000 0.133 0.9837 0.167 0.9352 0.2 0.7643 0.3 0.3153 0.4 0.1622 0.5 0.0944 0.7 0.0513 1 0.0182 Nonlin. 1.0000 0.9802 0.9282 0.7015 0.3305 0.1954 0.1300 0.0668 0.0359 a/L 0 0.1 0.167 0.2 0.3 0.4 0.5 0.7 1 T= 45° Linear 1.0000 0.9959 0.9115 0.6032 0.2314 0.1148 0.0666 0.0298 0.0139 Nonlin. 1.0000 0.9889 0.9411 0.7380 0.3358 0.2051 0.1344 0.0706 0.0366 The meaningless physical behaviour of linear results, for this configuration, allows us to state that only Fig. 10 can be considered reliable. For AP2 configuration global buckling occurs for a/L<0.167, mixed buckling for 0.167<a/L<0.5 and local buckling for a/L>0.5. For the laminates AP3 with 6 delaminations, the normalised buckling load with the relative crack length has been studied for T=30° and T=45°. Figure 11 shows that small differences existing between the two cases: the values of the critical loads for T =45° are a little bit smaller than those for T=30°. Global buckling occurs for a/L<0.067, mixed buckling for 0.067<a/L<0.3. Results are obtained with nonlinear approach only, because a/L 0 0.1 0.133 0.2 0.3 0.4 0.5 0.7 1 T = 60° Linear 1.0000 0.9987 0.9375 0.3666 0.1467 0.0776 0.0477 0.0270 0.0116 Nonlin. 1.0000 0.9893 0.9771 0.7726 0.3621 0.2196 0.1365 0.0703 0.0350 T = 90° a/L Linear 0 1.0000 0.1 0.9967 0.133 0.6297 0.2 0.2618 0.3 0.1129 0.4 0.0628 0.5 0.0400 0.7 0.0203 1 0.0106 Nonlin. 1.0000 0.9669 0.9594 0.7271 0.3475 0.2057 0.1191 0.0694 0.0350 linear analyses lead to interpenetration between sublaminates and a resulting underestimation of the critical loads. With respect to AP1 and AP2 configurations, we notice a significant influence of the delamination length on the relative buckling load confirmed by the small threshold value between global and mixed buckling. CONCLUSIONS A numerical procedure has been set up using 3D models within a FEM code to study the buckling phenomenon in composite specimens with multiple 187 3D Buckling Analysis of Multidelaminated Composite Specimens Vol. 14, No. 3, 2007 delaminations. The linear and nonlinear numerical approaches are compared. This procedure has been adopted to analyse different types of symmetric configurations: non-orthotropic as the [04//T4] laminate, balanced as the [04//+T2//-T2] laminate and the [02//+T2//T2//902] laminate. The assumption is that the delaminations are located at the interfaces between sublaminates with different angles of misalignment. The influence of the relative length of the delaminations and of the orientation of the laminae on the critical buckling load has been studied. The threshold values between global and mixed buckling and between mixed and local buckling have been evaluated. Results obtained with the linear and nonlinear analyses are in good agreement for the cases where sublaminates that exhibit instability at first are located in the surfaces of the laminate (unconstrained buckling); in the other cases (constrained buckling) the linear analyses lead to interpenetration between sublaminates and, most important, to the underestimation of the critical load. For cross-ply laminates 2D and 3D models give the same results; therefore, due to the advantage in terms of computational effort saved, 2D models are preferred. For all the configurations examined, 3D analyses confirmed that the absolute buckling load decreases with the angle of misalignment between the fibres and the compressive load. When a single configuration is considered, the relative buckling load is not influenced in a significant way by the orientation of the fibres. In the global and in the local buckling zones very small variations of the relative buckling load with the relative crack length can be observed, while in the mixed buckling zone we can observe a remarkable reduction of the buckling load with the relative crack length. This reduction is as much appreciable as the delaminations are located near the surfaces of the specimens. REFERENCES 1. Agarwal BD, Broutman L-J. Analysis and Performance of Fiber Composites. John Wiley & Sons, New York, 1980. 188 2. Shu D. Buckling of multiple delaminated beams. Int. J. Solids Structures 1998; 35(13): 1451-1465. 3. Chen HP. Shear deformation theory for compressive delamination buckling and growth. AIAA Journal 1991; 5: 813-819. 4. Anderson TL. Fracture Mechanics: Fundamentals and Applications. CRC Press, 2005. 5. Kouchakzadeh MA, Sekine H. Compressive buckling analysis of rectangular composite laminates containing multiple delaminations. Composite Structures 2000; 50: 249-255. 6. Kim HJ. Postbuckling analysis of composite laminates with a delaminations. Computers & Structures 1997; 62(6): 975-983. 7. Kim HJ, Hong CS. Buckling and postbuckling behavior of composite laminates with a delaminations. CompositeS Science and Technology 1997; 57: 557-564. 8. Kyoung WM, Kim CG, Hong CS. Buckling and postbuckling behavior of composite cross-ply laminates with multiple delaminations. Composite Structures 1999; 43: 257-274. 9. Hwang SF, Liu GH. Buckling behavior of composite laminates with multiple delaminations under uniaxial compression. Composite Structures 2001; 53: 235-243. 10. Hwang SF, Mao CP. Failure of delaminated interply hybrid composite plates under compression. CompositeS Science & Technology 2001; 61:1513–1527. 11. Cappello F, Tumino D. Numerical analysis of composite plates with multiple delaminations subjected to uniaxial buckling load. Composites Science & Technology 2006; 66: 264-272. 12. Zienkiewicz OC. The finite element method in engineering science. McGraw-Hill, London, 1991. 13. Kyoung WM, Kim CG. Delamination buckling and growth of composite laminated plates with transverse shear deformation. Journal of Composite Materials 1995; 29(15): 2047–2068.