Effects of Delamination on Guided Waves in a Symmetric Laminated Composite Beam

For successful structural healthmonitoring and structural integrity evaluation of a laminated composite structure, it is important to study the effects of delamination on the propagations of the guidedwaves in a delaminated composite beamby using an accurate and computationally efficient method.Thus, we developed a “frequency-domain” spectral element model for the symmetric composite beams. First-order-shear-deformation-theory (FSDT) based Timoshenko beam theory and Mindlin-Herrmann rod theory are adopted for the flexural (bending) waves and axial (extensional) waves, respectively. A spectral element model is derived from the governing equations of motion by using the variation method in the frequency domain. After validating the accuracy of the proposed spectral element model, the model is used to investigate the effects of delamination on the propagation of guided waves in examples of composite beams.


Introduction
Due to the many inherent advantages of composite materials over metallic materials such as high performance, stiffness, strength, and low weight, composite materials have been increasingly applied in the fields of aerospace, mechanical engineering, civil engineering, and other technologies over the past decades [1].However, composite materials are still susceptible to delamination and other modes of damage.Therefore, there is a growing need for structural health monitoring (SHM) of composite structures.
Recently, guided wave-based SHM techniques have received considerable attention due to the guided wave characteristics of long-range propagation and sensitivity to various modes of damage [2].The smaller wavelength of guided waves is more effective for detecting or monitoring incipient small damages; thus, high frequency guided waves are essential for most SHM techniques.Piezoelectric transducers have been used to excite or sense high frequency guided waves (on the order of MHz).
Knowledge of the effects of damage on guided wave propagation is very useful for successful structural health monitoring, and wave propagation analysis has been extensively reported in the literature.In practice, ultrasonic guided waves in damaged composite structures are very complicated due to anisotropic material properties, multilayered construction, and multiple reflections and mode conversions at material and geometric discontinuities.Thus it is essential to develop an efficient, low-cost, and accurate wave propagation analysis method to investigate ultrasonic guided waves in damaged composite structures.
As described in the literature, a number of computational methods have been used for the analysis of guided wave propagation.These include the finite element method (FEM) [3][4][5][6], strip element method [7], boundary element method [8], local interaction simulation approach (LISA) [9], and "timedomain" spectral element method [10][11][12].Among these approaches, the FEM, which is normally used in the timedomain, is certainly one of the most popular methods for wave propagation analysis due to its particular advantages.To successfully capture high frequency waves, the finite element mesh size must be smaller than the wavelengths involved, which can be very small at high frequencies.This can result in an extremely large finite element model.Hence, the FEM becomes increasingly inaccurate and inefficient as the frequency increases.This problem can be resolved by adopting the "frequency-domain" spectral element method (SEM), wherein frequency-dependent interpolation functions are used to generate an exact dynamic stiffness matrix (called the "spectral element matrix") by which high frequency waves can be successfully captured [13,14].The SEM provides solutions that are often referred to as "exact" due to the use of an exact dynamic stiffness matrix.It also has the advantage of reducing computational effort due to a significant reduction in the number of degrees of freedom (DOFs) by representing a uniform structural member as a single spectral element model, regardless of its dimension.It is worth mentioning that, though the "time-domain" SEM initiated by Patera [11] is also often called by the same name, "SEM, " in the literature, it is completely different from the "frequency-domain" SEM because it belongs to a -version FEM where the higher-order polynomials such as orthogonal Chebyshev polynomials are used as the shape functions.
In the literature, some researchers have used the "frequency-domain" SEM to study the propagations of guided waves in a cracked rod [15], cracked beam [16], delaminated beam [17], debonded damaged reinforced concrete beam [18], cracked composite beam [19], and delaminated composite beam [20,21].For composite beams, Palacz et al. [21] used a spectral element model based on the first-order shear deformation theory (FSDT), wherein the axial displacement and the lateral contraction in the thickness direction (i.e., the thickness contraction) were not considered.Nag et al. [20] considered axial-flexuralshear coupling in their FSDT-based spectral element model but ignored the effects of thickness contraction.The thickness contraction is known to be important for axial (or extensional) waves at high frequencies via the axialthickness contraction coupling, which can be represented by the Mindlin-Herrmann rod theory [22].Thus, it is very important to fully take into account the effects of axial-thickness contraction coupling in the spectral element model in order to successfully investigate the effects of delamination on axial and flexural (bending) waves in a delaminated composite beam at high frequencies.
Thus, the purposes of this study are (1) to derive governing equations of motion for a symmetric composite beam by adopting both Timoshenko beam theory and Mindlin-Herrmann rod theory, (2) to develop a "frequency-domain" spectral element model for a finite composite beam element, and (3) to investigate the effects of delamination on guided waves in some examples of composite beams.

