Structural damage detection using energy flow models

Abstract. The presence of a crack in a structure modifies the energy dissipation pattern. As a consequence, damaged structures can present high localized damping. Experimental tests have revealed that crack nucleation and growth increase structural damping which makes this phenomenon useful as a damage locator. This paper examines the energy flow patterns caused by localized damping in rods, beams and plates using the Energy Finite Element Method (EFEM), the Spectral Element Method (SEM) and the Energy Spectral Element Method (ESEM) in order to detect and locate damage. The analyses are performed at high frequencies, where any localized structural change has a strong influence in the structural response. Simulated results for damage detection in rods, beams, and their couplings calculated by each method and using the element loss factor variation to model the damage, are presented and compared. Results for a simple thin plate calculated with EFEM are also discussed.


Introduction
In the last few decades, many works were performed to develop non-destructive inspection methods based on vibration theory, which allow fault identifying and evaluating from modal parameters and dynamic response signal variations [1].In general, these techniques are well suited to detect large defects rather than small defects like a crack, which do not impose appreciable changes in the lower frequencies and the global structure behavior is unaffected.More recent damage detection methods are looking for variations in elastic wave propagation in structures at medium and high frequencies [2].They utilize the evidence that material discontinuities, such as a crack, generate changes in elastic waves propagating in the structure.Experimental tests have revealed that crack nucleation and growth increase damping in a vibrating system.The damping rate increases with increasing crack depth and in some cases its value can be 100% greater than the damping of the uncracked structure [3].
The Finite Element Method (FEM) and the Boundary Element Method (BEM) have an inherent characteristic that produces too large models at high frequencies.Therefore, a vast amount of research activity has been done for high frequency dynamics and one of the most commonly used techniques is Statistical Energy Analysis (SEA) [4].In SEA, the structure is divided into a set of subsystems that interact through energy exchange.However, the spatial variation of the response in each subsystem cannot be obtained as it provides one energy level for each subsystem.In this work, three other methods are used to explore the energy density and energy flow in the field of damage detection in structures at medium and high frequencies.All of them provide information on the spatial energy distribution within the various subsystems.The analytical approximated energy governing equations for rods and beams were developed by Wohlever [5], using the proportionality relationship between energy flow (intensity) and energy density, and an energy balance on a differential control volume.This method is known as Energy Flow Analysis (EFA).Another solution is to solve the EFA equations with standard FE approximations, which is called Energy Finite Element Method (EFEM) [6].EFA and EFEM were also applied to membranes, plates and acoustic cavities.Developed by Doyle [7], the Spectral Element Method (SEM) is the analytical solution of wave partial differential equations in the frequency domain tailored as a finite element method.The Energy Spectral Element Method (ESEM) is a new approach to solve the analytical energy equation of EFA using the same concepts of SEM [8].
In this paper a new approach to detect and localize damage is proposed, which uses the concept of changes in parameters related to the elastic waves propagating in the structure.The changes are indicated by differences in the energy flow pattern between the undamaged and the damage structures.The method utilizes energy flow of longitudinal and bending waves applied to rod, beam, and couplings calculated using the EFEM, SEM and ESEM methods.It was also verified for a single plate structure using the energy flow of bending waves calculated by EFEM.

Spectral element method
The examples treated in this paper were formulated using the elementary rod and Bernoulli-Euler beam elements.The dynamic stiffness relation for the rod spectral element is established via dynamic shape functions between element nodes, but instead of being simple polynomials, they are the exact displacement distributions of the used theory.A complete SEM formulation for elementary rod, beam and plate can be found in [7].The longitudinal displacement for a simple rod is given by: where, k is the wave number, L is the element length, and A and B are constants determined from the boundary conditions.The displacement end conditions for a two-node element are: where û1 and û2 are the nodal displacements.By solving for A and B in terms of nodal displacements, the longitudinal displacement at an arbitrary point in a finite rod is: where, ĝ1 (x) = e −ikx − e −ik(2L−x) 1 − e −i2kL and ĝ2 (x) = −e −ik(L+x) + e −ik(L−x) 1 − e −i2kL are the frequency-dependent rod shape functions.The element force at an arbitrary point in a finite rod can be obtained as where, is the first space derivative, and EA is the axial stiffness.The dynamic stiffness for a two nodes spectral element can be written in a matrix form as: For harmonic excitation the time-averaged energy density for longitudinal waves in a rod can be written as a summation of the potential and kinetic energy densities for longitudinal waves, where <> and * represent time-average quantity and complex conjugate, respectively.The time-averaged energy flow for longitudinal waves in a rod can be written as: where, R is the real part of a complex number.

