Wave Propagation Analysis in Composite Laminates Containing a Delamination Using a Three-Dimensional Spectral Element Method

A three-dimensional spectral element method SEM was developed for analysis of Lamb wave propagation in composite laminates containing a delamination. SEM is more efficient in simulating wave propagation in structures than conventional finite element method FEM because of its unique diagonal form of the mass matrix. Three types of composite laminates, namely, unidirectional-ply laminates, cross-ply laminates, and angle-ply laminates are modeled using three-dimensional spectral finite elements. Wave propagation characteristics in intact composite laminates are investigated, and the effectiveness of the method is validated by comparison of the simulation results with analytical solutions based on transfer matrix method. Different Lamb wave mode interactions with delamination are evaluated, and it is demonstrated that symmetric Lamb wave modemay be insensitive to delamination at certain interfaces of laminates while the antisymmetric mode is more suited for identification of delamination in composite structures.


Introduction
Owing to its superior mechanical properties and light weight, composite materials are finding more and more applications especially in aerospace industries 1 . However, composite structures still run a high risk of suffering from damage due to abrupt impact or growth of fatigue defects, which can result in catastrophic failure during their service life. It is, therefore, essential to develop techniques to inspect integrity and improve safety, reliability, and operational life of structures 2-7 .
Traditionally, nondestructive evaluation NDE techniques, such as C-scan and radiographic inspection, are used to evaluate the integrity and degradation of structures on a periodic basis. Now online structural health monitoring SHM techniques, for example, vibration-based and Lamb-wave-based techniques, are being developed to provide early warning and assure the performance of structures. Among them, Lamb-wave-based techniques for damage detection have received a considerable attention in the past decades due to its ability of long-distance propagation and sensitivity to a variety of defects.
For Lamb waves-based damage detection techniques, understanding wave propagation characteristics in structures is essential for their successful application. Therefore, a number of numerical methods have been applied to analyze propagation of elastic waves, such as finite difference method FDM 8,9 , finite element method FEM 10-12 , boundary element method BEM 13, 14 , finite strip elements FSE 15, 16 , mass-spring lattice model  , and local interaction simulation approach LISA 22,23 . In the literature, two kinds of spectral element method SEM were applied to wave propagation modeling, namely, fast-Fourier-transform-FFT-based SEM and the orthogonal polynomials-based SEM 24-26 . In the FFT-based SEM 24 , a single element is sufficient to model wave propagation in large uniform structures, which is suited for simple 1D and 2D problems 27 . On the other hand, the orthogonal-polynomials-e.g., Legendre and Cheybysev polynomials based SEM 25 is much more suitable for analyzing wave propagation in structures with complex geometry. The formulation of the spectral finite element SFE is similar to FE when assembly of element matrices and solution of equations are considered. Hence, the SEM can be used to handle the wave propagation in structures with complex geometry, and various types of defects can be modeled. The difference between the SEM and FEM comes in that orthogonal polynomials are used as approximation functions in SEM and, therefore, calculation can be more efficient because of the diagonal mass matrix. This method has been successfully applied to many fields, such as fluids, seismology and acoustics 28, 29 . More recently, the SEM was used to simulate wave propagation in structures for damage detection, for example, wave propagation in 1D structures, such as rod and beam 30, 31 . Numerical simulation results of the elastic wave propagation in a composite plate were presented by Kudela et al. 32 . A 2D spectral membrane finite element-based model, developed by Zak et al. 33 , was used to analyze wave propagation in a cracked isotropic panel. Wave propagation in 2D plate structures using a 3D SEM for damage detection was also discussed by Peng et al. 34 .
A wave propagation analysis in composite structures using 3D SEM for the purpose of damage detection has not been widely reported in the literature so far. Although characteristics of Lamb wave propagation in composite laminates and the damage evaluation using numerical simulation have been rigorously explored for a couple of years, in most of the related work, structures were simplified by either one-or two-dimensional models, resulting in approximate results especially for laminates of complicated layup or with damage. The SEM combines accuracy with flexibility in describing problems with complex geometries, which is highly desirable for modeling of elastic wave propagation. In this paper, multilayered composite laminates are modeled using the Legendre polynomials-based spectral finite element, elastic wave propagation characteristics are analyzed, and wave interaction with delamination is discussed.

