Wavelet based spectral finite element modelling and detection of de-lamination in composite beams

M Mitra, S Gopalakrishnan


In this paper, a model for composite beam with embedded de-lamination is developed using the wavelet based spectral finite element (WSFE) method particularly for damage detection using wave propagation analysis. The simulated responses are used as surrogate experimental results for the inverse problem of detection of damage using wavelet filtering. The WSFE technique is very similar to the fast fourier transform (FFT) based spectral finite element (FSFE) except that it uses compactly supported Daubechies scaling function approximation in time. Unlike FSFE formulation with periodicity assumption, the wavelet-based method allows imposition of initial values and thus is free from wrap around problems. This helps in analysis of finite length undamped structures, where the FSFE method fails to simulate accurate response. First, numerical experiments are performed to study the effect of de-lamination on the wave propagation characteristics. The responses are simulated for different de-lamination configurations for both broad-band and narrow-band excitations. Next, simulated responses are used for damage detection using wavelet analysis.


1. Introduction

Composites have several favourable properties, which are encouraging their use in different parts of an aircraft. However, the behaviour of composites under damage is very complicated and not well understood due to the anisotropic nature of the material. The transverse tensile and inter-laminar shear strengths of composites are much lower than the in-plane properties. This makes such structures very much prone to defects like matrix cracking, fibre fracture, fibre de-bonding, de-lamination/inter-laminar de-bonding, of which de-lamination is most common, vulnerable and can grow, thus reducing the life of the structure.

Wave propagation analysis is extensively used for structural health monitoring studies. Proper use of these techniques first requires good knowledge of the effects of damage on the wave characteristics. This needs accurate and computationally efficient modelling of the damaged composite structures. Wave propagation deals with loading of very high frequency content and finite element (FE) formulation for such problems is computationally prohibitive as it requires large system size to capture all the higher modes. These problems are usually solved in frequency domain using Fourier methods, which can in principle achieve high accuracy in numerical differentiation. One such method is FFT based spectral finite element method (FSFEM; Doyle 1999). In FSFEM, first the governing PDEs are transformed to ODEs in spatial dimension using FFT in time. These ODEs are then usually solved exactly, which are used as interpolating functions for FSFE formulation. This results in exact mass distribution and hence, in the absence of any discontinuity, one element is sufficient to handle a one-dimensional structure of any length and this reduces the system size substantially. Apart from this, another advantage of the spectral element method is that the wave characteristics can be extracted directly from such formulation.

The main drawback of FSFEM is that it cannot handle waveguides of short lengths. This is because the required assumption of periodicity results in wrap around problems, which totally distorts the response. It is in such cases, compactly supported wavelets, which have localized basis functions can be efficiently used for waveguides of short lengths. The wavelet based spectral finite element method (WSFEM; Mitra & Gopalakrishnan 2005a, in press) follows an approach very similar to FSFEM, except that Daubechies scaling functions are used for approximation in time for reduction of PDEs to ODEs. The approach removes the problem associated with ‘wrap around’ in FSFEM and thus requires a smaller time window for the same problem. Further, as a consequence, WSFEM can be used for finite length undamped structures, where FSFEM does not work well.

In this paper, the WSFE model of a composite beam with embedded de-lamination is formulated in a similar way to FSFE modelling presented by Nag et al. (2003). The technique used to model a through width de-lamination subdivides the beam into base-laminates and sub-laminates along the line of de-lamination. The base-laminates and sub-laminates are treated as structural waveguides and kinematics are enforced along the connecting line(s). These waveguides are modelled as Timoshenko beams with elastic and inertial coupling and the corresponding spectral elements have three degrees of freedom, namely axial, transverse and shear displacements at each node. The internal spectral elements in the region of de-lamination are assembled assuming constant cross-sectional rotation and equilibrium at the interfaces between the base-laminates and sub-laminates. Finally, the redundant internal spectral element nodes are condensed out to form two-noded spectral elements with embedded de-lamination.

Recently, wavelet analysis of wave propagation for damage detection has been widely reported in the literature (Yuan et al. 2003; Li et al. 2004; Sohn et al. 2004). In the present work, Daubechies wavelet based post processing of measured responses of damaged beams is done for simultaneous de-noising and arrival time based damage detection. Here, the simulated responses of a de-laminated beam modelled using WSFEM are used as substitutes for experimental results. Further, artificial white noise is added to the simulated response. The main aim of wavelet analysis here is to filter out the unnecessary frequency components from the response for better extraction of arrival time which otherwise cannot be obtained directly from the response. De-noising is done by universal soft thresholding at the wavelet levels mainly corresponding to higher frequencies (Donoho 1995).

