Finite element modeling of composite plates with internal delamination

Composite Structures 90 (2009) 21–27

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Finite element modeling of composite plates with internal delamination K. Alnefaie Mechanical Engineering Department, King Abdulaziz University, P.O. Box 80248, Jeddah 21589, Saudi Arabia

a r t i c l e

i n f o

Article history: Available online 27 January 2009 Keywords: Laminated composites Plates Finite element model Internal delamination Natural frequencies

a b s t r a c t A three-dimensional (3D) finite element model of delaminated fiber-reinforced composite plates is developed to analyze their dynamics. Natural frequencies and modal displacements are calculated for various case studies with different dimensions and delamination characteristics. Numerical results showed a good agreement with available experimental data. A new proposed model shows enhancement of the accuracy of the results. This study also introduces a method for detecting delamination in composite plates. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The use of laminated composite materials in space vehicles and various machine components has increased considerably over the past decades. In fact, a monograph written four decades ago showed very limited number of publications on laminated plate vibration. There is evidence that vibration of laminated plates was the subject of more than a 1000 publication in the last decade [1,2]. Researchers show significant interest in determining the natural frequencies of such important structures. Qatu published results for the frequencies of laminated plates having rectangular shapes under different boundary conditions [3,4] and triangular shapes [5]. Under repeated or impact loads these materials are subjected to various forms of damage, mostly delaminations and cracks [6–11]. Such damage becomes an obstacle to the more extensive usage of composite materials. Therefore, the monitoring of internal or hidden damage in composite material is critical in engineering practice [12]. The use of vibration based techniques as nondestructive testing methods for damage monitoring of laminated composite panels is a field attracting the interest of many researchers [13–19]. The effective damage monitoring for this kind of material or structure depends largely on the accurate prediction or estimation of mechanical or dynamic behaviors of both intact and damaged composite panels [20]. It is difficult to obtain accurate exact solutions for multi-layered panels having arbitrary lamination sequence and/or boundary conditions. This difficulty increases considerably when such structures are delaminated. Thus, computational approaches like finite element methods play an important role in detecting damage for laminated composites. There are

E-mail address: [email protected] 0263-8223/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2009.01.004

numerical and experimental investigations of delaminated multilayered composites [19–26]. Sankar [27] modeled a delaminated beam as two sublaminates by offsetting beam finite elements. Rikards [28] developed a model of finite superelements for sandwich composite beam and plate without delamination, each layer being considered as a simple Timoshenko beam. Later, Gadelrab [8] discussed the modal variation of delaminated beam under different boundary conditions. Ousset and Roudolff [9] analyzed the delaminated multi-layer composite plate based on Mindlin Reissner plate model. Zak et al. [29] and Ousset and Roudolff [9] developed models of finite elements for beams and plates with boundary delamination. Among most of the publications the prediction of material mechanical or dynamic behaviors is based on the classical laminated plate theory [3,30]. In this theory, the transverse shear deformation effect is ignored. Subsequently, this theory cannot provide accurate results for moderately thick laminated plates, for which the inplane elastic modulus is much higher than the transverse shear modulus [3,31]. Besides, the Poisson’s effect is significant for angle-ply laminated plates. Therefore, three-dimensional layer-wise theory is needed in order to obtain an accurate prediction of the dynamic response of multi-layered composite panels. In addition, potential dangers are often induced by hidden or internal delaminations in the in-service laminated composites. However, there is a lack of both numerical and experimental investigations on internal delaminations with different geometries in laminated composites. Yam et al. [34] discussed the relation between the strain mode and the displacement mode for vibrating elastic structure. They claimed modal strain analysis can be used for investigation of local variation of structural parameters. In this paper, a three-dimensional finite element model for multi-layered composites with internal delamination is established, and the fiber orientations of individual lamina as well as the transverse shear effect are taken into account. Numerical calculations

22

K. Alnefaie / Composite Structures 90 (2009) 21–27