Derivation of the Equations of Motion
Consider a symmetric composite beam as shown in Figure 1, where spatial coordinates  and  represent the axis and thickness directions of the composite beam, respectively.The composite beam has plane motion in the - plane.Axial and flexural displacements are represented by (, , ) and (, , ), respectively.Based on FSDT-based Timoshenko beam theory [1], these displacement fields can be represented by  (, , ) =  0 (, ) +  (, ) ,  (, , ) =  0 (, ) −  (, ) , where (, ) is the rotation of the cross-section of the composite beam and (, ) is the thickness contraction.The subscript "0" denotes a quantity at the midplane of the composite beam.The equations of motion for a composite beam element of length  can be derived from Hamilton's principle [1]: where , , and  are the kinetic energy, the potential energy, and the virtual work done by external forces, respectively.By using the displacement fields given by (1) and the constitutive relations for the composite laminates provided in Appendix A, the kinetic and strain energies can be given by where where  is the width,  ()  is the mass density of the th laminate of the composite beam, and material properties  *  and  *  are defined in Appendix A. Finally the virtual work  can be represented as where   (),   (), and   () ( = 1, 2) are, respectively, the resultant axial forces, bending moments, and the transverse shear forces applied at the boundaries  = 0 and .  (, ),   (, ), and   (, ) are externally applied distributed forces.
By substituting (3) and ( 5) into (2) and applying integration by parts, the governing equations of motion and associated natural boundary conditions can be derived simultaneously.The equations of motion are obtained in the following form: where w (, ) = { 0 (, ) ,  0 (, ) ,  (, ) ,  (, )}  , and D 0 , D 1 , D 2 , and  are four-by-four matrices defined in Appendix B. The natural boundary conditions are also obtained in the form where D 1 is the four-by-four matrix defined in Appendix B and

Spectral Element Model for a Composite Beam Element.
To formulate a spectral element model for a composite beam element of length , the equations of motion and the associated boundary conditions must first be transformed in the frequency domain.To that end, the time-domain displacement fields and forcing terms are all represented in spectral forms by applying the discrete Fourier transform theory [23] as follows: where  = √ −1 is the imaginary unit and the quantities with overbars are the spectral (Fourier) components of the corresponding time-domain quantities./2 is the number of spectral components up to the Nyquist frequency to be considered in the FFT-based spectral analysis.By using (10), the time-domain equations of motion given by ( 6) can be transformed in the frequency domain as follows: Similarly the boundary conditions can be also transformed into the frequency domain as In ( 11) and ( 12) and also in the following derivations, the subscript  is omitted for the sake of brevity.Consider the homogeneous frequency-domain equations of motion that can be reduced from (11) by setting p(; ) = 0 as follows: The general solutions of ( 13) are assumed in the form where  is a constant and  is the wavenumber.The vector r is defined by Substituting ( 14) into (13) yields the following eigenvalue problem: The wavenumber-frequency dispersion equation can be derived from (16) in the form By using eight wavenumbers   ( = 1, 2, 3, . . ., 8) computed from (17), the general solutions of ( 13) can be written as where The th eigenvectors r  = {1,   ,   ,   }  in (20) can be readily computed from ( 16) by specifying that  =   .For a finite composite beam element of length  as shown in Figure 2, the nodal degrees of freedom in the frequency domain (called spectral nodal DOFs) are defined as where By applying (18) to (22), the spectral nodal DOFs vector d() can be represented as Using (25), the constant vector a can be removed from (18) to write the general solutions in terms of d() as follows: where To formulate the spectral element model by using the variational method, the weak form of ( 6) is obtained from its weighted-integral statement, by integrating by parts and then applying the boundary conditions, as follows: where p (; ) = {  (; ) ,   (; ) ,   (; ) , 0}  , By substituting ( 26) into (28), the spectral element equation can be obtained in the following form: where and S() is the spectral element matrix defined in Appendix C.