The paper is organized as follows. In §2, a brief overview of the orthonormal bases of compactly supported Daubechies wavelets is presented. Next, the WSFE formulation for coupled composite Timoshenko beam with axial, transverse and shear degrees of freedom is described very concisely. This is followed by a section on modelling of embedded de-lamination in the composite beam. In §4, the wavelet based damage detection and de-noising techniques are explained. Finally, in the section on numerical experiments, responses simulated for composite beams with different de-lamination layout due to broad-band and narrow-band excitations are presented. Numerical examples of wavelet based damage detection and de-noising are also presented here. The paper ends with some important conclusions.

2. Daubechies compactly supported wavelets

In this section, a concise review of the orthogonal basis of Daubechies wavelets (Daubechis 1992) is provided. Wavelets, ψj,k(t) form the compactly supported orthonormal basis for L2(R). The wavelets and associated scaling functions φj,k(t) are obtained by translation and dilation of single functions ψ(t) and φ(t), respectively.Embedded Image(2.1)Embedded Image(2.2)The scaling functions φ(t) are derived from the dilation or scaling equation,Embedded Image(2.3)and the wavelet function ψ(t) is obtained asEmbedded Image(2.4)ak are the filter coefficients and they are fixed for a specific wavelet or scaling function basis. For compactly supported wavelets only a finite number of ak are non-zero.

The filter coefficients ak are derived by imposing certain constraints on the scaling functions which are as follows:

  1. The area under scaling function is normalized to one,Embedded Image(2.5)

  2. The scaling function φ(t) and its translates are orthonormal,Embedded Image(2.6)and

  3. wavelet function ψ(t) has M vanishing moments,Embedded Image(2.7)

The number of vanishing moments M denotes the order N of the Daubechies wavelet, where N=2M.

Let Pj(f)(t) be the approximation of a function f(t) in L2(R) using φj,k(t) as basis, at a certain level (resolution) j, thenEmbedded Image(2.8)where cj,k are the approximation coefficients. Let Qj(f)(t) be the approximation of the function using ψj,k(t) as a basis, at the same level j.Embedded Image(2.9)where dj,k are the detail coefficients. The approximation Pj+1(f)(t) to the next finer level of resolution j+1 is given byEmbedded Image(2.10)This forms the basis of multi-resolution analysis associated with wavelet approximation.

3. Wavelet based spectral finite element formulation

As mentioned earlier, the de-laminated beam is modelled by splitting the damaged region into sub and base laminates. Each of these laminates are modelled as Timoshenko beams considering elastic and inertial coupling. The details of WSFE modelling of a higher order coupled composite beam with axial, transverse, shear and contractional degrees of freedom are given by Mitra & Gopalakrishnan (2005b). Also, FSFE formulation for a similar composite beam has been presented by Mahapatra & Gopalakrishnan (2003). Here, the WSFE formulation of a composite Timoshenko beam is described very briefly for completeness.

(a) Reduction of wave equations to ODEs

The governing differential wave equations for a composite Timoshenko beam are obtained from those derived by Mahapatra & Gopalakrishnan (2003) for a higher order composite beam.Embedded Image(3.1)Embedded Image(3.2)Embedded Image(3.3)u(x, t), w(x, t) and ϕ(x, t) are the axial, flexural and shear displacements, respectively (see figure 1a).

Figure 1

(a) Beam cross-section and the displacements. (b) Composite beam element with nodal displacements and forces.

The stiffness coefficients are functions of ply properties, orientations, etc., and are obtained by integrating over the beam cross-section asEmbedded Image(3.4)where Embedded Image are the transformed stiffness coefficients of each ply. zi and zi+1 are the thickness coordinates of the ith layer and b is the corresponding width. Similarly, the inertial constants are obtained asEmbedded Image(3.5)where ρ is the mass density. The force boundary conditions associated with the governing differential equations areEmbedded Image(3.6)Embedded Image(3.7)Embedded Image(3.8)where P(x, t), V(x, t) and M(x, t) are the axial, transverse forces and moment, respectively.