using ABAQUS/CAE v6.3 package are carried out for different plates under free boundary conditions. Natural frequencies, modal displacements of the intact and damaged multi-layer composite plates are subsequently analyzed for various samples.

8 X ½Ni fdi g

Z

2

ee11 ; ee22 ; ee33 ; ee12 ; ee13 ; ee23

ð7Þ

t13 E1 t23 E2

0

0

0

0

1 E3

0

0

0

0

1 G12

0

0

0

0

1 G13

0

0

0

0

1 E2 t23 E2

½C ¼ 6 6 0 6 6 6 0 4 0

31

0

7 0 7 7 7 0 7 7 7 0 7 7 7 0 7 5

ð8Þ

1 G23

is the matrix of material constants with E1, E2, E3, G12, G13, G23, t12, t13 and t23 being the orthotropic elastic constants of individual lamina, and

6 6 6 6 6 ½A ¼ 6 6 6 6 4

fee g ¼

t12 E1

1 E 6 t1 6 12 6 E1 6 t13 6 6 E1

where {d} = (ui, vi, wi)T is the displacement vector at node i, [Ni] = Ni[I3], [I3] is the three-order unit matrix and Ni the shape function [32]. Then the strain vector of each element can be expressed in terms of displacement in the global coordinate system as 8 X ¼ ½Bi fdi g

½BT ½A½C½A1 ½BdV

where

2

T

cos2 h

2

sin h

0

0

0

2 sin h cos h

sin h 0

cos2 h 0

0 1

0 0

0 0

2 sin h cos h 0

0

0

0

cos h

sin h

0

0

0

0  sin h cos h

2

sin h cos h  sin h cos h 0

ð2Þ

0

where

ð3Þ

2

@ @x

6 0 6 6 6 0 6 ½ D ¼ 6 @ 6 2@y 6 6 0 4 @ 2@z

0 @ @y

0 @ 2@x @ 2@z

0

0

3

0 7 7 7 @ 7 @z 7 7 0 7 7 @ 7 2@y 5

ð4Þ

ð5Þ T

T T

with [B] = [[B1], [B2], . . ., [B8]] and {d } = ({d1} , . . ., {d8} ) . Therefore, the stresses of an element in global coordinate system can be expressed by the nodal displacement as

Vk

x2

y

p l

Fiber θ

o

U ek ¼

x

Fig. 1. The cross-ply laminated plate and the coordinate system.

N  e T  e  e  1X d K k dk 2 k¼1 k

ð11Þ

1 fdgT ½Kfdg 2

U¼

ð12Þ

where {d} and [K] are the global nodal displacement vector and stiffness matrix, respectively. Similarly, the kinetic energy of the kth element is

T¼

k

j

ð10Þ

1 2

Z Vk

n oT

n oT

n o

1 q d_ ek dV ¼ d_ ek ½Mek  d_ ek 2

ð13Þ

where fd_ ek g and ½Mek  represents the velocity vector and the mass matrix of the kth element, respectively.The total kinetic energy for a composite plate consisting of N elements can be expressed as N X k¼1

n

m

1  e T e e d ½K k fdk g 2 k

Therefore, after assembly of nodal displacements of all elements, the total strain energy of a multi-layer composite plate can be represented as

T ek ¼

Z

i

fee gT fre gdV ¼

k¼1

e

x3

Z

1 2

N X

@ 2@x

fee g ¼ ½Bfde g

2

where fdek g and ½K ek  represents the displacement vector and the stiffness matrix of the kth element, respectively. The total strain energy for a composite plate consisting of N elements can be expressed as

U¼

Thus, the strain vector in an element can be expressed in terms of nodal displacements as

0 cos2 h  sin h

is the transform matrix between the local and global coordinate systems.Then, the strain energy of the kth element is

U ek ¼

and

0

ð9Þ

i¼1

½Bi  ¼ ½D½N i 

ð6Þ

where [Ke] is the element stiffness matrix as

