Delamination Of Composites

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Finite Element Modelling of Delamination Buckling of Composite Panel Using ANSYS S.Rajendran and D.Q.Song Materials Technology Application Centre Singapore Productivity and Standards Board 1, Science Park Drive Singapore 118221 Abstract Delamination buckling analysis of composite panels is of considerable interest to aerospace industries. In this paper, finite element modelling of delamination buckling of composite panels is discussed. ANSYS 5.4 has been used for modelling the delamination buckling. A 3-D model with 8-node composite shell element is used. The panel is hypothetically divided into two sub-laminates by a plane containing the delamination. The two sub-laminates are modelled separately using 8-node composite shell element. Appropriate constraint conditions are added for the nodes in the undelaminated region using Coupled Nodes facility of ANSYS. The nodes in the delaminated region, whether in the top or bottom laminate, are left free. Using this modelling approach, a few typical test problems have been solved. The computed buckling loads and strain energy release rate values for the test problems tally closely with that of thoery and other researchers. In addition to the test problems, some results on delamination buckling of a woven-fabric carbon-epoxy composite panel are also presented. The two sub-laminate model discussed here provides a convenient approach to delamination buckling analysis.

1. Introduction Delamination in composite structures can be a serious threat to the safety of the structure. Delamination leads to loss of stiffness and strength of laminates under some conditions. This is particularly so in the case of compressively loaded structures as the loss of stiffness may lead to buckling, the consequences of which can be catastropic. Causes of delamination are many. In aerospace applications, this includes manufacturing defects, as well as operationally induced defects such as bird strikes, hits due to runway debris and tool drops.. The type of delamination that is dealt with in this report is the one that is already initiated by one of the above causes. When a laminate is subjected to in-plane compression, the effects of delamination on the stiffness and strength may be characterised by three sets of analytical results: a) Buckling load b) Postbuckling solutions under increased load c) Results concerning the onset of delamination growth and its subsequent development. Many of the analytical treatments deal with a thin near surface delamination. Such approaches are known as “thin-film” analysis in the literature. The thin-film analytical approach may involve significant errors in the post-buckling solutions. Naganarayana and Atluri (1995) have analysed the buckling behaviour of laminated composite plates with elliptical delaminations at the centre of the plates using finite element method. They propose a multi-plate model using 3-noded quasi-conforming shell element, and use J-integral technique for computing pointwise energy release rate along the delamination crack front.

2. Description of the problem considered for study A square panel with central circular delamination is considered. The radius of delamination (r) is 40 mm (Fig.1). Various panel sizes in the range 100 x 100 mm to 800 x 800 mm are considered. Exploiting the symmetry of the panel, a quarter of the panel is modelled for finite element analysis. Two locations of delamination through the thickness are considered (Fig.2) : (a) Between first and

1

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

second ply (b) At the middle plane of the laminate. Throughout this paper, the material properties and the layup details for the analysis shown below are used: Material: VICOTEX Woven fabric carbon-epoxy (Ciba Geigy) E1=48.285 GPa , E2=44.818 GPa, G12=2.563 GPa, γ12=0.0464 Number of layers: 6 Layup: [06] Laminate thickness: 1.85 mm The loading and boundary conditions correspond to uni-axial loading as shown in Fig. 3. The buckling load of the first mode and strain energy release rate under postbuckled state are the information we wish to compute using FEM.

Fig.1. Square composite panel with circular delamination

3. Finite element model The finite element delamination analysis is carried out using ANSYS finite element software (ver.5.4). As discussed in section 2, there are several ways in which the panel can be modelled for the delamination analysis. For the present study, a 3-D model with 8-node SHELL 99 composite element of ANSYS is used. The panel is divided into two sub-laminates (Fig.4) by a hypothetical plane containing the delamination. For this reason, the present finite element model would be referred to as two sub-laminate model. The two sub-laminates are modelled separately using 8-node SHELL 99 composite element, and then joined face to face with appropriate interfacial constraint conditions for the corresponding nodes on the sub-laminates, depending on whether the nodes lie in the delaminated or undelaminated region. Typically, a node in the undelaminated region of bottom sub-laminate and a corresponding node on the top sub-laminate are declared to be coupled nodes using Coupling/Constraint Equation facility of ANSYS. The nodes in the delaminated region, whether in the top or bottom laminate, are connected by Contact 52 element. This would mean that the two sublaminates are free to move away from each other in the delaminated region, and constrained to move as a single laminate in the undelaminated region.