The first step of formulation of WSFE is the reduction of the governing differential wave equations (equations (3.1)–(3.3)) to ODEs using Daubechies scaling functions for approximation in time. Let u(x, t) be discretized at n points in the time window [0tf]. Let τ=0,1, …, n−1 be the sampling points, thenEmbedded Image(3.9)where Δt is the time interval between two sampling points. The function u(x, t) can be approximated by scaling function φ(τ) at an arbitrary scale asEmbedded Image(3.10)where uk(x) (referred to as uk hereafter) are the approximation coefficient at a certain spatial dimension x. The other displacements w(x, t), ϕ(x, t) can be transformed similarly and equation (3.1) can be written asEmbedded Image(3.11)Taking inner product on both sides of equation (3.11) with the translates of scaling functions φ(τj), where j=0,1,2, …, n−1 and using their orthogonal properties, we get n simultaneous ODEs as,Embedded Image(3.12)where N is the order of the Daubechies wavelet and Embedded Image are the connection coefficients defined asEmbedded Image(3.13)Similarly, for first-order derivative, Embedded Image are defined asEmbedded Image(3.14)For compactly supported wavelets, Embedded Image, Embedded Image are non-zero only in the interval Embedded Image to Embedded Image. The details for evaluation of connection coefficients for different orders of derivative are given by Beylkin (1992).

It can be observed from the ODEs given by equation (3.12) that certain coefficients uj near the vicinity of the boundaries (j=0 and j=n−1) lie outside the time window [0tf] defined by j=0,1, …, n−1. These coefficients must be treated properly for finite domain analysis. Several approaches like capacitance matrix methods and penalty function methods for treating boundaries are reported in the literature. In this paper, a wavelet based extrapolation scheme (Williams & Amaratunga 1997) is implemented for solution of boundary value problems. This approach allows treatment of finite length data and uses a polynomial to extrapolate the coefficients lying outside the finite domain either from interior coefficients or initial/boundary values. The method is particularly suitable for approximation in time because of the ease of imposing initial values. The above method converts the ODEs given by equation (3.12) to a set of coupled ODEs given asEmbedded Image(3.15)where Γ1 is the first-order connection coefficient matrix obtained after using the wavelet extrapolation technique. It should be mentioned here that though the connection coefficients matrix, Γ2, for second-order derivatives can be obtained independently, here it is written as [Γ1]2 as it helps to impose the initial conditions (Mitra & Gopalakrishnan 2005a). These coupled ODEs are similarly decoupled using eigenvalue analysisEmbedded Image(3.16)where Π is the diagonal eigenvalue matrix and Φ is the eigenvectors matrix of Γ1. Let the eigenvalues be j, then the decoupled ODEs corresponding to equation (3.15) areEmbedded Image(3.17)where Embedded Image and similarly other transformed displacements areEmbedded Image(3.18)Following exactly the similar steps, the final transformed form of the equations (3.2) and (3.3) areEmbedded Image(3.19)Embedded Image(3.20)Similarly, the transformed form of the force boundary conditions given by equations (3.6)–(3.8) areEmbedded Image(3.21)Embedded Image(3.22)Embedded Image(3.23)where Pj and similarly Vj, Mj are the transformed P(x, t) and V(x, t), M(x, t), respectively.

(b) Spectral finite element formulation

The degrees of freedom associated with the element formulation is shown in figure 1b. The element has three degrees of freedom per node, which are Embedded Image, Embedded Image and Embedded Image. From the previous sections, we get a set of ODEs (equations (3.17), (3.19) and (3.20)) for composite beams with axial, transverse and shear modes, in a transformed wavelet domain. These equations are required to be solved for Embedded Image, Embedded Image, Embedded Image and the actual solutions u(x, t), w(x, t), ϕ(x, t) are obtained using inverse wavelet transform. For finite length data, the wavelet transform and its inverse can be obtained using a transformation matrix (Williams & Amaratunga 1994).

It can be seen that the transformed ODEs have a form which is similar to that in FSFEM (Doyle 1999) and thus, WSFE can be formulated following the same method as for FSFE formulation. In this section, the subscript j is dropped hereafter for simplified notations and all the following equations are valid for j=0,1, …, n−1.