Energy finite element method
From the EFA, the general governing energy equations for structural wave propagation in a steady state harmonically excited and with small structural loss factor (η << 1) is given by [5,6]: where, η is the structural loss factor, ω is the circular frequency, c g is the group speed, and Π is the input power.The time-average energy flow (intensity) is related to the time-average energy density by: When the Eq. ( 8) is solved as a FE approximation it is called EFEM [6].The formulation is illustrated here in a general form that can be applied to 1 or 2-D energy equations.The weak form variational statement of Eq. ( 8) can be written as or where ν is the trial function.For a Galerkin weighted residual using quadratic elements, the following approximation is made: By considering that ν = Φ j and substituting Eq. ( 12) into Eq.( 11) it produces the element approximation: where, Therefore, the element matrix form for the ith element can be written as: Waves propagating in built-up structures eventually encounter discontinuities, which comes from changes in the structural properties (e.g., changes in material, cross section, and geometry).In these discontinuities there are coupling relationships that take into account for the relative energy flow levels on either side of the joint.To introduce these discontinuities a joint element need to be inserted to connect coupled structural elements and Eq. ( 14) becomes There is one joint element for each type of structural configuration and the number of wave types propagating in the structure.They can be obtained for any coupled structural elements if the energy transmission between them can be described in terms of power transmission and reflection coefficients.

Energy spectral element method
Quoting [7], "The only way to efficiently handle wave propagation problems in structures with complicated boundaries and discontinuities is to develop a matrix methodology for use on a computer".Based on this statement, a new approach was proposed, which is obtained by applying the spectral formulation concept to the energy solution.This is named Energy Spectral Element Method (ESEM), and it allows to solve analytically the Eq. ( 8), using an element assembling scheme.The steady state behavior of the spatial and frequency-averaged longitudinal energy per unit length is given by The general solution of Eq. ( 16) can be written as where, L is the element length, H and G are arbitrary constants determined by the boundary conditions.The energy end conditions for a two-node element are where, < e 1 > and < e 2 > are the nodal energy densities.By solving for H and G in terms of nodal energy densities, the longitudinal energy density at an arbitrary point in a finite rod is: where, p 1 (x) = e ηk(2L−x) − e ηkx e 2ηkL − 1 and p 2 (x) = −e ηk(L−x) + e ηk(L+x) e 2ηkL − 1 are the frequency-dependent rod shape functions.The element energy flow at an arbitrary point in a finite rod can be obtained as The dynamic stiffness for a two nodes energy spectral element can be written in a matrix form as:

Damage indicator
The changes in dissipation properties that are due to local modifications of structural characteristics will produce variations in the energy density of structures.Since energy flow is linearly related with the space variation of the energy density, it may be used as an indicator for detecting and locating damage in structures.Therefore the damage indicator is defined as the difference of energy flow between the undamaged and damaged structure, or where q U and q D , are time-averaged energy flow of undamaged and damaged structure, respectively.For the damage indicator plots the δ values are scaled between 0 and 1.

Beam examples
The structural beam model consists of a traction and flexure element combination, excited at the left end by longitudinal and flexural harmonic forces (Fig. 1) with the following undamaged properties: Aluminum, F = 1.0 N, L = 6.0 m, A = 4 × 10 −4 m 2 , ρ = 2700 kg/m 3 , η = 0.01 and E = 71 GPa.The structure is discretized with 12 elements (13 nodes) and damage imposed at element 3 (nodes 3 and 4) and element 9 (nodes 9 and 10).The frequency-averaged longitudinal and flexural energy flows were calculated by ESEM, EFEM and SEM in a 1/12-octave band with center frequency of 12.5 kHz.  Figure 2 shows longitudinal energy flows (< q L >) and correspondent damage indicators (δ L ) calculated for damage with 1% change in the loss factor.It can be seen that damage is detected and localized by all methods in damage indicator plots (Fig. 2b), but it is not in energy flows plots (Fig. 2a). Figure 3 shows the same plots for damage with 50% change in the loss factor.In this case damage is detected and its position is correctly indicated by all methods in both energy flow and damage indicator plots.In general, previous results demonstrate that damage indicator and energy flow plots can be used successfully to detect and to locate damage.However, there are limits in using energy flow plots to detect the damage if the loss factor changes are low.Due to the characteristic oscillations in the flexural energy density calculated from the SEM solution, it is not possible to obtain a fair result for the flexural energy flow.Figure 4 shows an example with damage of 50% change in loss factor where the flexural energy flow plot is distorted and the flexural damage indicator plot detects the damage, but with false indications for the damage position.Figures 5 and 6 show a comparison of flexural energy flows (< q F >) calculated with ESEM and EFEM and the correspondent damage indicators (δ F ) with damages of 1% and 50% changes in the loss factor, respectively.Both methods present positive damage detection and location in the damage indicator plots, but in the energy flow  plots only for the 50% damage case it shows positive identification for detection and localization.