ð1Þ

i¼1



8 X ½K e fde g

Ve

The finite element model used for studying the dynamic behavior of multi-layered composite plates (Fig. 1) is an eightnoded rectangular solid thin plate element. For each node, there are three degrees of freedom, i.e. translations along the global coordinate axes of x, y, and z. The element thickness is assigned to be equal to that of the corresponding individual lamina. The local element coordinate system (x1, x2, x3) is arranged with the first axis being coincident with the fiber direction. All physical parameters throughout an element are assumed to be the same. For an eight-node finite element with three degrees of freedom per node, the displacement field over an element is given by

fdg ¼ ðu; v ; wÞT ¼



re11 ; re22 ; re33 ; re12 ; re13 ; re23 T ¼

i¼1

½K e  ¼

2. Finite element modeling for laminated composite plates



fre g ¼

x1

T ek ¼

N n oT 1X d_ e ½M ek fd_ ek g 2 k¼1 k

ð14Þ

Therefore, after assembly of nodal velocities of all elements, the total kinetic energy of a multi-layer composite plate can be represented as

T¼

1 _ T _ fdg ½Mfdg 2

ð15Þ

3 7 7 7 7 7 7 7 7 7 5

23

K. Alnefaie / Composite Structures 90 (2009) 21–27

_ and [M] are the global nodal velocity vector and mass where fdg matrix, respectively. If the composite plate experiences a harmonic motion with angular frequency x, the kinetic energy of the composite plate will be

T¼

400 is negligible. This leads to applying the 400 elements per layer for each lamina model in the subsequent computations to achieve high accuracy and minimum computations time Table 2. 2.2. Model verification

1 2 x fdgT ½Mfdg 2

ð16Þ

Therefore, the modal parameters such as natural frequencies xi, mode shapes, modal strains, etc., can be obtained. Pairs of nodes are modeled at the same location where a delamination is to be induced. For an arbitrary intact laminated plate, in order to ensure the material continuity, the displacements and their variations for each pair of coincident nodes on the upper and lower adjacent laminae have to be equal in the whole process of computation. When the plate is delaminated, the pairs of nodes are detached and their displacements on the upper and lower surfaces within the delamination region are not to be connected to each other.

To verify the results of the developed model, the results obtained here are compared with the numerical results of Yam et al. [33] and both the numerical and experimental work of Lin et al. [35]. The present analysis shows closer results to the reported experimental results than the other reported analytical results. In addition, another laminated composite plate is considered for model verification. This is a 16-layer plate with an area of 240  180 mm2 and a total thickness of 2.08 mm. The ply orientations are of [0°/0°/90°/90°/0°/0°/90°/90°]s and the material constants are E1 = 125 GPa, E2 = E3 = 8.5 GPa, G12 = G13 = 4.5 GPa, G23 = 3.27 GPa, t12 = t13 = t23 = 0.3 and q = 1550 kg m3, Wei et al. [36]. The first six natural frequencies of this plate are compared with the numerical and experimental work by Wei et al. [36]. Again the present model demonstrates better agreement with the experimental results reported by Wei et al. [36] than the analytical results reported.

2.1. Model refinement

3. Results and discussions

A sensitivity analysis is carried out to study the effect of element size on the accuracy of the finite element results. In this study, an 8-layer laminated square composite plate with a side length of 178 mm and a total thickness of 1.58 mm is considered. All the ply orientations are equal to 0° and the material constants are E1 = 172.7 GPa, E2 = E3 = 7.2 GPa, G12 = G13 = 3.76 GPa, G23 = 2.71 GPa, t12 = t13 = 0.3, t23 = 0.33 and q = 1566 kg m3, Yam et al. [33] and Lin et al. [35]. Table 1 shows the natural frequencies of the first six modes for different generated mesh sizes. The attained results are representing variable number of elements per layer. A remarkable change in results is obtained when increasing the total number of elements from 25 to 100 elements per layer. Meanwhile, the change in natural frequency obtained with 900 elements in comparison with