2

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Fig.2. Delamination locations through the thickness

Fig.3. Boundary conditions and loads

3

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

While building the two sub-laminate model in ANSYS, manual declaration of coupled nodes, pair by pair, is time consuming and is troublesome, as this involves selecting the nodes using the cursor. The option of Coincident Nodes under using Coupling/Constraint Equation facility is useful for automatic generation of coupled nodes data. In order to enable this, the top and bottom sub-laminates of undelaminated region are first off-set by +δ and −δ, respectively, from the delamination plane. The value of δ is chosen to be very small real number. Then, the top and bottom sub-laminates of delaminated region are off-set by +(δ +ε) and −(δ +ε), respectively, where ε is another small real number. The nodes are then declared coincident using a tolerance value in the range 2δ
Fig.4. Two sub-laminate delamination model When we use the shell element, it is important to decide where to position the nodes across the thickness of the element. In the two sub-laminate model described above, it is essential that the nodes of the individual laminates lie in the plane passing through the delamination. This is necessary for a smooth transition of bending behaviour of sub-laminates from the undelaminated region to delaminated region. In ANSYS, this can be achieved using the node shifting option. Using this option, for all the elements at the top laminate, the nodes are shifted to bottom face of the element, and vice versa.

The finite element model (Fig.5) typically has 2468 SHELL 99 elements and 7618 nodes. The region close to delamination front has a finer mesh to facilitate computation of strain energy release rate later.

4

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Fig.5. Typical finite element mesh. No. of elements=2468 SHELL 99, No. of nodes=7618 The following are the important assumptions in the present finite element analysis: a) The material properties are homogeneous. b) There is no material nonlinearity. c) The geometry is perfectly square and symmetric. d) The loading is symmetric, i.e. the loads always act at the mid-plane location, unless otherwise specified.

4. Validation of finite element model 4.1 Test problem no. 1: Buckling load computation The aim of this test problem is as follows: a) To check if the mesh refinement is sufficient for accurate buckling load computation. b) To check if the coupled nodes facility of ANSYS is useful for the present problem. Geometry: 0.09 m x 0.09 m square plate of thickness 0.00185 m. Material: Isotropic with Young’s modulus=2.1x1011 N/m2, Poisson’s ratio=0.3 Finite element mesh: Using symmetry, only one fourth of the plate is modelled. The plate is considered to be made of two layers of equal thickness and of isotropic material properties listed above. A two sub-laminate model is considered. The corresponding nodes at the top and bottom laminates are declared as coupled nodes using the procedure already described in Section 3. The finite element mesh is created using SHELL 99 element of ANSYS. Edge D

5

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Boundary condition: (1/4 th symmetry model)

y

Edge A

Edge B

x