The exact interpolating functions for an element of length L, obtained by solving equations (3.17), (3.19) and (3.20), respectively, areEmbedded Image(3.24)where [Θ] is a diagonal matrix with the diagonal terms Embedded Image and [R] is a 3×6 amplitude ratio matrix for each set of k1, k2 and k3.Embedded Image(3.25)k1, k2 and k3 are obtained by substituting equation (3.24) in equations (3.17), (3.19) and (3.20) and solving the characteristic equation. The characteristic equation is obtained by equating the determinant of the 3×3 companion matrix to zero. The corresponding [R] is obtained using singular value decomposition of the matrix. This method of determining wavenumbers and corresponding amplitude ratios was developed to formulate FSFE for graded beam with Poisson's contraction by Chakraborty & Gopalakrishnan (2004). k1, k2 and k3 corresponds to the three modes, i.e. axial, transverse and shear, respectively, and as explained in Mitra & Gopalakrishnan (in press), these are the wavenumbers but only up to a certain fraction of Nyquist frequency.

Here, {a}={A, B, C, D, E, F} are the unknown coefficients to be determined from transformed nodal displacements Embedded Image, where Embedded Image and Embedded Image, Embedded Image, Embedded Image and Embedded Image, Embedded Image, Embedded Image, (see figure 1b for the details of degree of freedom the element can support). Thus, we can relate the nodal displacements and unknown coefficients asEmbedded Image(3.26)From the forced boundary conditions, (equations (3.21)–(3.23)), nodal forces and unknown coefficients can be related asEmbedded Image(3.27)where Embedded Image and Embedded Image, Embedded Image, Embedded Image and Embedded Image, Embedded Image, Embedded Image (see figure 1b). From equations (3.26) and (3.27), we can obtain a relation between transformed nodal forces and displacements similar to conventional FEEmbedded Image(3.28)where Embedded Image is the exact elemental dynamic stiffness matrix. After the constants {a} are known from the above equations, they can be substituted back to equation (3.24) to obtain the transformed displacements Embedded Image at any given x.

(c) Modelling of de-lamination

Here, the modelling of embedded de-lamination is done according to the method presented by Nag et al. (2003). Figure 2a, shows the de-laminated beam broken down into two base-laminates, i.e. laminates 1, 2 on either sides of the de-lamination and two sub-laminates, i.e. laminates 3, 4 on its bottom and top, respectively. Each of these laminates are considered as individual waveguides and modelled as a coupled composite Timoshenko beam using WSFEM described in the last sub-section. These internal spectral elements with the nodes are shown in figure 2b and are assembled assuming constant cross-sectional slope/rotation and force equilibrium at their interfaces. The first assumption of constant and continuous slope at the two interfaces between the sub-laminates and base-laminates gives the following kinematic relation:Embedded Image(3.29)where Embedded Image represents the nodal displacement vector of the ith node in figure 2b. Thus, all the nodal displacements of the sub-laminate elements 3 and 4 are represented in terms of nodal displacements of base-laminate elements 1 and 2. S1 and S2 are the 3×3 transformation matrices in terms of top and bottom sub-laminate thicknesses 2h1 and 2h2 (see figure 2a) and are given asEmbedded Image(3.30)Next using force equilibrium at the two interfaces of sub and base-laminates, the following relations can be derived:Embedded Image(3.31)where Embedded Image represents the nodal displacement vector of the ith node (see figure 2b). The elemental equilibrium equation given by equation (3.28) for jth element with p and q nodes can be rewritten using 3×3 sub-matrices asEmbedded Image(3.32)Using the transformation equations, equations (3.29) and (3.32), the four internal spectral elements are assembled. Further, using equation (3.30), we getEmbedded Image(3.33)Finally, by condensing the redundant degrees of freedom associated with internal nodes 4 and 7 in equation (3.33), the reduced equation can be obtained asEmbedded Image(3.34)where Embedded Image is the elemental dynamic stiffness matrix with embedded de-lamination.

Figure 2

(a) Modelling of an embedded de-lamination with base-laminates and sub-laminates. Waveguides 1, 2: base laminates and 3, 4: sub-laminates. (b) Representation of the base-laminates and sub-laminates by spectral elements.

4. Damage detection and de-noising using wavelet analysis

In this paper, a wavelet based damage identification technique using time of arrival information is implemented for damage detection from wave propagation responses in de-laminated composite beams. Here, the responses simulated using the developed WSFE model are used to represent the experimental response for the inverse problem of predicting damage location. The procedure involves measuring the time of arrival of waves reflected from the damage and knowing the distance travelled by the wave from the speed of the wave (from dispersion curves) at the excitation frequency, the position of damage can be established. The wave speeds at required frequencies can be obtained from the dispersion relation which, as mentioned earlier, are directly available from WSFE formulation.