Wave Propagation in Composite Plate
Composite laminates are commonly fabricated by stacking unidirectional lamina with a certain layup configuration. After a composite is properly cured, a multilayered structure is formed with all the layers bonded together. can be regarded as a quasihomogeneous orthotropic or transversely isotropic material with the main principal axis parallel to the fibres. Therefore, in numerical calculations, a composite laminate is modeled as a multilayered medium with different elastic and anisotropic material properties 35 . In order to form a numerical model of composite structure, a Cartesian coordinate system is firstly defined by the z-axis normal to the central plane of a composite laminate spanned by the x-axis and the y-axis for modeling of the wave propagation, as shown in Figure 1.
For each layer of the composite laminate, the stress-strain relations for arbitrary direction have the following form 36 : where σ is stress vector, ε is strain vector, and D is the flexibility matrix. In order to study wave propagation, the elastic constants of all the layers must be expressed in the global coordinate system. For those layers, the principal material coordinate system does not coincide with the global coordinate system; this can be achieved by using a transformation matrix method.
In a homogeneous media, the elastic wave propagation is described by the governing equation 37 : where u is the displacement vector. ρ is mass density and f is external force vector.

Formulation of 3D Spectral Element Method
For the Legendre polynomials-based 3D spectral finite element, it requires that the domain Ω in three dimensions is decomposed into a number of nonoverlapping hexahedrons, Ω e , as in the conventional FE method. In SEM, the equations of wave propagation is 37 where the Ω denotes the physical region of interest and w is an arbitrary test vector.
In SEM, the nodes are defined into two steps: 1 each element in its physical domain is mapped to a reference domain Λ −1, 1 3 using an invertible local mapping f; a set of basis functions consisting of Legendre polynomials of degree N are introduced; 2 a set of 4 Mathematical Problems in Engineering nodes are defined as ξ i ∈ −1, 1 , i ∈ 1, . . . , N 1. These Gauss-Lobatto-Legendre GLL points are the N 1 roots of 29 where P N ξ is the derivative of the Legendre polynomial of degree N. The ξ i are different from the classical FE method in which the nodes are uniformly spaced. As an example, a 108-node spectral element in the local coordinate is shown in Figure 2. The Lagrange interpolation function, u e N , supported by the GLL points can be expressed as where Ψ mnr is defined as the orthogonal shape functions in 3-D h m ξ denotes the mth 1D Lagrange interpolation at the N 1 GLL points defined above. The property of h m ξ is where δ mn denotes the Kronecker delta and n i , i 1, 2, 3, are the numbers of GLL points in each direction in the local coordinate.
Therefore, the element matrices, M e mass matrix , K e stiffness matrix , and F e force vectors , are calculated numerically on the reference coordinate:

3.5
where ρ is the mass density, D e is termed material stiffness matrix, and P is a distributed load. Ψ e are the shape functions based on the Legendre polynomials. The matrix B e is the strain-displacement matrix calculated by where L denotes a differential operator matrix: Mathematical Problems in Engineering J e is the Jacobian matrix associated with the mapping f from the physical domain Ω e to the reference domain Λ: The weights ω i are defined by Therefore, wave equation 3.1 can be rewritten into matrix form and wave propagation modeling is transformed to an ordinary differential equation in time. Let U denotes the global vector of unknown displacement in the medium. Then the ordinary differential equation can be written as where M denotes the global mass matrix, K is the global stiffness matrix, and F is the vector of time-dependent excitation force. In this study, the differential equation 3.10 is solved using a central difference time integration scheme. Initial conditions of displacement and velocity are U 0 andU 0 at the time t 0, and the central difference time integration scheme is implemented as where the symbol t denotes time and Δt denotes the time step of integration. In the central difference time integration scheme, the stable condition is Δt ≤ Δt cr L/c, where L is the minimum distance between two adjacent nodes and c is the wave speed in elastic medium.
In comparison with FE, less computation effort is required for SFE because of the choice of Lagrange interpolation supported on the GLL points in conjunction with the GLL integration rule. The efficiency of SEM was demonstrated using two 3D models based on SFEs and FEs with the same degrees of freedom, and a reduction of about 65% in CPU time for calculation is used in comparison with the FEM.