Coupled beam examples
Two same length (L A = L B = 3.0 m) coupled beams structure were evaluated.In the first structure the beams have collinear coupling, with cross section area ratio of four (A B = 4A A and A A = 4 × 10 −4 m 2 ), and are excited at the left and right ends, each one at a time, by longitudinal harmonic forces (Fig. 7a).In the second structure the beams are connected at an angle α = 30 • , with same cross section area (A B = A A = 4 × 10 −4 m 2 ), and are excited at the left end by longitudinal and flexural harmonic forces (Fig. 7b).For both examples all others properties: discretization, frequency range, damage positions, and damage level are the same used to the examples of Section 3.1.
Figure 8a shows longitudinal damage indicator by ESEM, EFEM and SEM, to the case of different cross section coupled beams.The structure was excited on the left end and a loss factor change of 1% was applied to simulate the  damage.Figure 8b shows similar results to the different cross section coupled beams except that a loss factor change of 50% was applied to simulate the damage and the excitation point is on the right end of the structure.The results show positive damage detection and localization by all methods.Also it can be seen an expected effect of closeness between input power and damage positions, where the damage indicator becomes more evident (higher peaks) at the damage closest to the input power.For the angle coupled beams with excitation on the left end, longitudinal and flexural damage indicators computed with ESEM, EFEM and SEM for a damage level corresponding to 1% are shown at Fig. 9. Figure 10 shows similar results for the angle coupled beam with 50% change in the loss factor.As seen before, in both cases damage indicators present positive detection and localization with longitudinal waves by all methods and with flexural waves by ESEM and EFEM, but fails with flexural waves by SEM.

Thin plate example
The structural model is a thin plate excited out-of-plane by flexural harmonic forces (Fig. 11) with the following undamaged properties: Aluminum, L x = L y = 1.0 m, h = 0.01 m, ρ = 2700 kg/m 3 , η = 0.01, E = 71 GPa.The plate is discretized with 20 × 20 elements of 9 nodes, the damage is imposed on elements 24 (center node coordinates, x = 0.175 m, y = 0.075 m) and 158 (center node coordinates, x = 0.875 m, y = 0.375 m), and the input power at the element 305 (center node coordinates, x = 0.225 m, y = 0.775 m).The frequency-averaged flexural energy flow is calculated by EFEM in a one-octave band with center frequency of 8.0 kHz.Figures 12a and  12b show flexural energy flows (< q F >) without and with damages that cause 1% (element 24) and 50% (element 158) change in the loss factor, respectively.It can be seen in this case that, even with a very strong damage of 50% change in the loss factor, the energy flow plots are completely unable to detect and locate the damages.

Conclusion
A new approach for structural damage detection and localization using the energy flow pattern changes caused by localized damping in rods, beams and plates was presented.The proposed procedure was verified using simulated examples calculated by the Energy Finite Element Method (EFEM), the Spectral Element Method (SEM) and the Energy Spectral Element Method (ESEM).The formulation of the methods (SEM, EFEM and ESEM) for rods, beams and thin plates was briefly reviewed.In general, all methods produced a positive detection and localization of the damage, exception made for the flexural energy flow calculated by SEM, which failed in this task.However, it is known that flexural energy density calculated by SEM presents an oscillatory behavior, and with the high frequency and low spatial discretization used here, it is not possible to obtain a fair result for the flexural energy flow.Physical cracks were not modeled it is assumed that element loss factor increases with increasing of damage (crack).This assumption is corroborated by experimental results reported in literature.The presented results show that the proposed damage indicator is very sensitive and efficient for all analyzed cases.It presents greater sensitivity to detect and locate the damage than the energy flow, whose damage detection sensitivity seems to be damage level  and dimension problem dependent.Numerical results clearly indicate that the methods can be used to detect and localize the damage regions, but it must be admitted that this study represents only a simulation test of the concepts.Further experimental investigations are needed and they are currently under way.

Fig. 1 .
Fig. 1.Beam model excited at the left end by longitudinal and flexural harmonic forces.
(a) Energy flow without and with a damage level causing a 1% change in the loss factor.(b) Damage indicator for a damage level causing a 1% change in the loss factor.
(a) Energy flow without and with a damage level causing a 50% change in the loss factor.(b) Damage indicator for a damage level causing a 50% change in the loss factor.
(a) Energy flow without and with a damage level causing a 50% change in the loss factor.(b) Damage indicator for a damage level causing a 50% change in the loss factor.
(a) Energy flow without and with a damage level causing a 1% change in the loss factor.(b) Damage indicator for a damage level causing a 1% change in the loss factor.
(a) Energy flow without and with a damage level causing a 50% change in the loss factor.(b) Damage indicator for a damage level causing a 50% change in the loss factor.
Figures 13a and 13b show the corresponding damage indicator (δ F ) in a 3-D bar graph and a 2-D stream line plot, respectively.It can be seen that the damage are precisely detected and localized in the 3-D bar graph and only damage regions can be indicated by the 2-D stream line plot.
(a) damage 1% and excitation on left end (b) damage 50% and excitation on right end

Fig. 9 .
Fig. 9. Longitudinal and flexural damage indicator by ESEM, SEM and EFEM for beams angle connection coupling with damage level of 1% change in the loss factor.

Fig. 10 .
Fig. 10.Longitudinal and flexural damage indicator by ESEM, SEM and EFEM for beams angle connection coupling with damage level of 50% change in the loss factor.

Fig. 11 .
Fig. 11.Thin plate model with center node coordinates (x, y) for the excitation and damage elements.
(a) Energy flow of the undamaged plate (b) Energy flow of the damaged plate
(a) Damage indicator 3D bar plot (b) Damage indicator 2D stream line plot