However, the implementation of the above technique may not be so direct in many cases. In problems, where the length of the waveguides are very small and if the damage is present very near to the discontinuities, there many be several reflected waves in a very short time interval. In such cases, the different reflected waves are not very differentiable and their time of arrivals cannot be measured directly from such responses. This is illustrated with several examples later in the section on numerical experiments. Wavelet analysis is used here primarily to derive, from such responses, the time of arrival of waves generated by the defect.

Wavelet analysis or filtering gives the time–frequency map of a signal. In other words, the wavelet filtering decomposes a signal into different frequency bands or wavelet levels and within each the time information is retained with a certain resolution. Here, first a fast Daubechies wavelet transform (Mallat 1988) of the response is performed. The wavelet coefficients of the level containing the frequencies of the loading frequencies are retained. In cases where more than one of such levels exists, the lower frequency level is considered for further analysis. These retained wavelet coefficients give a very good estimation of the arrival time of different waves, as has been shown through numerical experiments in §5. Another advantage of this method is that since the higher frequency components are filtered out and not considered in the analysis, the presence of noise, which has higher frequencies, in the response does not affect the results.

It should also be mentioned that, though here wavelet analysis is done for the time domain response, for model based damage estimation methods, the wavelet analysis can be incorporated in the WSFE solver prior to the inverse wavelet transform and will result in large computational savings. Such model based damage detection techniques may be required for more precise estimation of both extent and position.

Experimentally measured responses are always associated with noise and thus de-noising is an important issue. Here, a wavelet based de-noising is done to the simulated responses with added white noise. Since the de-noising method is based on wavelet analysis it can be done along with the damage detection procedure. The universal threshold λ used isEmbedded Image(4.1)where σ is the standard deviation and n is the length of the signal. The soft thresholding (Donoho 1995) to wavelet coefficients di at certain levels is done as follows:Embedded Image(4.2)Similar to the wavelet based damage detection, for simulation of responses from noisy experimentally obtained loading or excitations, de-noising can be done after forward wavelet transform. Thus, this wavelet based de-noising can also be included in the WSFE solver without much increase in the computational cost.

5. Numerical experiments

Here, first the formulated WSFE model of a de-laminated composite beam is used to simulate responses for different fixed-free de-laminated AS4/3501-6 graphite–epoxy beams. The beam configuration is shown in figure 3, where L and Ld is the length of the beam and de-lamination, respectively, and L1 is the distance of the de-lamination from the free end of the beam.

Figure 3

Fixed-free beam configuration with mid-plane de-lamination.

The cross-sectional dimension is 0.01×0.01 m2 with depth 2h=0.01 m and width 2b=0.01 m. Numerical examples are presented for different values of these lengths, ply-lay-up sequences and positions of de-laminations along the thickness of the beam. In addition, the wave propagation responses are studied for both broad-banded impulse and narrow-banded modulated pulse loadings. The unit impulse load shown in figure 4a has a duration of 50 μs and a frequency content of 44 kHz. Similarly, the modulated pulse loading with central frequency of 70 kHz is shown in time and frequency domains in figure 4b.

Figure 4

(a) Broad-band impulse load and (b) narrow-banded pulse modulated at 70 kHz in time and frequency (inset) domain.

The WSFE model is formulated with N=22 and a time sampling rate Δt=2 μs. As mentioned earlier, only a single spectral element is used to simulate the responses.

Next, the responses of the de-laminated beams, simulated with two-dimensional FE, are used for the inverse problem of damage detection using wavelet analysis. As discussed in §4, the position of de-lamination is estimated from the arrival time of waves reflected from de-laminations. Numerical experiments are performed to predict positions of damage from simulated responses, where these responses contain artificially added noise, which are de-noised using wavelet filtering.

(a) Responses of de-laminated beams

First, the WSFE model of a de-laminated beam is validated with responses simulated using two-dimensional FE. The transverse tip velocities of the beam shown in figure 3 with ply-lay-up [0]8, L=0.5 m and L1=0.25 m are plotted for centreline de-lamination length of Embedded Image and compared with the response obtained using two-dimensional FE in figure 5a. The FE results are obtained using 400, four-noded quadrilateral plane stress elements and Newmark's time integration with time step 1 μs. It can be seen that results compare well. However, the small difference in the wave speeds predicted by the two methods can be further reduced by refining the FE mesh.