Spectral Element Model for a Delaminated Composite
Beam Element. Figure 3 shows a composite beam element that contains a delamination inside.Following the approaches introduced in [20,21], a spectral element model for the delaminated composite beam element can be obtained by assembling spectral element equations for the four segments shown in Figure 3.The spectral element equations for four segments can be written in partitioned form as where and the superscripts  = 1, 2, , and  represent the left, right, upper, and lower segments in Figure 3, respectively.The kinematic connectivity conditions at the vertical interfaces between four segments can be represented as where By using (34), the spectral element equations for four segments represented by (32) can be assembled in the following form: S (1)   12 0 0 S (1)  21 S (1)  22 + S (+) + S (2)  11 S (2)   12 0 0 S (2)   21 (2) 1 where It is very often to use piezoelectric (PZT) transducers to excite or measure the high frequency dynamic responses of a structure.To that end, the PZT transducers are usually bonded on a surface of the structure.The spectral element model for the composite beam element whose upper surface is bonded with a PZT transducer is available from [24].

Assembly and Spectral Element Analysis.
The assembly of the spectral elements for a composite beam structure is quite straightforward because it can be readily accomplished with the same approach that is commonly used in a standard FEM.The spectral element analysis methods for both frequencydomain solutions (i.e., frequency response functions and natural frequencies) and time-domain solutions (i.e., time histories of dynamic responses and waves) have been well developed in the literature (see [13,14]).

Numerical Results and Discussion
The

Validation of Spectral Element Model.
To evaluate the accuracy of the present spectral element model, a cantilevered cross-ply symmetric composite beam [0/90]  without delamination was considered as the first example problem.The examples of composite beam has thickness of 5 mm and length of 600 mm.The natural frequencies and frequency response functions (FRFs) obtained using the present spectral element model and the ANSYS commercial finite element (FE) analysis package [26] are compared in Table 1 and Figure 4, respectively.A one-element model was used to obtain the SEM results.The number of finite elements used in the ANSYS FE analysis was increased step by step to investigate the convergence of the ANSYS results.For the ANSYS FE analysis, eight-node quadratic plane-stress elements were used.Table 1 and Figure 4 show that the ANSYS results converged to the SEM results as the number of finite elements (i.e., the number of DOFs) was increased.Small deviations are probably due to differences in the kinematic assumptions and constitutive relations adopted in the present SEM and ANSYS formulations.
As the second example problem, a cantilevered crossply symmetric composite beam [0/90]  with a delamination was considered.It was assumed that the length () of the composite beam is 600 mm.The delamination is located in the middle ( = 300 mm) of the composite beam, as shown in Figure 5, and it has length   = 5 mm.The natural frequencies obtained using the present spectral element model and ANSYS are compared in Table 2.It is obvious from Table 2 that the ANSYS results converged to the SEM results as the number of finite elements was increased.These observations made from two example problems certainly validate the high accuracy of the present spectral element model.Note: ( ) is the number of degrees of freedom used in the analysis.