Edge C Symmetry boundary condition are applied on edges A and C; uz is suppressed on edge B. The edge D is left free. θx is suppressed for all nodes so as to simulate a beam-like behaviour of the plate. This is done to facilitate comparison with theretical prediction which employs this assumption. Load: Compressive load at edge B. ANSYS Analysis Procedure: Linear buckling analysis is carried out. This is a two run analysis in ANSYS: 1) A static analysis run with Prestress effects on: For this analysis a unit in-plane compressive load intensity (1 N/m) is applied on edge B. This run essentially creates the matrices required for the eigenbuckling analysis run. 2) Eigenbuckling analysis run: This run computes the buckling load multiplier, λ. If unit in-plane compressive load intensity is applied in the first run, then λ simply represents the buckling load. If the in-plane compressive load intensity is P (N/m), then the buckling load is given by λP. Results of Analysis: The buckling load (without delamination) computed using ANSYS is 13330.0 N. The theoretical buckling load can be computed using the expression (Bruhn, 1973): σcr=(π2E/12(1-γ2)) (t/b)2 (1) Where E is the Young’s modulus, γ is the Poisson’s ratio, t is the thickness and b is the side of the square plate. For this problem, E=2.1x1011 N/m2, γ=0.3, t=0.00185 m and b=0.09 m. This yields a buckling load of 13352.7 N. The difference between theory and FEM is only 0.2% which is quite good. 4.2 Test problem no. 2: Strain energy release rate computation ANSYS has built-in facility for strain energy release rate computation for some 2D analysis. However, such a facility is not available for 3D analysis. Hence, a FORTRAN program has been developed to postprocess ANSYS results in order to compute strain energy release rate. Modified virtual crack closure technique (Rybicki and Kanninen, 1977; Davidson and Schapery, 1990) is used for computing strain energy release rate. For plate/shell elements, Robinson et al. (1994, 1995) have successfully employed a particular implementation of modified crack closure technique in which the element resultants are multiplied with the corresponding element displacements. The accuracy of such an implementation has already been verified earlier (Rajendran et al., 1998) for ANSYS SHELL 99 element and hence is used for the present work. The details of the test problem considered are as follows: Geometry: 0.2 m x 0.2 m square plate of thickness 0.0185 m. Radius of delamination is 0.04 m. Material: Isotropic with Young’s modulus=50 GPa, Poisson’s ratio=0.3 Finite element mesh: Using symmetry, only one fourth of the plate is modelled. Here again, the plate is considered to be made of two layers. The thickness of the thinner layer is taken as 0.00185 m (and hence that of the thicker layer is 0.0185-0.00185 m). Both the layers are assumed to be of isotropic material properties listed above. A two sub-laminate model is considered. The corresponding nodes at the top and bottom laminates in the undelaminated region of the plate are declared as coupled nodes using the procedure already described in Section 3. The finite element mesh is created using SHELL 99 element of ANSYS. Boundary condition: (1/4 th symmetry model)

6

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Edge D y

Edge B

Edge A

x

α Edge C Symmetry boundary condition are applied on edges A and C; uz and θx are suppressed on edge B. uz and θy are suppressed on edge D. Load: Compressive load at edges B and D. ANSYS Analysis Procedure: Linear buckling analysis is carried out in the same way as in the previous test problem. Results of Analysis: This problem is a typical thin film delamination problem. The thickness of the film is 0.00185 m which is 0.1 times the thickness of the laminate, 0.0185. Hence, the thicker laminate may be considered to be infinitely stiff compared to thinner laminate. Under this condition, the delaminated portion of the thinner laminate may be considered as a clamped circular plate with uniform radial compressive stress, in which case, the buckling stress is given (Evans and Hutchinson, 1984) by

σc=1.2233[E/(1-γ2)](t2/a)2

(2)

where E is the Young’s modulus, γ is the Poisson’s ratio, t2 is the thickness of the thinner laminate, and a is the delamination radius. This, therefore, also represents the buckling stress of local delamination. The buckling stress computed by this formula is 3993765 N/m2 whereas the finite element analysis gives a value of 3989838 N/m2. The difference is less than 0.1% which is quite good. When the postbuckling deformation is axi-symmetric and nearly linear in the neighbourhood of buckling point, the pointwise energy release rate is given (Evans and Hutchinson, 1984) by G=(σ2-σc2) (1-γ2) t2/ ((1.8285+γ)E)

(3)

where σ is the applied stress at which the energy release rate is required to be computed. Fig. 6 shows the energy release rate distribution plotted against the angle, α, computed using finite element analysis for σ = 1.2 σc. The plot shows a constant value of G as expected from Eq.(3). However, the average value of G computed by FEM is 0.02069 N/m whereas the value given by Eq.(3) is 0.01847 N/m. The FEM value is about 12% higher. Such a deviation has also been observed by Naganarayana et al. (1996). The reason for deviation is attributed to the assumption in Eq.(3) that the thicker laminate is infinitely rigid whereas in finite element modelling it is still flexible, and hence the deviation of 12% is not of concern. The important point is that the strain energy distribution obtained by finite element analysis is constant as predicted by Eq.(3). Thus, the model is adequate with respect to prediction of delamination buckling load and strain energy release rate.

7

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

2.50E-02 2.00E-02 1.50E-02 1.00E-02 5.00E-03 0.00E+00 0

10

20

30

40 50 60 Angle (Deg.)

70

80

90

Fig. 6. Distribution of strain energy release rate at the delamination front for the test problem