Figure 5

Transverse tip velocity of fixed-free graphite–epoxy beam due to tip impulse load applied in transverse direction with (a) de-lamination length Ld=20 mm (validated with two-dimensional FE), (b) different de-lamination lengths Ld=10, 20 and 30 mm, (c) de-laminations at different heights above centreline and (d) different ply-lay-up [0]8, [45]8, [60]8.

In figure 5b, the transverse tip velocities of the same beam configuration are plotted, but for different lengths, Ld, of de-laminations and are compared with the undamaged response. The de-laminations are along the centreline of the beam and are of lengths Ld=10, 20 and 30 mm. The excitation is the unit impulse load (see figure 4a) applied at the tip in the transverse direction. It can be seen that in addition to the reflection from the fixed end, the damaged responses show early reflections generated from the de-laminations and amplitudes of these reflected waves increases with increase in de-lamination lengths, as expected. Similar tip transverse velocities are presented in figure 5c, except that, here, the de-lamination length is kept fixed at Ld=20 mm while the position along the thickness directions are varied from h1=h, h/2 and h/4 (see figure 3). As in the previous plot, even here, the damaged responses show reflection from the de-lamination and it can be observed that their amplitudes increases as h1, i.e. depth of the de-lamination from the top surface of the beam decreases. In figure 5d, the transverse velocities due to tip transverse impulse load are plotted for beams with different ply orientation sequences. In all the cases, the beam configuration is similar to figure 3 with L=0.5 m, L1=0.25 m and centreline de-lamination of Ld=20 mm. The three ply-lay-ups used are [0]8, [45]8 and [60]8. Different ply-lay-up changes the stiffness of the beam and hence the wave speed also changes as can be seen from figure 5d, where the responses show different amplitudes and time of arrival of reflections. The [0]8 beam has the lowest amplitude and the time of arrival as it has the highest flexural stiffness and hence highest group speed.

Next, numerical experiments are performed using narrow-banded sinusoidal load (see figure 4b) with central frequency of 50 kHz as input excitation. For such loadings, the waves traverse non-dispersively and are widely used for damage detection applications. The load is again applied in a transverse direction at the tip of a [0]8 beam shown in figure 3 with L=0.5 m and L1=0.25 m. The responses studied are the transverse velocities measured at the tip. In figure 6a–c, the velocities of undamaged and de-laminated beams with Ld=10 and 20 mm are plotted, respectively. In either case, the de-lamination is along the centreline of the beam. Similar to the responses due to impulse loading in previous examples, the damaged responses here show an additional reflection from the de-lamination and their amplitude increases with increase in the de-lamination length. From these plots, the position of damage can be obtained directly using the time of arrival and wave speed.

Figure 6

Transverse velocities of fixed-free graphite–epoxy beam due to narrow-banded load at 50 kHz applied in a transverse direction (a) undamaged, (b) de-laminated Ld=10 mm and (c) Ld=20 mm.

The wave speeds in terms of group speed at a required frequency can be obtained from the dispersion plot shown in figure 7 for the given [0]8 graphite–epoxy beam. The method of extraction of such plots directly from WSFE formulation is given in Mitra & Gopalakrishnan (in press). In this figure, the normalized group speeds Cg/C0 are presented, where C0 is the axial wave speed given by Embedded Image, E, A and ρ being the Young's modulus, cross-sectional area and mass density, respectively. In addition to the transverse bending mode, the dispersion plot shows the shear mode, which propagates only after the cut-off frequency of approximately 105 kHz. The flexural wave speed at 50 kHz is approximately 2.08×103 m s−1. In figure 6b,c, the time of arrival of wave reflected from de-lamination is nearly 0.263 ms. The position of the de-lamination predicted from using the time of arrival and wave speeds is 0.273 m from the tip, which is very near to the actual distance of 0.25 m.

Figure 7

Dispersion relation of graphite–epoxy [0]8 beam.