3.1. Samples

Then, using the Lagrange’s principle the equation of motion for free vibration of the composite plate is reduced to the eigenvalue problem of

ð½K  x2 ½MÞfdg ¼ 0

ð17Þ

Table 1 Natural frequencies (Hz) of the intact [0°/0°/0°/0°]s free plate with a side length of 178 mm and a total thickness of 1.58 mm using different meshes. Mode

1 2 3 4 5 6

Number of elements for each lamina 25

100

200

400

900

78.51 100.17 277.31 316.42 489.29 581.06

81.23 107.20 207.72 294.01 422.93 523.91

81.36 109.50 200.84 298.13 396.95 530.36

81.48 109.20 199.50 300.46 391.40 533.89

81.56 109.58 199.51 301.73 391.12 535.80

Five square plates are considered made of multi-layer carbon fiber-reinforced epoxy composites. Every sample plate has a side length of 225.5 mm and a thickness of 2.05 mm and consists of 8-layer in orientation of [0°/90°/0°/90°]s. One of the plates is intact and named O, the other plates are named A, B, C, D, respectively, are damaged with delamination of different areas. The delamination location is simulated as an interspace, or actual separation, of 0.02 mm between the third and the fourth layers counting from the top of the plate. Every delaminated plate has only one square delamination located at the position with the center as shown in Fig. 2. The delamination center coordinates are x = 163.5 mm, y = 163.5 mm and z = 1.28125 mm. The delamination areas of plates A, B, C and D are 11.275  11.275, 33.825  33.825, 56.375  56.375 and 78.925  78.925 mm2, respectively. The material constants of the samples for FEM computations are E1 = 37.78 GPa, E2 = E3 = 10.9 GPa, G12 = G13 = G23 = 4.91 GPa, t12 = t13 = 0.3, t23 = 0.11 and q = 1813.9 kg m3, Yam et al. [33]. 3.2. Delamination effect on the natural frequencies The natural frequencies are computed for the first ten modes of the plates. Table 3 lists the natural frequencies for the plates with different sizes of delamination. It can be seen that with the increase of delamination area, the natural frequency decreases. The

Table 2 Natural frequencies (Hz) of the intact [0°/0°/0°/0°]s free plate and the intact [0°/0°/90°/90°/0°/0°/90°/90°]s free plate for different methods. Mode

[0°/0°/0°/0°]s plate Present

1 2 3 4 5 6

81.48 109.20 199.50 300.46 391.40 533.89

[0°/0°/90°/90°/0°/0°/90°/90°]s plate Yam et al. [33]

Lin et al. [35]

Present

Numerical

Numerical

Experimental

82.26 113.10 207.29 325.28 408.51 539.92

83.57 118.42 207.79 329.41 419.83 546.93

81.50 107.40 196.60 285.50 382.50 531.00

89.16 278.97 330.35 354.92 393.26 574.50

Wei et al. [36] Numerical

Experimental

90.52 279.17 333.59 354.22 397.62 583.71

90 289 318 354 386 570

24

K. Alnefaie / Composite Structures 90 (2009) 21–27

Delamination Center y x 225.5 163.5

225.5

z

163.5

Fig. 2. The finite element model of the laminated plate and the delamination region center.

Table 3 Numerical results of the natural frequencies (Hz) for the free plates. Delamination area (mm2) 11.275  11.275

33.825  33.825

56.375  56.375

78.925  78.925

Mode

Plate O

Plate A

Plate B

Plate C

Plate D

1 2 3 4 5 6 7 8 9 10

70.09 133.87 168.73 195.86 218.30 360.34 373.41 426.28 460.30 508.41

70.08 133.84 168.67 195.82 218.23 360.27 373.41 426.22 460.30 508.42

70.08 133.57 168.12 195.51 217.71 359.62 373.05 425.51 459.69 507.67