Numerical Calculation
Lamb wave propagation in 8-ply T300/F593 composite laminates is analyzed in this study. The material properties of unidirectional lamina are listed in Table 1. Three types of composite laminates, namely, unidirectional 0 8 , symmetric cross-ply 0 2 /90 2 s , and quasi-isotropic 45/45/0/90 s laminates are investigated. All the laminates have the same geometric configurations of 500 mm × 500 mm × 1.72 mm, as shown in  was meshed using the three-dimensional spectral finite elements with 108-node shown in Figure 2 .
In the literature, piezoelectric PZT wafer is one of the frequently applied transducers to excite and capture Lamb waves propagating in structures for damage detection. The PZT transducers can be modeled by adding force to the incident points or using electromechanical coupling for elastic wave modeling in structures 38, 39 . In the present study, in order to obtain relatively simple Lamb wave mode, two forces perpendicular to the plane of the plate are applied at point A in Figure 3 on the upper and the lower surfaces of the plate, respectively. When the two forces are in phase, antisymmetric modes are activated, and when the two forces are out of phase, symmetric modes are excited. The excitation force is a 5-cycle sinusoidal signal modulated by Hanning window with a center frequency of 100 kHz and absolute maximum magnitude is 1N, as shown in Figure 4. The displacement responses at point B and C on the upper surface of the composite laminate are used to investigate wave propagation characteristics.

Wave Propagation in a Multilayered Composite Plate
As aforementioned, three types of the composite laminates are analyzed in this study. plotted in Figure 5 a . It can be observed that the wave front of the A 0 mode has ellipselike shape, and the group velocity in the fibre direction is greater than that in the direction perpendicular to the fibre, indicating that group velocity of this type of wave mode in the composite depends on the orientation of wave propagation. Generally, group velocity in composite laminate is a function of the direction of wave propagation and direction of the fibres. The group velocity in composite laminates can be calculated analytically using the transfer matrix method In the appendix . For such an 8-ply unidirectional laminate, the analytical group velocity is plotted in a polar coordinate, as shown in Figure 5 b . It presents similar feature as Figure 5 a . Group velocities in the direction along the fibers x-axis and perpendicular to the fibres y-axis are further calculated from displacement responses of SEM simulations at the point B and C, as shown in Figure 6. Based on the peak of the received responses, time of flight ToF from point A to B can be defined. The calculated group velocity c g of the A 0 mode of Lamb waves propagating in the x direction is 1794 m/s, which is 2.6% smaller than the analytical value. In a similar way, the group velocity of the A 0 mode in the y-direction is 1319 m/s and the one calculated analytically is 1245 m/s, which gives a relative error of 5.9%. It can be concluded that the simulation results of SEM model, the proposed model, are in good agreement with the analytical results thus validating the effectiveness of the model. Under the two excitations of out-of-phase forces, the fundamental symmetric wave modes, the S 0 mode and the SH 0 mode, are excited. According to the simulation results, the displacement components in the x-direction and in the y-direction are dominated. Hence, only those two displacement components are plotted, as shown in Figures 7 a and 7 b . It can be seen that S 0 mode and the SH 0 mode are excited simultaneously, which is because the S 0 mode and the SH 0 mode are coupled in multilayered composite laminate. Analytical group velocities of the S 0 mode and the SH 0 mode are also calculated, as shown in Figure 7 c , and the simulation results of SEM modeling are in good agreement of the analytical results.
In case of cross-ply laminates 0 2 /90 2 s , the laminate was meshed to 50 × 50 × 4 spectral finite elements. The symmetric mode and antisymmetric modes are excited using the above-mentioned method. In case of antisymmetricmode, the displacement component in the z-direction is shown in Figure 8. On the other hand, in case of symmetric mode the displacement in the x-direction and the y-direction are plotted in Figure 9. The complexity or 225 • directions because outer lamina are orientated in these directions, which dominates the bending properties related to A 0 mode. In case of angle-ply laminate, the wavelength of the A 0 mode is shorter than that of the S 0 and the SH 0 modes in composite laminate from the simulation results and analytical results. It can be expected that the A 0 mode is possibly more sensitive to small damage than the S 0 mode and the SH 0 mode because of its short wavelength.
It is demonstrated that characteristics of the wave propagation in multilayered composite is complex due to the nature of anisotropic of the material constants and the multilayered configurations, which leads to the group velocity of Lamb waves depending on the laminate layup and the direction of wave propagation. Good agreement between the simulation results based on SEM and the analytical results demonstrates that the proposed model provides an effective tool to investigate the wave propagation in composite structures.