In figure 8a, the response shown in figure 6b is plotted up to 1.0 ms and in figure 8b, the corresponding noisy signal with signal to noise ratio, SNR=30 obtained by artificially adding white noise is presented. It can be seen that, in the noisy signal, the waves reflected from de-lamination are not distinguishable visually. The de-noised signal is shown in figure 8c and is obtained through universal soft thresholding as described earlier. The thresholding is done only in the first three wavelet levels corresponding to higher frequencies and the order of Daubechies basis used for wavelet decomposition is N=22. Though the first reflected wave from the de-lamination is slightly distorted in the de-noised response, the second reflection from de-lamination can be clearly identified.

Figure 8

Transverse velocities of fixed-free graphite–epoxy de-laminated (Ld=10 mm) beam due to narrow-banded load at 50 kHz applied in a transverse direction (a) simulated response, (b) noisy response and (c) de-noised response.

Figure 9a–c shows the tip transverse velocities due to different positions of de-lamination with Ld=20 mm along the centreline, i.e. L1=0.2, 0.3 and 0.4 m, respectively. The beam configuration is otherwise similar to that used in previous examples. The loading is narrow-banded pulse (see figure 4b) modulated at 50 kHz and is applied at the tip in a transverse direction. Here, the distance of damages from the tip predicted using time of arrival and wave speed are 0.23, 0.33 and 0.44 m, which match well with the actual positions (L1+Ld/2) 0.21, 0.31 and 0.41 m, respectively.

Figure 9

Transverse velocities of fixed-free graphite–epoxy de-laminated (Ld=20 mm) beam due to narrow-banded load at 50 kHz applied in a transverse direction (a) L1=20 mm, (b) L1=30 mm and (c) L1=40 mm (see figure 3).

It should be mentioned here, that, though in the above examples with narrow-banded modulated excitation, the position of the de-lamination can be obtained directly by measuring the time of arrivals from the plots, it may not be so straightforward in many other cases. Particularly, in problems dealing with short length structures or where the damages are very near to boundaries/discontinuities, the incident waves, waves reflected from damages and boundaries/discontinuities are cluttered together. In such situations, it is very difficult to extract the time of arrivals directly from the plots. Wavelet analysis of measured response helps to alleviate this problem to a great extent. This is explained with numerical examples in §5b.

(b) Damage detection using wavelet analysis

A wavelet based damage detection technique as discussed earlier is adopted here for predicting the positions of de-laminations. First, responses are simulated for a [0]8 graphite–epoxy beam with embedded de-lamination of Ld=20 mm along the centreline. Simulations are done using two-dimensional FE, with 400, four-noded quadrilateral plane stress elements and Newmark's time integration with time step 1 μs. The beam is similar to that shown in figure 3 and used in previous experiments except that the length L is shorter and is equal to 0.15 m. The excitation is the narrow-banded pulse with central frequency 70 kHz (see figure 4b) applied at the tip in a transverse direction. In figure 10a,b, the transverse tip velocities of the undamaged and de-laminated beam with de-lamination at L1=0.04 m from the tip are presented, respectively. It can be seen that in the damaged response shown in figure 10b, the incident wave and wave reflected from de-lamination are not distinguishable. Similarly, in figure 10c, the response of the de-laminated beam with de-lamination at L1=0.09 m from tip, i.e. near to the fixed end of the beam is presented. Again in this case, the waves reflected from the de-lamination and the fixed boundary are cluttered together.

Figure 10

Transverse velocities of fixed-free graphite–epoxy beam due to narrow-banded load at 70 kHz applied in a transverse direction (a) undamaged, (b) de-laminated Ld=20 mm and L1=40 and (c) de-laminated Ld=20 mm and L1=90.

To determine the positions of damages from these simulated responses, first wavelet analyses are done using Daubechies basis of order N=22. This generates wavelet coefficients of the responses at different frequency bands or wavelet levels. Here, the wavelet level corresponding to the highest frequency band is referred to as the level one. Since the frequency content of the applied load with central frequency of 70 kHz is 20–120 kHz, the responses will also have the same frequency band. Thus, wavelet decomposition is done only up to first three levels corresponding to frequencies, 125–250, 62.5–125 and 31.25–62.5 kHz, respectively. In figure 11a, the normalized absolute values of the wavelet coefficients at level three are plotted for the undamaged and damaged responses shown in figure 10a,b. The plots show few peaks, which represent different waves. The wavelet coefficients of the damaged response show two additional peaks due to the waves reflected from de-lamination. The time of arrival of the incident wave and the wave reflected from the fixed end are ti=48 μs, tr=192 μs, respectively, for both undamaged and damaged responses and they match exactly. The two additional peaks in the damaged responses arrive at td1=96 μs and td2=240 μs. The speed of the flexural mode at 50 kHz derived from the dispersion relation presented in figure 7 as mentioned before is 2.08×103 m s−1. This speed is used to calculate the position of de-lamination from the arrival time since 50 kHz is approximately the average of the frequency band in wavelet level three. Using this speed, the distance of the fixed end predicted from ti and tr is 0.15 m which is exactly equal to that used in simulations of the responses. Similarly, the position of the damage obtained using both ti, td1 and ti, td2 is 0.05 m from the tip which matches exactly with actual distance of Embedded Image.