70.07 132.71 166.71 194.53 216.15 357.12 370.37 422.44 454.77 502.99

70.03 130.88 162.47 192.51 213.32 351.20 358.13 413.77 432.76 488.42

delamination effect on the natural frequencies is however very small for the first six modes and higher for the other modes. The relationship between the frequency and delamination size is investigated. Figs. 3 and 4 show the percentage changes in the natural frequencies as a function of the delamination areas, where the

height of the column represents the absolute value of the percentage change of natural frequency, i.e.

j¼

jxdelaminated  xintact j

Mode 1 Mode 2 Mode 3 Mode 4 Mode 5

5

ð18Þ

4

3

3

2

2

1

1

A

Mode 6 Mode 7 Mode 8 Mode 9 Mode 10

5

κ (%)

4

0

 100

6

6

κ (%)

xintact

B

C

D

Fig. 3. Percentage changes of the natural frequencies for plates with delamination of different areas for modes 1–5.

0

A

B

C

D

Fig. 4. Percentage changes of the natural frequencies for plates with delamination of different areas for modes 6–10.

25

K. Alnefaie / Composite Structures 90 (2009) 21–27

0.8

0.018 Mode 1 Mode 2 Mode 3

0.016 0.014

0.7 0.6

0.012

0.5

0.01

δφ

δφ

Mode 7 Mode 8 Mode 9 Mode 10

0.008

0.4

0.006

0.3

0.004 0.2

0.002 0

1

2

3

4

5

6

0.1

1

2

3

4

5

Fig. 5. Displacements along z-direction for points on the upper and lower surfaces within delaminated region of Plate C for modes 1, 2 and 3.

It is noticed that these absolute values increase with the delamination area. It is also seen that the decrease of the natural frequencies is not the same for different modes. The delamination-induced decreases of natural frequencies are relatively large for mode 3 and mode 9, and the change of the natural frequencies is almost negligible in mode 1. However, the variation manner of the values is not the same for each plate. When the delamination area is 11.275  11.275 mm2, the delamination-induced changes of the natural frequencies are nearly zero for all the considered cases, which indicates that the delamination-induced frequency change is insignificant for small delamination. Therefore, it is indispensable to analyze the delamination-induced changes of other parameters such mode shapes, modal strains, etc., for effective detection of delamination in composite plates.

Fig. 7. Displacements along z-direction for points on the upper and lower surfaces within delaminated region of Plate C for modes 7, 8, 9 and 10.

Fig. 8. The first mode shape of Plate C (70.07 Hz).

3.3. Delamination effect on mode shapes The above analysis demonstrates that the delamination-induced changes of plate parameters are mode-dependent. This may imply that the delamination region exerts specific effects on the relevant modes. In order to further investigate the relationship

0.045 Mode 4 Mode 5 Mode 6

0.04 0.035

Fig. 9. The second mode shape of Plate C (132.71 Hz).

δφ

0.03 0.025 0.02 0.015 0.01 0.005 0

1

2

3

4

5

6

Nodes within delamination region

Nodes within delamination region

6

Nodes within delamination region Fig. 6. Displacements along z-direction for points on the upper and lower surfaces within delaminated region of Plate C for modes 4, 5 and 6.

Fig. 10. The third mode shape of Plate C (166.71 Hz).

26

K. Alnefaie / Composite Structures 90 (2009) 21–27

Fig. 11. The fourth mode shape of Plate C (194.53 Hz).

Fig. 14. The seventh mode shape of Plate C (370.37 Hz).

occurs somewhere in a composite plate, there may be an interactive motion or impact within the delamination region during the vibration of the plate. These phenomena cause the variations of energy dissipation in the plate, and they are modal dependent. Thus, the delamination can be detected according to the variation of energy dissipation in the plate during vibration. The effect of contact on the delaminated stiffened plate was discussed by Chen et al. [37]. 4. Conclusions The delamination problem for multi-layer composite plates has been analyzed using finite elements method and modal analysis. To achieve accurate results for delamination detection of multilayer composite plate, a relatively fine mesh is considered for the finite element computation. The following conclusions may be drawn from the results of numerical simulation in this study.