5. Delamination buckling analysis of composite panel In the present study, we are concerned with only the first buckling load and the corresponding buckling mode. Using ANSYS, buckling analysis can be carried out by two approaches: (a) Eigenvalue (or linear) buckling analysis (b) Nonlinear (buckling) analysis. The eigenvalue buckling analysis is simple and faster than the nonlinear analysis. It yields the buckling load straightway. The nonlinear buckling analysis is essentially a nonlinear static analysis considering large deflections. However, this analysis does not straightway yield the buckling load. The buckling load is inferred rather than calculated. The load deflection data from the nonlinear analysis is plotted and the load at which the deflection rapidly rises is taken as buckling load. 5.1 Linearised (eigenvalue) buckling analysis For all the panels with mid-plane delamination analysed, the buckling is predominantly global, and the buckling modeshape is similar to Fig.7. However, in the case of panels with close-to-surface delamination, the buckling modeshape is predominantly local as seen in Fig.8 for small r/a values. However, for low values of r/a, the modeshape involves both global and local buckling as seen in Fig.9.

Fig.7. Typical global buckling mode

8

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

In Fig. 10, the horizontal axis represents the ratio of delamination radius to half the panel side, r/a, while the vertical axis represents the ratio of buckling load of panel with delamination to the buckling load of panel without delamination, Pcr/Pocr. Figure 10 shows that as r/a increases, the nondimensional buckling load, Pcr/Pocr, decreases for both the cases, t2/t1=1.0 and 0.2. The decrease is more pronounced for t2/t1=0.2 than for t2/t1=1.0.

Fig.8. Typical local buckling mode

Fig.9. Typical buckling mode which involves both global and local buckling

9

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

1 0.9 0.8

P cr /P o cr

0.7 0.6 0.5

Delamination at t2/t1=1.0

0.4

Delamination at t2/t1=0.2

0.3 0.2 0.1 0 0

0.2

0.4

0.6

0.8

1

r/a Fig.10. Non-dimensional buckling for various delamination size and location

5.2 Nonlinear analysis 5.2.1 Comparison of buckling loads from linear and non-linear buckling analysis Figures 11 and 12 summarise the results of nonlinear analysis. The solution path has been traced for three values of eccentricity, viz. one-, two- and three-ply-thickness. For values of eccentricity lower than one ply eccentricity, the solution had convergence problems particularly near the bent portion of the curve (Fig.11), and hence these solutions are not available. The trend of the curve suggests that the nonlinear buckling load is around 4000 N/m which is fairly close to the linear buckling load of 4042.2 N/m. Absolute comparison of nonlinear and linear buckling load in terms of percentage difference is less meaningful as the nonlinear buckling load, 4000 N/m, is only an estimate. It may be concluded that the agreement is quite good. 5.2.2 Post buckling solution and strain energy release rate at the delamination front For the computations reported in this section, the same problem as summarised in Table 1 is considered. The postbuckling solutions and the strain energy release rates have been computed for various eccentricity and load factor values. The load factor is defined as the ratio of actual applied load, P, to the linearised buckling load, Pcr. The results are summarised in Fig.12. The strain energy release rate (G) shown in this figure is the total strain energy release rate considering all the three modes. As expected, the strain energy release rate values increase with load factor value. For low load factor values, i.e. the applied load close to the buckling load, the strain energy release rate is close to zero, particularly for low eccentricity of loading. It is also seen from Fig.12 that, as the eccentricity of loading increases, the strain energy release rate values also increase.

10

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Table 1 Summary of Input Data for nonlinear buckling analysis Lay-up: [O6] Laminate thickness = 0.00185 m; Ply thickness = 0.00185/6 m = 0.00030833 m Ratio of thickness of thinner sub-laminate to that of thicker sub-laminate, t2/t1 = 1/2 Material properties: Ex=48.285 GPa , Ey=44.818 GPa, Gxy=2.563 GPa, γxy=0.0464 Panel size: 0.2 m×0.2 m Radius of delamination=0.04 m

y x

Finite element model: One fourth of the panel is modelled as shown below: D A

B

Displacement at panel centre (m)

C Boundary Condition: a) Symmetry boundary condition on edges A and C. b) uz = θx = 0 on Edge B; Load: Compressive pressure is applied on edge B. The maximum intensity of pressure applied (P) is 1.1 times the buckling pressure calculated by linear buckling analysis (Pcr). For this problem, Pcr= 4042.2 N/m.