Figure 11

Normalized wavelet coefficients of level corresponding to frequency content 31.25–62.5 kHz of undamaged and damaged responses, (a) figure 10a,b, respectively, and (b) figure 10a,c, respectively.

Similar wavelet coefficients for undamaged and damaged responses shown in figure 10a,c, respectively, are presented in figure 11b. Even in this case, the arrival time of the incident wave and the wave reflected from the fixed end are ti=48 μs and tr=192 μs, respectively, which are exactly the same as those obtained in the previous example. The time of arrival of the wave resulting from de-lamination is td1=144 μs and thus, using the wave speed at 50 kHz, the predicted position of the damage from the tip is 0.1 m, which is the actual distance used in simulation.

It should be mentioned here that similar analysis could have been done using the coefficients corresponding to the wavelet level two with frequency content 62.5–125 kHz. However, the level with lower frequencies is used for analysis with the motive of eliminating effects of noises which are generally of higher frequency content and are very often present in experimentally measured responses.

Next, similar wavelet analysis is done for the prediction of damage location from noisy response. The noisy response presented in figure 12a is obtained by adding white noise to the response shown in figure 10b to produce a SNR=15. In figure 12b, the corresponding de-noised response is plotted. The wavelet based de-noising, described earlier, is done at the first two wavelet levels and the order of Daubechies basis used is N=22. In figure 13, the normalized absolute wavelet coefficients at level three are plotted for the undamaged response, noisy damaged response and the corresponding de-noised response. The arrival times, ti, tr, td1 and td2 in this case is exactly similar to those in the previous example shown in figure 11a for the response without any added noise. In addition, it can be seen that the coefficients for the noisy and de-noised responses overlap. These two observations indicate that in the wavelet-based damage detection technique, the presence of noise in the response does not effect the prediction of damage location. This can be explained as the wavelet coefficients used for analysis have low-frequency content and the higher frequency contents including noises are filtered out.

Figure 12

Transverse velocities of fixed-free graphite–epoxy de-laminated beam due to narrow-banded load at 70 kHz applied in a transverse direction (see figure 10b), (a) noisy response and (b) de-noised response.

Figure 13

Normalized wavelet coefficients of level corresponding to frequency content 31.25–62.5 kHz of undamaged, noisy and de-noised damaged responses (figures 10a and 12a,b, respectively).

6. Conclusions

In this work, a WSFE is formulated for a de-laminated composite beam. Spectral element method is an efficient alternative for wave propagation problems and decreases the computational cost substantially. In addition, unlike FSFE, the wavelet based method can be used for accurate modelling of finite length undamped structures. This is primarily because the use of Daubechies scaling functions as a basis allows imposition of initial/boundary conditions in the finite domain analysis. The developed WSFE model of a de-laminated beam is used for simulation of wave propagation in beams with different de-lamination lengths and positions. The responses are studied for both broad-banded and narrow-banded loadings. These numerical experiments performed help in understanding the effects of de-lamination on the wave propagation characteristics.

The simulated responses including responses with added white noise are then used for the problem of damage identification using the wavelet based technique. Particular importance is given to simulated responses where the different incident and reflected waves are not visually distinguishable and thus the damage location cannot be directly predicted from such responses. The positions of the damages detected from these responses using the wavelet analysis, match very well with the their actual positions. In addition, since both de-noising and damage detection involves wavelet analysis, they can be done simultaneously and thus result in computational efficiency. It is also shown that as in this wavelet based damage detection method, the higher frequency components present due to noise are automatically filtered out and their existence does not effect the prediction of damage positions.


    • Received August 22, 2005.
    • Accepted December 21, 2005.


View Abstract