Fig. 12. The fifth mode shape of Plate C (216.15 Hz).

Fig. 13. The sixth mode shape of Plate C (357.12 Hz).

between the delamination and the mode-dependent variations of energy dissipation in the plate, the unit-normalized local displacements of points within the delamination region of the plates are computed. The relative displacements are analyzed for points along the line (y = 163.5 mm, z = 1.28125 mm) just on the upper and lower surfaces within the delamination region of Plate C (the plate with delamination area of 56.375  56.375 mm2). Figs. 5–7 show the differences of displacement along z-direction between the upper and lower points

du ¼ jwupper  wlower j

ð19Þ

which are assumed to be coincident with each other before the plate vibration starts. It can be seen that significantly higher displacements appear in modes 7–10, when compared with the modes with lower frequencies. Therefore, energy dissipation variation will be larger in the modes for the delaminated plates during vibration. Figs. 8–14 show the mode shapes of Plate C. The delamination appears clearly in mode 7 because of the high energy dissipation in this mode as it is seen in Fig. 7. Therefore, when delamination

(1) The finite element model proposed in this paper can predict accurately the dynamic behaviors of a multi-layer composite plate with internal delamination at arbitrary location. (2) Local internal delamination has slight effect on the natural frequencies of a multi-layer composite plate although the extent of the natural frequency variation increases with both the delamination dimension and the order of the natural frequency. (3) Delamination effect is more clear on mode shapes than it is on the frequencies. The changes of the displacement are mode-dependent. The location of the delamination in the plate is critical. If the delamination region occurs in the plate where the displacement of a certain mode is high, this mode is highly impacted by the delamination. (4) The results of numerical analysis in this study can be taken as guidance for arrangement of displacement measurements on the surface of specimen when experimental modal analysis is carried out to investigate local changes of structural parameters for damage detection. The finite element model, strategy and numerical method provided in this paper can be used for dynamic response analysis of damaged engineering structures, especially for multi-layer composite plates.

References [1] Leissa AW, Vibration of plates, NASA SP-160, US Government Printing Office., Washington, DC, 353 pp., republished 1993, The Acoustical Society of America. [2] Qatu MS. Vibration of laminated shells and plates. Elsevier; 2004. [3] Qatu MS. Free vibration of laminated composite rectangular plates. Int J Solids Struct 1991:941–54. [4] Qatu MS, Leissa AW. Vibration studies for laminated composite cantilever plates. Int J Mech Sci 1991:927–40. [5] Qatu MS. Natural frequencies for cantilevered laminated composite right triangular and trapezoidal plates. Compos Sci Technol 1994:441–9.

