Abstract
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 mode may be insensitive to delamination at certain interfaces of laminates while the antisymmetric mode is more suited for identification of delamination in composite structures.
1. 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 (MSLM) [17–21], 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.
2. 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. For analysis of wave propagation, each lamina 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 -axis normal to the central plane of a composite laminate spanned by the -axis and the -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 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 is the displacement vector. is mass density and is external force vector.
3. 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,, as in the conventional FE method. In SEM, the equations of wave propagation is [37] where the denotes the physical region of interest and 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 using an invertible local mapping ; a set of basis functions consisting of Legendre polynomials of degree are introduced; (2) a set of nodes are defined as , . These Gauss-Lobatto-Legendre (GLL) points are the roots of [29] where is the derivative of the Legendre polynomial of degree . The 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,, supported by the GLL points can be expressed as where is defined as the orthogonal shape functions in 3-D denotes the th 1D Lagrange interpolation at the () GLL points defined above. The property of is where denotes the Kronecker delta and , , are the numbers of GLL points in each direction in the local coordinate.
Therefore, the element matrices, (mass matrix), (stiffness matrix), and (force vectors), are calculated numerically on the reference coordinate: where is the mass density, is termed material stiffness matrix, and is a distributed load. are the shape functions based on the Legendre polynomials. The matrix is the strain-displacement matrix calculated by where denotes a differential operator matrix:is the Jacobian matrix associated with the mapping from the physical domain to the reference domain: The weights 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 denotes the global vector of unknown displacement in the medium. Then the ordinary differential equation can be written as where denotes the global mass matrix, is the global stiffness matrix, and 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 and at the time , and the central difference time integration scheme is implemented as where the symbol denotes time and denotes the time step of integration. In the central difference time integration scheme, the stable condition is , where is the minimum distance between two adjacent nodes and 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.
4. 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 , symmetric cross-ply , 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 Figure 3. The composite plate 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 , 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.
4.1. Wave Propagation in a Multilayered Composite Plate
As aforementioned, three types of the composite laminates are analyzed in this study. Firstly, Lamb wave propagation in the unidirectional composite laminate is investigated. The plate is meshed to 50 × 50 × 1 spectral elements. Under the two excitation of in-phase forces, the fundamental antisymmetric mode A0 is excited. According to the simulation results, the displacement component in the -direction is dominant. Therefore, only the displacement component in the -direction is analyzed here. Responses of the laminate at 0.097 ms are plotted in Figure 5(a). It can be observed that the wave front of the A0 mode has ellipse-like 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 (-axis) and perpendicular to the fibres (-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 of the A0 mode of Lamb waves propagating in the direction is 1794 m/s, which is 2.6% smaller than the analytical value. In a similar way, the group velocity of the A0 mode in the -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.
(a)
(b)
(a)
(b)
Under the two excitations of out-of-phase forces, the fundamental symmetric wave modes, the S0 mode and the SH0 mode, are excited. According to the simulation results, the displacement components in the -direction and in the -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 S0 mode and the SH0 mode are excited simultaneously, which is because the S0 mode and the SH0 mode are coupled in multilayered composite laminate. Analytical group velocities of the S0 mode and the SH0 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.
(a)
(b)
(c)
In case of cross-ply laminates , 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 -direction is shown in Figure 8. On the other hand, in case of symmetric mode the displacement in the -direction and the -direction are plotted in Figure 9. The complexity of Lamb wave propagation in composite is also demonstrated, and the effectiveness of the proposed model is validated by comparison of simulation results and the analytical results.
(a)
(b)
(a)
(b)
(c)
In case of angle-ply laminates , the composite laminate is meshed by using 50 × 50 × 8 spectral finite elements. The displacement components in the -direction at 0.086 ms are plotted in Figure 10(a), while displacement components in the -direction and in the -direction at 0.069 ms are plotted in Figures 11(a) and 11(b). The group velocities of analytical solutions for those Lamb modes are also plotted in a polar coordinate for comparison, as shown in Figures 10(b) and 11(c). A good agreement found in both symmetric and antisymmetricmodes. It can be observed that the angular dependence of Lamb nodes in the laminates becomes weaker because of its quasi-isotropic layup. Under the S0 and the SH0 modes, the group velocities are approximately independent of direction of wave propagation, but it still can be discerned that the A0 mode has the maximum in the 45° (or 225°) directions because outer lamina are orientated in these directions, which dominates the bending properties related to A0 mode.
(a)
(b)
(a)
(b)
(c)
In case of angle-ply laminate, the wavelength of the A0 mode is shorter than that of the S0 and the SH0 modes in composite laminate from the simulation results and analytical results. It can be expected that the A0 mode is possibly more sensitive to small damage than the S0 mode and the SH0 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.
4.2. 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.
(a)
(b)
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 - and the -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 - and the -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.
(a)
(b)
(a)
(b)
Under the antisymmetricmode, the displacements at the point B in the -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 -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.
5. 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 , cross-ply , and quasi-isotropic 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.
Appendix
The analytical wave front can be calculated as follows.
Consider wave propagation solutions in the following form: where is the wave number, is the circular frequency, is still an unknown parameter, and and are ratios of the displacement amplitudes of and , respectively. The choice of the solution leads to the three coupled equations that can be written as where the elements of the matrix are The existence of nontrivial solutions for , , and demands the vanishing of the determinant of the matrix and yields the sixth-degree polynomial equation: There are six roots of this equation, which correspond to the three sets of mode pairs. For each , , the displacement ratios and can be expressed as The formal solutions for the displacements and stresses in the expanded matrix form where
Acknowledgments
The authors are grateful for the support received from the National Natural Science Foundation of China (NSFC nos. 11072148 and 11061160491), Research Project of State Key Laboratory of Mechanical System and Vibration (MSV201110), and the National High Technology Research and Development Program of China (no. 2009AA044800).