3.50E-02 Thicker sublaminate, half ply eccentricity Thinner sublaminate, half ply eccentricity Thicker sublaminate, one ply eccentricity Thinner sublaminate, one ply eccentricity Thicker sublaminate, three ply eccentricity Thinner sublaminate, three ply eccentricity

3.00E-02 2.50E-02 2.00E-02 1.50E-02 1.00E-02 5.00E-03 0.00E+00 0

1000

2000

3000

4000

Applied compressive load intensity (N/m) Fig.11. Nonlinear solution and the effect of eccentricity of loading

11

5000

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

500 Half ply eccentricity, load factor=1.100

450

Half ply eccentricity, load factor=1.056

400

Half ply eccentricity, load factor=1.012 Half ply eccentricity, load factor=0.968

350

One ply eccentricity, load factor=1.100

G (N/m)

300

One ply eccentricity, load factor=1.056 250

One ply eccentricity, load factor=1.012 One ply eccentricity, load factor=0.968

200

150

Three ply eccentricity, load factor=1.012

100

Three ply eccentricity, load factor=1.056 Three ply eccentricity, load factor=1.012

50

Three ply eccentricity, load factor=0.968

0 0

20

40

60

80

100

Angle (Deg.)

Fig.12. Strain energy release rate distribution along the crack front 7. Concluding remarks Buckling analysis of square composite panels with embedded circular delaminations has been carried out using ANSYS 5.4. The following are some important observations: 1.

The two-sublaminate model developed in this work provides a convenient method to model delaminated composite panels. The accuracy of the model in predicting the buckling loads and strain energy release rates has been verified using a typical test problem. The agreement of results with theory is quite good.

2.

Whenever the delamination is located at the mid-plane of laminate, the panel exhibits only global buckling, i.e. there is no local buckling of delaminated region. On the other hand, whenever the delamination is close to the surface (typically t2/t1=0.2), the buckling mode is predominantly local, i.e. only the delaminated region buckles for small r/a values (typically, 0.3). For higher r/a values, both global and local buckling co-exists.

3.

The reduction of buckling load is more significant for close-to-surface delamination than for midplane delamination.

4.

The eccentricity of loading significantly influences the buckling load as well as the strain energy release rate distribution at the delamination front. Hence, delamination buckling studies must consider eccentricity of loading particularly in the case of thin panels such as the ones considered in this report.

12

Proceedings of 2nd Asian ANSYS User Conference, Nov 11-13, 1998, Singapore

Acknowledgement - The funding by National Science and Technology Board (NSTB) of Singapore for this project is greatly acknowledged. The authors also would like to thank other members of the team, A.Lee and C.K.Tan of SPSB for their invaluable help in experiments, and S.P.Phua and H.F.Lee of ST Aero for their comments and discussions during the project.

References ANSYS User’s Manual Vol. I (1994), Swanson Analysis Systems, Inc. Bruhn, E.F. (1973), Analysis and Design of Flight Vehicle Structures, Tri-state offset co. (USA). Davidson, B.D. and Schapery, R.A. (1990), A technique for predicting mode-I energy release rates using a first order shear deformable plate theory, Engng. Fracture Mech., 36, 157-165. Evans, A.G. and Hutchinson, J.W. (1984), On the mechanics of delamination and spalling in compressed film, International Journal of Solids and Structures, 20, 455-466. Naganarayana, B.P. and Atluri, S.N. (1995), Strength reduction and delamination growth in thin and thick composite plates under compressive loading, Computational Mechanics, 16, 170-189. Rajendran, S. and Song, D.Q. (1998), Finite element delamination buckling analysis of composite panels, Part 1 - Some modelling considerations, a part of Quarterly Report of ATP R&D Project 9. Robinson, P. and Song, D.Q. (1994), The development of an improved mode III delamination test for composites, Composite Science and Technology, 52, 217-233. Robinson, P., Javidrad, F. and Hutchings, D. (1995), Finite element modelling of delamination growth in the DCB and edge delaminated DCB specimens. Composite Structures, 32, 275-285. Rybicki, E.F. and Kannenen, M.F. (1977), A finite element calculation of stress intensity factor by modified crack closure integral, Engrg. Fracture Mech., 9, 931-938.

13