K. Alnefaie / Composite Structures 90 (2009) 21–27 [6] Tsai SW, Hann HT. Introduction to composite materials. Westport, CT: Technic Publishing Company; 1980. [7] Voyiadjis GZ. Damage in composite materials. Amsterdam, New York: Elsevier; 1993. [8] Gadelrab RM. The effect of delamination on the natural frequencies of a laminated composite beam. J Sound Vib 1996;197(3):283–92. [9] Ousset Y, Roudolff F. Numerical analysis of delamination in multi-layered composite plates. Comput Mech 2000;20(1/2):122–6. [10] Campanelli R, Engblom J. The effect of laminations in graphite/PEEK composite plates on modal dynamic characteristics. Compos Struct 1995;31:195–202. [11] Tenek L, Henneke E, Gunzburger M. Vibration of delaminated composite plates and some applications to non-destructive testing. Compos Struct 1993;23:252–3. [12] Gerardi TG. Health monitory aircraft. J Intell Mater Syst Struct 1990;1:375–85. [13] Salawu OS. Detection of structural damage through changes in frequency: a review. Eng Struct 1997;19:718–23. [14] Rytter A, Kirkegaard PH. Vibrational-based inspection of a steel mast. In: Proceedings of the 12th international modal analysis conference, Hawaii; 1994. p. 1602–8. [15] Gomes AJMA, Silva JMME. On the use of modal analysis for crack identification. In: Proceedings of the eighth international modal analysis conference, FL USA; 1991. p. 1108–15. [16] Cawley P, Adams RD. The location of defects in structure from measurements of natural frequencies. J Strain Anal 1979;4:49–57. [17] Sanders D, Kim YI, Stubbs RN. Non-destructive evaluation of damage in composite structures using modal parameters. Exp Mech 1992;32:240–51. [18] Ceravolo R, Stefano AD. Damage location in structure through a connectivistic use of FEM modal analyses. Int J Anal Exp Modal Anal 1995;10(3):178–86. [19] Crema LB, Mastroddi F. Frequency-domain based approaches for damage detection and localisation in aeronautical structure. In: Proceedings of the International Modal Analysis Conference; 1995. p. 1322–30. [20] Tappet PM, Snyder TD, Robertshaw HH. Attacking the damage identification problem. In: Proceedings of SPIE – the international society for optical engineering smart structures and materials: smart sensing, processing and instrumentation 2443, 1995. p. 286–94. [21] Raghuram PV, Murty AVK. A high precision coupled bending-extension triangular finite element for laminated plates. Comput Struct 1999;72(6):763–77.

27

[22] Lee J. Free vibration analysis of delaminated composite beams. Comput Struct 2000;74(2):121–9. [23] Rikards R. Interlaminar fracture behaviour of laminated composites. Comput Struct 2000;76(1–3):11–8. [24] Borovkov A, Palmov V, Banichuk N, Saurin V, Barthold F, Stein E. Macro-failure criterion for the theory of laminated composite structures with free edge delaminations. Comput Struct 2000;76(1–3):195–204. [25] Lingen FJ, Schipperen JHA. An efficient parallel procedure for the simulation of free edge delamination in composite materials. Comput Struct 2000;76(5):637–49. [26] Riccio A, Perugini P, Scaramuzzino F. Modelling compression behaviour of delaminated composite panels. Comput Struct 2000;78(1–3):73–81. [27] Sankar BV. A finite element for modeling delaminations in composite beams. Comput Struct 1991;38(2):239–46. [28] Rikards R. Finite element analysis of vibration and damping of laminated composites. Compos Struct 1993;24:193–204. [29] Zak A, Krawczuk M, Ostachowicz W. Numerical and experimental investigation of free vibration of multi-layer delaminated composite beams and plates. Comput Mech 2000;26:215–309. [30] Jones RM. Mechanics of Composite Materials. New York: McGraw-Hill Book Company; 1975. [31] Matsunaga H. Assessment of a global higher-order deformation theory for laminated composite and sandwich plates. Compos Struct 2002;56:279–91. [32] Bathe KJ. Finite element procedures in engineering analysis. Englewood Cliffs, NJ: Prentice-Hall; 1982. [33] Yam LH, Wei Z, Cheng L, Wong WO. Numerical analysis of multi-layer composite plates with internal delamination. Comput Struct 2004;82:627–37. [34] Yam LH, Leung TP, Li DB. Theoretical and experimental study of modal strain analysis. J Sound Vib 1996;191(2):251–60. [35] Lin DX, Ni RG, Adams RD. Prediction and measurement of the vibrational damping parameters of carbon glass fiber-reinforced plastics plates. J Comp Mater 1984;18:132–52. [36] Wei Z, Yam LH, Cheng L. Detection of internal delamination in multi-layer composites using wavelet packets combined with modal parameter analysis. Compos Struct 2004;64:377–87. [37] Chen H, Wang M, Bai R. The effect of nonlinear contact upon natural frequency of delaminated stiffened composite plate. Compos Struct 2006;76:28–33.