Effects of the Delamination Length and Location.
To investigate the effects of the delamination length and location on the dynamic characteristics and waves of a composite beam, we reconsidered the cantilevered cross-ply symmetric composite beam [0/90]  shown in Figure 5.The composite beam has the length () of 1400 mm and a delamination of length   is located in the middle ( = 700 mm) of the composite beam.The location of the delamination in the thickness direction is represented by  (the ratio of the thickness of the upper part of the delamination to the composite beam thickness).It was assumed that examples of composite beams considered in this paper are all initially at rest.Table 4 shows the effects of delamination length on the natural frequencies of the composite beam.Four delamination lengths (  = 0, 5, 50, and 100 mm) were considered and it was assumed that they were located at  = 0.25.Table 4 certainly shows that the natural frequencies tend to decrease as the delamination length gets increased.
To investigate the effects of the delamination length on the guided waves when  = 0.25, the flexural (bending) and axial (extensional) waves that propagate forward and backward through the upper and lower parts of the delamination were observed at  = 500 mm for four values of delamination length (  ).The guided waves were excited by applying an input point force in the form of a Morlet wavelet signal with a 100 kHz center frequency and having five cycles at the free end ( = 1400 mm) of the composite beam.The results are compared in Figure 6.The center frequency of the excitation force was chosen to be 100 kHz.Accordingly, wavelengths of the flexural wave  0 (, ) and the axial wave  0 (, ) are 18.6 mm and 74.4 mm, respectively.The axial displacement and flexural displacement are not coupled for a symmetric composite beam.Only the flexural waves are initially excited by the transverse loading.This is why only flexural wave packets were observed at  = 412 s when no delamination existed (i.e.,   = 0).On the other hand, Figure 6 shows that axial wave packets can be observed when delamination exists (i.e.,   ̸ = 0).Axial wave packets are, in fact, generated by delamination-induced axial-bending coupling.From Figure 6, we note the followings.(1) When the delamination length is larger than the wavelength of a flexural wave, a part of the flexural wave that propagates across the delamination is repeatedly reflected backward and forward within the delamination zone.This seems to result in tail waves that appear after the flexural wave packet at  = 412 s.On the other hand, when the delamination length is smaller than the wavelength of a flexural wave, the aforementioned wave reflections seem to be negligible.(2) The axial wave packets that appear before the flexural wave packets at  = 412 s are generated in part by the mode conversion that may occur via delamination-induced axial-bending coupling when a flexural wave arrives at the delamination and in part by multiple wave reflections within the delamination zone.
Details of the flexural waves that propagate through the upper and lower parts of a delamination located at  = 0.25 are shown in Figure 7 for two cases: (1) when the delamination length is 5 mm and (2) when the delamination length is 50 mm.When the delamination lengthis smaller than the wavelength of the flexural wave (i.e., 18.6 mm), wave reflections within the delamination zone are so weak that they are not visible, as shown in Figure 7(a).On the other hand, Figure 7(b) shows that there are strong multiple wave reflections within the delamination zone when the delamination length is larger than the wavelength of the flexural wave.Comparing the wave reflections in the upper part and lower part of the delamination shown in Figure 7(b), we note that the wave reflections seem to be stronger in the thinner (upper) part than in the thicker (lower) part of the delamination.
To investigate the effects of the delamination location in the thickness direction on the guided waves, we assumed that 5 mm long delaminations are located at  = 0.25, 0.5, and 0.75.We note that  is defined as the ratio of the thickness of the upper part of the delamination to the total thickness of the whole composite beam.For example,  = 0.5 implies that a delamination is located at the midplane of the composite beam.The guided waves to be observed at  = 0.5 m were simulated for three different delamination locations and the results are compared in Figure 8.The large wave packets at  = 412 s are the flexural waves that start from the excitation point at the free-end tip and propagate directly to the observation point at  = 0.5 m. Figure 8 shows that the flexural wave packets are very similar in both amplitude and phase, regardless of the delamination location.The wave packets observed at  = 350 s are axial waves that are converted from flexural waves via the delaminationinduced axial-bending coupling.Axial wave packets can be observed (see Figure 8) when the delamination is located at  = 0.25 and 0.75 but not when the delamination is located at  = 0.5.This is because the delamination-induced axial-bending coupling diminishes when  = 0.5 due to the material and geometric symmetry of the composite beam under consideration.From Figure 8, we also note that axial wave packets when the delamination is located at  = 0.25 and at  = 0.75 are out of phase with each other.This is because when a symmetric composite beam is subjected to an axial force, the delamination at  = 0.25 tends to deform the composite beam downward, while a delamination at  = 0.75 tends to deform the composite beam upward, and vice versa.3 shows that, regardless of which part of the delamination a flexural wave propagates through, the group velocity of the flexural wave tends to decrease as the fly orientation angle  is increased from 0 ∘ to 60 ∘ .Thus, Figure 9(a) shows that the first arrival of the flexural waves is delayed longer at the larger value of .Both the wavenumber dispersion curves and the group velocity dispersion curves shown in Figures 9(b) and 9(c), respectively, demonstrate that the flexural waves at the larger value of  are more dispersive in the frequency range of excitation.As a result, Figure 9(a) shows that the widths of flexural wave packets at the larger value of  are larger than those at the smaller value of .Compared with the flexural waves observed when no delamination exists in the symmetric composite beam (blue solid lines), the flexural waves observed when a delamination exists (red dotted lines) are found to have tail waves that appear after the main flexural wave packets.The tail waves are generated via the multiple wave reflections within the delaminated zone, as discussed previously.Table 3 shows that  the difference between the group velocity of a flexural wave that propagates through the upper part of the delaminated zone and through the lower part of the delaminated zone becomes larger at a larger value of .Faster waves that travel through the upper part and slower waves that travel through the lower part of the delaminated zone combine to form the flexural waves observed at  = 500 mm.This is why more complex tail waves appear at larger value of .