Wave Propagation in a Composite Plate Containing Delamination
The typical damage forms in composite laminate are transverse microcracking, fiber-breakage, and delamination. Typically, the transverse microcracking through the thickness of the ply occurs as the first-ply failure, and then delamination damage follows. The fiber breakage usually happens at the last stage of the failure. However, a catastrophic failure can occur only with the microcracking and delamination damage without the fiber breakage. Delamination is known to happen because of excessive interlamina normal and shear stress at the ply boundaries, which not only causes reduction in stiffness, but also affects the strength and integrity of the structure, leading to failure.
The wave propagation in composite laminates containing a delamination is investigated in this study. The delamination in the laminates is modeled using nodes separation method. Laminate without delamination is initially meshed. At the interface between the adjacent elements where the delamination occurs, nodes that are affected by the delamination are separated, as shown in Figure 12.
Lamb wave modes interaction with delamination is analyzed in quasi-isotropic laminates 45/−45/0/90 s using the proposed 3D SEM model. The size of the delamination is 30 mm × 30 mm in a square shape, as shown in Figure 3. The wave interaction of the symmetric mode and antisymmetricmode with delamination is investigated.
The effects of the delamination at different interfaces in the composite laminate are addressed. Under the symmetric mode, the displacement responses at the point B in the xand the y-directions are plotted in Figure 13. Responses of the intact composite laminate are also provided for comparison. It is evident that the scattered waves from the delamination are not so obvious in comparison with the intact laminate, on the captured responses at the point B, when the delamination is located in different interfaces. The displacement responses at 0.069 ms in the x-and the y-directions of the composite laminates are plotted in Figure 14, when there is a delamination at the interfaces between 3-4 layers. However, when the delamination is at the interface between 4-5 layers, responses of the point B are same as those from intact composite laminate, indicating that the symmetric mode is insensitive to the delamination in the symmetric plane of the plate, attributed to the reason that the layup of the composite laminate is symmetric about the central plane between 4-5 layers. Therefore, the symmetric wave mode related to extensional wave mode can travel through the delamination without any scattering from it.
Under the antisymmetricmode, the displacements at the point B in the z-direction are plotted in Figure 15. It can be seen that the antisymmetricmode is more sensitive to  delamination located at all interfaces of the composite laminate. The displacement responses of the laminate in the z-direction at 0.173 ms are plotted in Figure 16. Under the same excitation frequency, the scattered waves from delamination become clearer in the captured response at the point B under the antisymmetricmode than that under symmetric mode. Hence, the antisymmetricmode can be used to detect smaller delamination than the symmetric mode, since the wavelength of the antisymmetricmode is less than that of the symmetric mode. In addition, according to the simulation results, the symmetric mode may be less sensitive to delamination located at certain interfaces. It can be expected that the antisymmetricmode is more suitable for identification of the delamination.  The scattered waves are very weak compared with the incident wave, resulting in that the reflected wave packet is difficult to identify delamination in such composite materials with a high attenuation ratio. Fortunately, the transmitted wave has been affected a lot, as shown in Figure 16.

Conclusions
A three-dimensional spectral element method is developed to investigate the wave propagation characteristics in composite laminates in the present study. Three types of laminates, namely, unidirectional 0 8 , cross-ply 0 2 /90 2 s , and quasi-isotropic −45/45/0/90 s laminates are modeled using 3D spectral elements, and the laminates are excited using two forces in-phase and out-of-phase to generate the symmetric mode and antisymmetricmode respectively. It is demonstrated that the proposed 3D spectral element method can be efficiently and effectively used to simulate wave propagation in composite laminates. Complexity of wave propagation characteristic is also demonstrated even for a single Lamb wave mode in composite laminates. Finally, interactions between the Lamb wave mode and a delamination are analyzed. It is concluded that symmetric mode of Lamb wave may be insensitive to delamination in certain interfaces in the laminate. And, therefore, it is essential to understand wave propagation characteristics in composite laminate when Lamb-wave-based structural health monitoring strategy is carried out.
The existence of nontrivial solutions for U 1 , U 2 , and U 3 demands the vanishing of the determinant of the matrix K and yields the sixth-degree polynomial equation: A.4 There are six roots of this equation, which correspond to the three sets of mode pairs. For each α q , q 1, 2, . . . , 6, the displacement ratios V q U 2q /U 1q and W q U 3q /U 1q can be expressed as V q K 11 α q K 23 α q − K 13 α q K 12 α q K 13 α q K 22 α q − K 12 α q K 23 α q , W q K 11 α q K 23 α q − K 12 α q K 13 α q K 12 α q K 33 α q − K 23 α q K 13 α q . A.5 The formal solutions for the displacements and stresses in the expanded matrix form where E q e iξα q x 3 , D 1q iξ C 13 C 36 V q C 33 α q W q , D 2q iξ C 55 α q W q C 45 α q V q , D 3q iξ C 45 α q W q C 44 α q V q , q 1, 2, . . . , 6. A.7