Conclusions
In this study, a frequency-domain spectral element model was developed for symmetric composite beams.The governing equations of motion were derived by adopting Timoshenko beam theory and Mindlin-Herrmann rod theory for flexural waves and axial waves, respectively.By using a spectral element model formulated using the variation method, the following phenomena were numerically investigated.
(1) Delamination-induced axial-bending coupling can generate axial waves even in a symmetric composite beam.The axial waves in a delaminated symmetric composite beam are generated by mode conversion via the delamination-induced axial-bending coupling.
(2) When the delamination length is larger than the wavelength of a flexural wave that is propagating through a delamination zone, the multiple reflections of the flexural wave within the delamination zone were found to be more significant.
(3) The wave reflections are in general more severe in the thinner (upper) part of the delamination zone than in the thicker (lower) part.
(4) When a transverse loading is applied to a symmetric composite beam, no matter where the delamination is located in the thickness direction, the flexural waves are very similar in both amplitude and phase, which is not true for the axial waves.
(5) For a composite beam with a lay-up of [/90]  , the flexural waves become more dispersive when the fly orientation angle  is increased, with the increase in their group velocities.This makes wave propagation in a delaminated composite beam more complicated.

Figure 3 :
Figure 3: Finite composite beam element with a delamination.

Figure 4 :
Figure 4: Comparison of the frequency response functions (FRFs) obtained by SEM and ANSYS for the composite beam shown in Figure 5:  = the number of finite elements used in the analysis.

Figure 6 :
Figure 6: Flexural and axial waves observed at  = 0.5 m versus the delamination length   .

4. 3 .
Effects of the Lay-Up of Composite Laminates.

8 (Figure 7 :
Figure 7: Flexural waves propagating through the upper and lower parts of a delamination that is located at  = 0.25.

Figure 8 : 5 Figure 9 :
Figure 8: Effects of delamination location () in the thickness direction on the guided waves in the symmetric composite beam shown in Figure 5.

Table 1 :
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the composite beam without delamination.

Table 2 :
Comparison of the natural frequencies (Hz) obtained by SEM and ANSYS for the delaminated composite beam shown in Figure5(  = 5 mm,  = 0.25).

Table 3 :
Group velocities (m/s) of flexural waves versus the fly orientation angle  for the lay-up [/90]  .